Discovery prospects for heavy dark matter in KAGRA

We discuss the discovery prospects for kg-scale dark matter with a Yukawa-like long-range interaction in Kamioka Gravitational Wave Detector (KAGRA). We consider the interaction range to be in the order of kilometers, thus the dark matter may interact with multiple KAGRA mirrors simultaneously. The induced signal strain is in many cases dominated by frequencies below about 130 Hz. Due to the geometry of the detector there is some strong directional sensitivity. It turns out that KAGRA would be able to constrain such kind of dark matter within a few years of operation. With our set of assumptions we expect that KAGRA could detect a few events per year for 10 kg dark matter, with a Yukawa interaction with a range of 10 km and an effective coupling constant around $10^6$.


INTRODUCTION
Despite dominating the matter density of the universe, little is known about the nature of dark matter (DM). The lower bound of the dark matter mass is about 10 −22 eV while the upper bound is around five solar masses [1]. Given that there is an 80 orders of magnitude difference between the lower and upper bound of this mass range, many different possible experiments need to be considered in unraveling the properties of dark matter. Many current direct detection experiments have explored dark matter in the WIMP mass window around 100 GeV and axion dark matter with masses well below 1 eV. These mass regions are motivated by the solutions to other problems. Primordial black holes with masses heavier than asteroid masses have also recently gained attention as heavy dark matter candidates. But all direct searches have not found any conclusive evidence so far. Large regions of model parameter space have been constrained (see Refs. [1][2][3][4]), but there are still an enormous mass and parameter range left which have been weakly explored, if at all. Hence, new ways to look for dark matter focusing on other mass regions and detection techniques have recently gained more attention and we will also go in this direction here.
To be precise, we want to study here the possibility of using a gravitational wave detector as a DM detector. There have been proposals going in a similar direction using either gravitational wave detectors or generally speaking optomechanical experiments . However, the mass range we consider is very different from many of these proposals and hence we will use a different formalism.
We consider here one of the possible DM candidates in the kg-to ton-scale mass range. The DM in this mass range interacts significantly with the proposed detector while we can still have a detectable rate. This class of DM can be motivated by several models such as primordial extremal black holes [41], dark quark nuggets [42], Fermi balls [43] and electroweak symmetric DM balls [44]. Currently, there is no dedicated direct detection search for this class of the dark matter to our knowledge. There have been some proposed searches sensitive to such DM using an array of force sensors [45,46] or precision displacement sensors [5]. It is also possible that dark matter can have an additional Yukawa-like interaction with the normal matter besides the gravitational interaction. What will be of particular importance here is that for a longer range interaction, the dark matter might interact simultaneously with two or more different objects located in some distance. Therefore correlations between the different objects need to be taken into account. This can lead to some interesting phenomena, such as signal enhancement or suppression and directional sensitivity. This long-range interaction also implies that dark matter can have very large cross sections for interactions with ordinary matter compensating the very low number density.
Our work here closely follows the model and general proposal in [6] which discussed the prospects for looking for heavy DM in LIGO. But we will go significantly beyond that adding details on the modeling of the signal-to-noise (SNR) ratio, including a more realistic noise model, and discussing more explicitly various parameter depend-encies. The work here relies heavily also on the insights from [35] by some of us where we modeled DM scattering in KAGRA. The difference is that in this previous work we considered light DM with a very short-range interaction. In our work here instead, we will consider a kg-to ton-scale DM with O(km)-range interaction with the normal matter. This suggests that we look for a small, but correlated force on objects which are also about O(km) away from each other. Ideal detectors for that are terrestrial gravitational wave detectors which can detect minuscule relative displacements of their mirrors arranged in an interferometer. We consider here as an explicit example the currently running KAGRA observatory. KAGRA is the most recent, large scale gravitational wave detector with the most modern technology making it an ideal benchmark for the near future as well. We use the realistic KAGRA suspension setup and noise estimated by the KAGRA collaboration in this work to provide a reasonable benchmark.
Our paper is organized as follows: In Sec. 2, we explain how to calculate the observable for a DM signal in KAGRA and similar experiments, starting from the potential which consists of a Yukawa and a Newtonian term. With this information and having a noise model at hand we discuss how to calculate the SNR. In Sec. 3, we then show the dependence of the SNR on the DM phase space and other dark sector parameters. In Sec. 4, we estimate the event rate per year for various DM masses, the Yukawa range, and the coupling strengths. We summarize and conclude this work in Sec. 5.

MODELING HEAVY DM SCATTERING IN KAGRA
In this paper we assume that DM has a gravitational and a long-range Yukawa interaction with the mirrors in KAGRA. The potential for this interaction with one of the test masses is, cf. [6], where x T , x DM and m T , m DM are position and mass of the target and the passing DM object respectively. The constants δ SM and δ DM are effective coupling constants. G N is Newton's constant and λ parametrizes the range of the Yukawa interaction. The factor s = +1 for a scalar and tensor exchange and −1 for vector exchange. Throughout the paper the Yukawa part will be dominant and we set s = +1. The differences to s = −1 are marginal in all numerical results in this paper. We assume where r DM and r T measure the size of the DM object and the target, respectively, and L is the length of the interferometer arms. While d DM is defined as d −3 DM ≡ n DM in which n DM is the local DM number density. This implies that we can approximate the DM object and the targets as pointlike objects.
In this work we will focus on DM masses between 0.1 and 1000 kg, ranges of Yukawa forces between 0.1 and 10 km, and an effective SM-DM coupling |δ SM δ DM | between 10 3 and 10 7 . One could be concerned at this point with constraints from fifth force searches. Indeed, if this would be an SM-SM force, that would be ruled out. There the Yukawa coupling |δ SM | 2 is constrained to be less than about 5×10 −4 [47] for the considered range in λ.
But as it was also discussed, for instance, in [6] and [37], this bound does not directly apply here but it can be combined with a limit on the cross section of DM self-interactions, σ T /m DM 1 cm 2 /g at v DM ∼ 1000 km/s. In our notation and |δ DM | 1.7 × 10 10 kg 1/4 m where we have used log (λ/r DM ) = 5. Hence, after combining the fifth force and DM selfinteractions constraint |δ DM δ SM | 10 8 and our considered range |δ SM δ DM | between 10 3 and 10 7 is below current bounds.
Since the realistic forces are tiny, we neglect the back-reacting force on the DM object and hence its velocity v DM is constant, where x 0,DM = x DM at t = 0 is the initial position of the DM particle. Of course, what time we choose to be t = 0 is arbitrary and we will choose it such that | x 0,DM | is minimal in the frame with the corner of the interferometer at the origin. If there would be only one target mass at the origin | x 0,DM | would be very close to the impact parameter of the scattering process. Please note that in the following we assume the DM velocity to be defined in the lab frame and write v DM = v lab DM . This will become important when we discuss the assumed DM velocity distribution which is usually defined in the galactic rest frame.
KAGRA is an L-shaped interferometer with four mirror suspension systems located at both ends of each interferometer arm, cf. Fig. 1 of [48]. We choose one of the arms to be along the x-direction of our lab frame and the other one to be along the y-direction. Their respective positions are labeled as EX, IX, EY, and IY [external (E)/internal (I) and x-/ydirection (X/Y)]. Each suspension system can be described as a system of three coupled, damped, harmonic oscillators, which are called intermediate mass (IM), blade springs (BS), and test mass or target mirror (TM) [35,49]. As shown in [50] and based on the experimental setup information given in [49] in each system the TM is at the bottom, the BS are on top 30 cm above the TM and the IM is 2 cm below the BS. From now on, we will use indices i = 1, 2, 3, to refer to IM, BS, and TM respectively.
Let us first consider the inner suspension system in the arm pointing in x-direction (IX). Hence, we consider here the horizontal external force F ext,h to point in x-direction as well. The equations of motion for this suspension system are given as (2.8) Note that the three components stand here for the components of the suspension system, not the three directions of the lab frame. The displacementsx ih are defined as the deviations in the lab frame x-direction from the respective equilibrium positions depending on the position of the mounting components. That means the total displacement of the IM in the lab frame in x-direction isx 1h , the total displacement of the BS isx 1h +x 2h and for the TM it is x 1h +x 2h +x 3h . The spring constants K ih ≡ k ih (1 + i φ ih ) are complex, for their values please see Ref. [49]. The masses of the suspension system components are m 1 = 20.5 kg, m 2 = 0.22 kg and m 3 = 22.8 kg.
One also has to consider the forces in vertical direction (the z-direction of our lab frame). The form of the equation of motion for vertical direction would be the same but the values for the spring constants are different. For more details, see Refs. [35,49].
In both cases, we consider as external force the long-range force induced by a passing DM object, i.e.,F ext,ih (t) = −∂V DM,i /∂x ih for the horizontal x-direction andF ext,iv (t) = −∂V DM,i /∂x iv for the vertical z-direction, where we have to use the respective positions of the suspension components in the potential.
We can take the Fourier transform of the equations of motion above and find where we label Fourier transforms of quantities with a tilde and we have also absorbed a factor of L to make the displacement in frequency-domain dimensionless.
Since the coordinates x ih in Eq. (2.6) are defined relative to its predecessor, the Fourier transform of the TM x-position is given bỹ The other directions are completely analogous. The so-called horizontal direction points in the y-direction of the chosen frame for the other interferometer arm. It is estimated that the vertical modes cross-couple to the horizontal modes due the curvature of the earth and mechanical imperfections [51]. Hence, the relevant quantity from a single suspension system is where VHC = 1/200 is the vertical-horizontal coupling. The strain amplitude for just one dynamical suspension system is then given by In Fig. 1, we plot it as a red dashed line, where we chose for illustration, for the dark sector parameters, m DM = 10 kg, x 0,DM = (0, 2 km, 2 km), v DM = (220 km/s, 0, 0), λ = 1 km, δ SM δ DM = 1 × 10 5 and s = +1. The effective coupling constant value chosen is within the limit stated in Eq. (5) from [6]. In this figure we also show the estimated KAGRA noise in the detuned resonant side-band extraction configuration taken from [49] for comparison. We see that the DM signal picks up the characteristic resonance frequencies like the suspension thermal noise. We have seen this before also for an instantaneous, short-range force; see Ref. [35] for detailed explanations.
Since we consider here a long-range force which will act on all four suspension systems simultaneously we have to extend this formalism to that case. As mentioned before in KAGRA the four suspension systems are arranged in an L-shaped interferometer and labeled as EX, IX, EY, and IY. The rest positions of the four TM mirrors in our lab frame coordinate system are (2.14) Note that the horizontal mode for EX and IX is in x-direction while for EY and IY it is in ydirection. The vertical mode for all four suspension systems is in z-direction. The calculations for each suspension system and mode is in complete analogy to the case discussed above. The observable strain amplitude by a DM mass passing by for the full system is then In Fig. 1 we compare the strain amplitude for the one-suspension setup with the foursuspension setup with the same dark sector parameters. We see that there is a difference between the DM signal with the four-mirror setup and the one-mirror setup. For example, That there is such a significant difference between the one-and four-mirror setups should not come as a surprise. Here we consider here a long-range Yukawa interaction together with Newtonian gravity which implies that there are configurations where there is a significant force on more than one of the mirror strings. That is different to our previous work [35] where we considered only an instantaneous, very short-range force in one of the mirror components. Hence, we need to consider here the quantity in Eq. (2.14) and simulate the force on all mirror components. Please note that the enhancement from one mirror to four mirrors is not trivial. As we will see later, the enhancement is not uniform across DM parameters and cancellations can happen as well. Hence, the detailed studies done in the following sections are necessary.
We define the signal-to-noise (SNR) for the rest of the paper as where we choose f min = 3 Hz and f max = 133 Hz. Here, we need to make a trade-off. Naively, one can get the largest SNR by integrating from zero to infinity which is numerically impossible and we have to choose a finite range. At low frequencies one is quickly dominated by seismic noise and it makes sense to cut off the integration at O(1) Hz. For large frequencies the choice is less obvious. In Fig. 1 we see that the DM signal can bend over at a certain frequency and then quickly go to zero. But at which frequency the signal starts to bend depends strongly on the chosen parameters, cf. Fig. 2 where we chose a number of different strains for some plausible dark sector parameters. At larger frequencies quantum noise dominates apart from some resonances from the thermal noise which we neglect for now. In KAGRA the quantum noise has roughly a V-like shape and we can split it into two frequency ranges. Between about 10 Hz and f int ≈ 93 Hz and for frequencies larger than f int where we neglected the resonance peaks from the thermal noise. In many cases the DM signal on the other hand scales roughly as where we again neglect (anti-)resonances and also ignore that at some frequency the signal can drop much faster. Then (2.20) What we see from this approximation is that the frequency range below f int is much more important. Going to higher frequencies also helps but the contributions are comparatively suppressed. And this is under the assumption that the signal falls slowly even for large frequencies which is often not the case. In fact, in our analysis later we rarely care about the precise value of ρ 2 KAGRA . We are mostly interested in relative and qualitative features and the region where ρ 2 KAGRA > 1. This explains why we fix f max = 133 Hz which is quite a bit larger than f int to be on the safe side while keeping computation time feasible. We will also show that, for instance, the expected observable event rate at KAGRA depends only weakly on f max for that large frequency later on when we discuss our numerical results. Nevertheless, we want to stress that in the following sections we will simulate the full SNR including the resonances and not use the above approximations. They are used only to justify the chosen frequency range for the numerical integration.

DIRECTIONAL AND PARAMETER DEPENDENCE OF THE SNR
In this section we want to have a closer look at the directional and dark sector parameter dependence of the SNR which we just defined. Since KAGRA measures the relative dis- placement of the detectors in two directions there is some apparent intrinsic directionality related to the initial position (the closest point to the mirrors) and velocity of the passing DM particle. It is challenging to discuss the possible six-dimensional DM phase comprehensively. Therefore we will first fix the initial position and just change the direction of the DM velocity and then change the initial position and consider a velocity orthogonal to the initial position vector to give some impression of some of the important features.
For the first case, we set x 0,DM to be the origin of the coordinate system which is at the intersection of the two interferometer arms. We also fix the absolute DM velocity to be v 0 ≡ | v DM | = 220 km/s and then only vary its direction. Hence, we are left with two parameters (angles) and we set v DM = v 0 (cos φ cos θ, sin φ cos θ, sin θ) . (3.1) We then scan log 10 ρ 2 KAGRA over the polar angle from 90 • to −90 • and the azimuthal angle from 0 • to 360 • and plot the result in a Mollweide projection in Fig. 3. For the other dark sector parameters we chose m DM = 10 kg, λ = 1 km and δ SM δ DM = 10.
The four brightest spots around θ = 0 • , which are φ = 0 • , 90 • , 180 • and 360 • , correspond to when the DM particle is basically going along the x-and y-direction, that means it is getting very close to the suspension systems and the exerted force is large. The force and hence the SNR is actually getting so large that we took the logarithm of the SNR in the figure and restricted the color scale to keep all relevant features visible.
Apart from the four very bright spots there is another feature. If the DM passes by in the direction for φ ≈ 45 • and 225 • the SNR is suppressed. That is when the DM passes through the right angle spanned by the L. In this case the force exerted on both arms is very similar and the effect cancels out when taking the difference between the two arms, cf. Eqs. (2.14). On the other hand, for the direction orthogonal to that case, the effect adds up when taking the difference and the signal is enhanced. What we can also see is the mirror symmetry about the x-y-plane leading to an intrinsic ambiguity if one would try to reconstruct the DM velocity.
In the second case in Fig. 4, we fix the dark sector parameters as before in Fig. 3 where m DM = 10 kg, λ = 1 km and δ SM δ DM = 10. We also fix the DM velocity to be 220 km/s and the Mollweide projection maps its different directions as before. That means the two angles again have the same meaning. What changes now is the initial position. Here the DM does not pass through the origin but it passes through a ring in the x-y plane with a radius of 3.5 km. The initial position vector is chosen orthogonal to the velocity vector which still leaves a sign ambiguity. This is resolved here such that when v DM was pointing in negative ydirection, i.e., v DM = (0, −| v DM |, 0), the position vector was given by x 0,DM = (+3.5 km, 0, 0). Similarly for v DM = (| v DM |, 0, 0) the position vector was given by x 0,DM = (0, +3.5 km, 0). We have made a sketch in Fig. 5 to make this clearer. That explains why these two directions are brighter (have a larger SNR) in the Mollweide projection in Fig. 4 since here the DM comes closest to one of the suspension systems. Note that the angle φ labels the direction of the DM velocity. And in this figure φ = 0 • means it passes by close to the outer detector of the Y arm and for φ = 270 • it passes by close to the outer detector of the X arm.If the DM velocity points in the opposite direction it also passes by at the other side of the circle and it is much further away from any suspension system. And hence, the effect on KAGRA and the SNR are smaller.
At this point we can also already study how the boundary ρ 2 KAGRA = 1 changes when changing the dark sector parameters. This boundary is important since it will be used in the next section to classify which events are considered observable.
In Fig. 6 we assume the DM velocity to be in z-direction with a value of 220 km/s. The In all three plots the dents in the lines are the cancellations which we have already partially seen before in Fig. 3. Apart from that we can also clearly see the forces on the mirror and hence the SNR scales. As the DM mass, δ SM δ DM and λ get larger, the forces felt by the suspension systems become larger as well. Hence, the area where ρ 2 KAGRA > 1 becomes larger. For example the ratio between the area for m DM = 1000 kg and m DM = 10 kg is around 2.57 while the ratio for m DM = 10 kg and m DM = 0.1 kg is around 4.75. While the mass scales by a factor of 100, the area grows only by an order one factor. So we see that the scaling is nonlinear and that it does not even follow a simple power law. Since the number density of the dark matter is inversely proportional to the mass of the dark matter, this observation will have an impact on the sensitivity as discussed in the next section.
In the cases of changing δ SM δ DM and λ, the area for ρ 2 KAGRA > 1 grows larger without any reduction on the DM density. Hence, we expect that the observable rate will grow larger as any of the parameters increase. In the case of varying δ SM δ DM , we found the area ratio to be A 10 6 /A 10 5 = 1.68 and A 10 5 /A 10 4 = 1.98 where the subscripts denote the values of δ SM δ DM . For the λ dependence there is another remarkable feature. For a small value of λ, the DM interacts significantly with one mirror only as shown in the case of λ = 0.1 km. Hence we found that while the area ratio A 10 /A 1 = 5.17, the ratio A 1 /A 0.1 becomes 50.4. Here the subscripts denote the values of λ in km.
Finally we note that in all the plots in this section we use the dark matter velocity to be 220 km/s, which is a typical value in the standard halo model. For different DM velocities, the shapes and numbers might change in some plots. However the overall trends do not change so we will not discuss these changes in detail here. Besides, we used a much lower effective coupling constant value when plotting the Mollweide projections to make sure that the numerical values of the SNR do not blow up when the DM particles pass through along the beamline. With that we can present a much clearer picture of the SNR variation. Similar is the case for the DM velocity, where if the coupling constant value is changed, the overall outcome will remain the same. In the later sections, we will scan over the coupling constant when we calculate the event rate.

EXPECTED EVENT RATES
To provide potential discovery prospects of KAGRA for the type of DM we consider here that we have to calculate the expected observable event rate, R obs . That is the number of expected events in a unit time with an SNR value greater than a threshold value we will set to one.
To derive this rate we begin with the distribution of SNR values which is given by g(ρ 2 ) d ρ 2 is a normalization factor and v DM and x 0,DM are the threedimensional velocity and initial position of the DM particles, cf., Sec. 2. The functions f v and f x are normalized distribution functions for velocity and position. This equation is difficult to evaluate since it is nontrivial to find the initial conditions for a given value of ρ 2 . In other words it is difficult to find the phase space values where the Dirac δ-function does not vanish.
But since we are only interested in the number of events with an SNR value greater than a certain threshold value, ρ 2 th , we can calculate instead which is much easier to calculate than g(ρ 2 ). The p(ρ 2 th ) is the fraction of DM events which would give an SNR greater than ρ 2 th . Since the total fraction of events with ρ 2 > ρ 2 th and ρ 2 ≤ ρ 2 th must be one we can also easily calculate, the normalization factor We assume that all initial positions are equally likely, that means f x = 1/V where V is the considered volume. Naively one might want to take the solar system as the volume or something of similar size. But this is not necessary since the volume where the DM can have a significant impact on the detector is much smaller and using a smaller volume is numerically more efficient. To estimate the maximum radius of positions we have to consider to get a detectable signal we define F min as the minimal force on the KAGRA mirrors which could be resolved experimentally. From we can then calculate a maximum distance between the DM particle and the KAGRA mirrors where the DM might still be detected. For DM to pass by in a further distance the force would be weaker and a detection would be impossible. We can estimate F min ≈ m 3 g ∆d min l , (4.5) where g is the gravitational acceleration and l = 0.35 m is the length of the wire that connects the blade springs and the test mirror with mass m 3 . The minimum displacement that can be resolved is given by ∆d min = h min L √ ∆f where we use the strain sensitivity of KAGRA h min = 2 × 10 −24 Hz −1/2 at f = 100 Hz [52] and L is the interferometer arm length. The bandwidth ∆f is taken as 10 Hz which is a typical value used for analysis in Gravitational Wave (GW) detectors [53]. With these numbers, we get F min ≈ 1.21 × 10 −17 kg m/s 2 .
Putting this value back into Eq. (4.4), we can then calculate the corresponding maximum distance between the DM particle and one of the KAGRA mirrors for a given set of dark sector parameters which we need to consider in our scans. For instance, for λ = 1 km, δ SM δ DM = 5 × 10 5 and M DM = 10 kg we find r max = 15.8 km. But this is the distance from one of the mirrors. To take into account all suspension systems, and as a simple and conservative choice, we consider as volume V a sphere with radius centered on the origin of our lab frame. We checked numerically that this estimate is quite robust and could not find any parameter points where we should have chosen a larger radius. In fact, in many cases this radius is significantly overestimated but since it is simple to implement we use this volume in our numerical calculation. Apart from that volume restriction we used another observation in our calculations to scan the initial conditions for DM more efficiently. In our analysis we ignore the backreaction of the detector on the DM particle. We assume it follows a straight line with a constant velocity to a good approximation. That implies we get the same SNR value no matter what point on the line we choose as initial position.
The SNR depends just on the absolute value |x 4m (f )|. Choosing different initial points on a DM trajectory would lead to a different global phase ofx 4m (f ) which disappears when taking the absolute value. We can therefore represent each DM trajectory by any point on it and we choose the point closest to the origin of our coordinate system. This point has the advantage that the position vector and the DM velocity are orthogonal to each otherv where we introduce the normalized DM velocity vector,v DM , for later convenience. In fact, this defines a plane, not just a trajectory but our argument remains true for every trajectory on that plane as well. Note that the SNR value is in general different for different directions on that plane due to the directional sensitivity of the detector. The length of any trajectory within the considered sphere is 2 R 2 max − x 2 0,DM . We can then introduce this insights into our integral as a δ-function multiplied with the length of the trajectory where These restrictions do not only limit the coordinate integral to a meaningful volume and get rid of one integration; they have another useful side effect. For the Fourier integral to get x 4m we only have to integrate over the time interval from −T to T where T = R max /| v DM | since only during that time the DM could have a measurable impact on the detector.
What is left to discuss before we can evaluate the integral is the DM velocity distribution. In this work, we employ for it the Standard Halo Model (SHM) given by a Maxwellian distribution in the galactic rest frame. The DM velocity in this frame is labeled as v gal DM . Its distribution has a dispersion σ v and is truncated at the escape velocity v esc [54][55][56]. But we are working here in the lab frame. The relationship between the DM velocity in galactic rest frame and the lab frame is given by where v lab (t) is the time-dependent velocity of the lab in the galactic rest frame. The timedependence here leads to the well-known daily and annually modulation of the expected DM rate, cf. [54][55][56].
The DM velocity distribution in the lab frame is then given as for | v gal DM | < v esc and zero otherwise. The normalization constant is given by (4.12) In our numerical calculations we performed the velocity integration actually over v gal DM which is trivial since we fixed v lab (t) to be a constant value. For notational simplicity, we will remain to use v DM as the DM velocity in lab frame from now on even in integrals. At this point it is important to point out that this assumed DM velocity distribution should be taken with a grain of salt. Not only that this velocity distribution might be incorrect within our astronomical neighborhood, the gravitational and Yukawa backreaction and other new forces might modify this via friction, attenuation, focusing, clustering, etc. This could lead to small deviations but it could even imply that a significant fraction of DM could not penetrate the Earth to reach KAGRA.
Additional new interactions are highly model dependent and a discussion of them is beyond the scope of the current work. Considering only the Yukawa interaction, the velocity distribution can still be modified. In [37] they considered a similar scenario but with a much larger Yukawa range λ ∼ O(AU) and they found deviations induced by interactions with the sun. In our case due to the much shorter range we expect deviations mainly induced by the Earth itself. This can lead to an enhancement of the expected signal rates due to local clustering or a depletion due to shielding and repulsion effects. Since our work here focuses on the modeling of DM interactions in KAGRA we leave the modeling of a more realistic velocity distribution and number density including possible new interactions for a future work. Now we have everything together to numerically calculate p(ρ 2 th ) for any given value of dark sector parameters using the Vegas Monte Carlo routine in Python [57]. As a warmup we want to pick up on a previous issue and study how p(ρ 2 th ) depends on the integration boundaries for the SNR, f min and f max , cf. Eq. (2.16). In Fig. 7 we show the dependence of p(ρ 2 th = 1) for these two boundaries independently leaving the other boundary fixed to our chosen value. We can see that p(ρ 2 th = 1) depends only weakly on f max as long as it is above f int justifying our choice of f max = 133 Hz. We also see that the difference for f min of O(1) is marginal.
We now turn finally to the expected observable event rate where R total is the rate of DM particles entering (or leaving) the sphere given by the radius R max in a given time and ρ 2 th is the threshold value for observability. To calculate it though we need to first calculate the total rate of DM particles passing by. That is the number flux of DM particles Φ DM going through the surface of the sphere, i.e., the integral where the factor 1/2 is to avoid double counting (we only count DM particles entering OR leaving the sphere) and S is the normalized normal vector on the surface of the sphere. Now we evaluate this formula We will first evaluate the integral d 2 S. For this integral we can keep the direction of the DM velocity fixed and the scalar product is just the angle between the two directions, θ S For September 1, 2022, we find numerically R total = (13000 km/s) × n DM R 2 max . For this paper we choose corresponding to a DM energy density of 0.3 GeV/cm 3 . The annual and daily modulation induced by the motion of the Earth relative to the galactic frame would modulate this number by a few percent and we will neglect the modulation for our final results. In Fig. 8, we finally show the expected rate for various values of m DM , λ and δ SM δ DM . The number of points and iterations is chosen such that the Monte Carlo integration error is smaller than the dot size in the figure.
From the plots we can see that lighter DM results in a bigger rate which is driven mostly by the larger number density of DM. For example, in the case of m DM = 0.1 kg, the rate can be larger than once per year for λ 1 km and δ SM δ DM 10 4 . If the DM is heavy enough, e.g., m DM = 1000 kg, the number density is too low to generate enough of a rate. This can be seen from Fig. 6 where we plot the ρ 2 KAGRA = 1 boundary for various DM masses. In the plot, the ratio between the area for m DM = 1000 kg and m DM = 10 kg is around 2.57. However the number density of the heavier dark matter drops by a factor of 100. Therefore the larger volume of the ρ 2 KAGRA = 1 region cannot compensate for the decrease of the number density in the rate calculation.
For the other parameters (δ SM δ DM and λ), the boundary gets larger for larger values without any decrease of the number density. Hence the rate for those cases are higher for a larger value of parameters. For many parts of parameter space, we found that the rate is one event per year or higher. In those cases, KAGRA could be able to see or constrain dark matter given a few years of operation.
In Fig. 8, we see that for a lower value of δ SM δ DM the distance between the points for different λ are smaller compared to larger values of δ SM δ DM . For smaller values to get ρ > 1 the DM has to get very close to the suspension systems and the exponential factor from the Yukawa potential, cf. Eq. (2.1), is practically one. For larger values of δ SM δ DM on the other hand it is possible to get ρ > 1 for distances where the exponential factor becomes relevant. That affects how fast the fraction of significant DM events grows.
We also note here that KAGRA could not observe any events for DM with purely Newtonian interactions, i.e., for δ SM δ DM = 0 in the considered DM mass range.

SUMMARY AND CONCLUSIONS
In this work we discuss how a GW detector can be utilized as an alternative technology to look for a certain kind of DM. Here, we assume DM has a gravitational and Yukawa longrange interaction with ordinary matter exert-ing an external force on the mirror suspension systems in the laser interferometer of a GW detector. The Yukawa force is assumed to be proportional to the DM and the target mass like the gravitational force. It is then possible to get relatively high strain amplitudes for heavy DM which can potentially overcome the background noise in the detector. Another important advantage of that assumption is that the long-range interaction leads to a large cross section which can overcome the small local DM number density and we can still expect reasonable event rates for a large portion of the possible parameter space.
We first model in detail the DM strain amplitude for KAGRA considering all four sets of suspension systems, one at each end of each interferometer arm. We compare it with the outcome from a single-mirror setup which is very similar to the case considered in our previous work [35]. We see from our results that it is crucial to model all four suspension systems as the strains of the one-and four-mirror setup are considerably different depending on the initial conditions. To calculate the SNR from the DM signal strain one should naively integrate over a large frequency range which is numerically challenging. But we could show with a simple estimate and numerical checks that low frequencies especially below about 100 Hz in most cases dominate the SNR, in particular for SNR values around one.
It is obvious from its geometry that KAGRA cannot measure forces in all directions equally well. For that reason we had a closer look at its directional sensitivity for DM and studied how the SNR is affected by different DM velocities and initial positions. The SNR is obviously large when the DM would pass by close to one of the suspension systems. But we could also see that the effect of the force on the two arms can add up or cancel in the observable phase difference. Hence, KAGRA can even act as a directional DM detector although there are some degeneracies between different directions.
We then calculate the expected observable event rate which we define here as the rate of events with an SNR greater than one. We assume DM follows the standard halo model velocity distribution and then find that it is in principle possible for KAGRA to detect a few events per year for 10 kg DM, a Yukawa interaction with a range of 10 km, and an effective coupling constant around 10 6 . For lighter DM one could even expect more events since the number density increases faster than the cross section decreases. Hence, KAGRA and other GW detectors can serve as DM detectors for kg-scale DM with a long-range interaction, a possible scenario which is so far rather poorly constrained.