Theoretical Radio Signals from Radio-band Gravitational Waves Converted from the Neutron Star Magnetic Field

Now, we show the probability of converting GWs into electromagnetic waves in a typical neutron star magnetic field. In this paper, we only consider the simpler, rapidly rotating magnetic field structures. In our subsequent paper, we will simulate the entire neutron star’s magnetic field using the particle-in-cell method and refine the simulated field with current observational data (W. Hong et al. 2025). In the context of rapidly rotating magnetic fields, pulsars and magnetars can often be approximated as rotating magnetic dipoles, with the magnetic dipole axis aligned to the rotation axis, represented as , where r₀ = 10 km is the radius of neutron stars, and the magnetic dipole axis is aligned with the rotation axis where θ represents the polar angle from the rotation and magnetic axis, and B₀ is the surface magnetic field at the poles. While the Goldreich–Julian model (GJ model) was first suggested for aligned neutron stars with θ_m = 0, it works just as well for oblique neutron stars (P. Goldreich & W. H. Julian 1969). It gives a charge density of , where Ω is the neutron star’s spin period, and θ is its polar angle with respect to the axis of rotation. We shall use the charge density as an approximate measure of the electron number density: . Therefore, the plasma frequency is denoted as where is the component of the magnetic field along the perpendicular direction of the traveling direction +x_s of GWs, and depends on time due to the rotation of the neutron stars. As we only consider the cold plasma scenario in order to simplify the model, relativistic effects on the plasma frequency correction are not relevant. Electromagnetic wave propagation in plasma is based on the premise that the frequency of the electromagnetic wave is higher than that of the plasma. This allows us to calculate the minimum spatial scale of the inverse GZ effect occurring radially using the formula

It is important to note that, when the GW frequency is higher, the calculated result of r_occur will be less than the radius of the neutron star r₀, which is not physical. Therefore, when the calculated result of r_occur is less than the neutron star radius r₀, we fix r_occur to r₀ = 10 km. Obviously, when the radius is greater than r_occurs, the GZ effect still exists, but the conversion probability decreases as the neutron star’s magnetic field strength decreases.

GWs are converted into photons in neutron stars’ magnetic fields. However, some neutron stars’ magnetic fields are stronger than the critical magnetic field (D. Lai & E. E. Salpeter 1995; D. Lai 2001), which will lead to some additional physical processes, where ∣q_f∣ denotes the elementary charge. In the presence of an external field, the interaction of observables must be analyzed due to the resummation of higher-order diagrams, leading to nonlinear dependencies known as “nonlinear QED.” Within the frequency range of 10⁶–10¹¹ Hz, the photon’s wavelength is λ_γ ≈ 3 × 10⁻⁴–300 m, which is much larger than the electron’s Compton wavelength: . Additionally, their energy E_γ ≈ 10⁻¹⁰–10⁻⁵ eV, significantly smaller than the electron’s mass m_e = 0.511 MeV, thereby validating the Heisenberg–Euler effective Lagrangian (W. Heisenberg & H. Euler 1936; J. S. Schwinger 1951). Consequently, the action of GW–photon conversion within the proper time integral is (S. L. Adler 1971; W.-y. Tsai 1974; W.-y. Tsai & T. Erber 1974, 1975; D. B. Melrose & R. J. Stoneham 1976; L. F. Urrutia 1978; W. Dittrich & M. Reuter 1985; W. Dittrich & H. Gies 2000)

where R represents the Ricci scalar, and g denotes the determinant of the metric g_μν. The effective Lagrangian is given by L_eff = L⁽⁰⁾ + L⁽¹⁾, where is the original Maxwell Lagrangian, and L⁽¹⁾ is equal to , and represents the effective action for the gauge field A_μ. In , τ represents the proper time,  > 0 is an infinitesimal parameter used for the resum of the effective action, tr signifies the remaining trace over the Dirac spinor space, and denotes the Hamiltonian. D^μ =∂^μ + iq_f A^μ(x), , and indicates the free Hamiltonian (W.-y. Tsai & T. Erber 1974, 1975; K. Hattori & K. Itakura 2013a, 2013b). Under the one-loop effective action assumption and the Bianchi identity with the relation , we can expand the Lagrangian L⁽¹⁾ to the first order in the form with the heat-kernel method , where the ’s denote local functions of the field strength tensor F_μν, and the background gauge is assumed (V. P. Gusynin & I. A. Shovkovy1996, 1999; D. V. Vassilevich 2003). Moreover, the two Lorentz invariants of a mass-dimension of four are and . Then, the remaining trace term in L⁽¹⁾ can be obtained as , where is . In the term , F_λ denotes the matrix form of the tensor ∇_λF_μν, and are functions of , , and τ, and a summation over i, j, and k is assumed. Finally, we can express the gauge potential in terms of the electromagnetic fields from the field strength tensor as . Combining with the electrodynamics equations in curved spacetime, we can calculate the linearized equations of motion for GWs and photons (D. Boccaletti et al. 1970). At the same time, some enlightening calculations of the GZ effect can also be referred to in these papers that discuss the conversion of axions into photons (A. Hook et al. 2018; U. Bhura et al. 2024; L. Walters et al. 2024).

Given that the inverse GZ effect is most pronounced when the direction of the magnetic field is perpendicular to the direction of the GW propagation, as is well known, GWs can be expressed with the traceless transverse gauge, which means h_0i = 0, ∂^jh_ij = 0, and . The Euler–Lagrange equations of motion from the Heisenberg–Euler effective Lagrangian (2) for the propagating photon and graviton fields components, A^μ and h_ij propagating in the external magnetic field, are generally read as ∇²A⁰ = 0, , and . When we consider the case where GWs pass through the rotational axis, the above equation can be further simplified. The linearized equations of motion for GWs and photons propagating along the direction x_s, perpendicular to the background magnetic field, are derived as follows as the Coulomb gauge condition makes ∂_iAⁱ = 0, and we can also choose A⁰ = 0 (G. Raffelt & L. Stodolsky 1988; M. Maggiore 2000; A. D. Dolgov & D. Ejlli 2012),

where TT means transverse and traceless, so and , similar to the vector potential A_i. The total frequency term involves the plasma and QED effect frequencies (W. Heisenberg & H. Euler 1936; L. D. Landau & E. M. Lifshitz 1960; S. L. Adler 1971; C. Itzykson & J. B. Zuber 1980). Moreover, we also consider the cyclotron frequency ω_cyc within magnetar. In this paper, we consider cold plasma, and its frequency is where is the fine structure constant, and n_e is the electron number density. It is worth noting that the frequency of photons ω_γ needs to be greater than the plasma frequency ω_plasma to propagate in the plasma, indicating a minimum detectable frequency ω_GW > ω_plasma. The QED effect frequency is derived from the vacuum polarization tensor and corresponds to the refractive indices. For an ordinary magnetic field, the QED effect frequencies are , and . For a strong magnetic field whose B > B_critical, the Coulomb force on electrons is a small perturbation compared to the magnetic force (D. Lai 2001). Therefore, we calculate the QED effect frequency under a wrenchless field state with a pure magnetic field (J. S. Heyl & L. Hernquist 1997a, 1997b; G. V. Dunne 2004; J. O. Andersen et al. 2016). These frequencies are denoted as and , where , , ρ = B_critical/∣B∣, and ζ(z, a) =  is the Hurwitz zeta function (E. Elizalde 1986; see also a brief description in Appendix A). It is worth noting that, for weaker magnetic fields and lower detection frequencies, it is calculated that the QED effect can be ignored. But with an increase of magnetic field intensity and detection frequency, the QED effect gradually becomes significant. In the analysis of the inhomogeneous wave equation system (Equation (3)), one can typically determine the eigenmodes of the associated homogeneous problem by employing the method of separation of variables in conjunction with Duhamel’s principle. Accordingly, we first assume that the photon propagates a short distance s along the +x_s direction, which is perpendicular to the background magnetic field. Under this assumption, the magnetic field within this region can be regarded as slowly varying and nearly uniform. That is, the temporal variation of the neutron star’s magnetic field at this short distance is much slower than the characteristic frequency of the photon. The solution of Equation (3) with a single frequency mode ω can be expressed as where λ = × and +. Hence, Equation (3) is restated as

where is the mixing mass matrix (G. Raffelt & L. Stodolsky 1988; D. Ejlli 2020; V. Domcke & C. Garcia-Cely 2021; T. Liu et al. 2024)

where B_eff represents the effective magnetic field strength related to the magnetic field distribution at the local level. Similar to the WKB limit (G. Wentzel 1926; C. Eckart 1948), Equation (4) is reduced to a linearized system with for j = 1, 2. This system yields an exact solution , where m_j represents the eigenvalue of , and U denotes the eigenvector matrix that is associated with the scattering cross section from the particle collider. Then, the GW–photon conversion process is obtained by solving with the condition , so

where .

Furthermore, the distance over which GWs propagate through the magnetosphere of a neutron star is evidently much larger than the infinitesimal displacement s, rendering the approximation of a homogeneous magnetic field invalid. The spatial inhomogeneities in both the magnetic field and the electron density lead to position-dependent coefficients in the wave equation, necessitating a perturbative treatment of Equation (3). This can be achieved by transforming to the interaction picture via , allowing the full solution to be systematically constructed through iterative expansion from an usual iteration. In practice, a first-order approximation using the distorted wave function approach proves sufficient, wherein conversion probabilities between regions of locally homogeneous magnetic field and electron density are incorporated. Finally, the probability that a GW with the frequency ω traveling a distance L at the polar angle θ of a neutron star is converted into photons is

For readability, please refer to Appendix C for the detailed calculation of Equation (7); Figure 16 shows the variation of the real and imaginary parts of this integral with the traversal distance L. In theory, the typical distance L should be influenced by the absorption and scattering of matter in the magnetic field. However, given that the interaction cross section of gravitons is with the four-dimensional Planck length l_p (S. Boughn & T. Rothman 2006; A. Palessandro & M. S. Sloth 2020), and this value is generally unaffected by the internal structure of gravitons (R. F. Sawyer 2020), coupled with the absence of extreme gravitational objects like black holes in our observation path, we infer that the absorption and scattering coefficients of gravitons are minimal and can be disregarded. Therefore, the distance L is solely determined by the graviton’s energy and the intensity, orientation, and magnitude of the magnetic field. Moreover, due to the absorption and scattering effects of the medium on photons, the conversion probability from a single GW to a photon is not equivalent to the conversion probability from a single photon to a GW. So we need to calculate the length of this typical distance L. In Figure 1, we show the minimum radius r_occur at which GWs cross the magnetar equator and radio-observable converted photons. Starting with the simplest, we can now determine the conversion probability of neutron stars for a particular path at the location of r_occur at their equators θ = 0.

Zoom In Zoom Out Reset image size

For electromagnetic waves traveling through plasma, the propagation speed of the wave energy is described as the group velocity . Obviously, in the neutron star magnetosphere, the speed of the GW and the generated electromagnetic wave are not consistent, so the electromagnetic wave generated by the GZ effect will cause interference in the path of the GW. Instead of calculating complex interference phenomena, we can turn to calculating the coherence lengths L_coherence of electromagnetic waves in different regions of the neutron star’s magnetic field (D. Fargion 1995)

In the magnetosphere and the surrounding magnetic field of a neutron star, the coherence length L_coherence is always significantly smaller than the travel distance L of the GW. Therefore, in the travel distance of GWs with the length L, the conversion probability considering the electromagnetic wave interference phenomenon can be approximated by the total number of coherent domains η_coherence = L/L_coherence and the conversion probability P_g→γ(L, ω, θ) without considering interference conditions

To estimate the order of magnitude of the conversion probability for a single GW traveling along the equator, we preselected neutron star data from the ATNF Pulsar Catalog (R. N. Manchester et al. 2005; A. Staff 2005) and the McGill Online Magnetar Catalog (S. A. Olausen & V. M. Kaspi 2014; M. Staff 2014). They are listed in Table 1. Our selection of these pulsars or magnetars is based on the distance of the neutron stars and the strength of the magnetic field. Estimating the order of magnitude of the conversion probability can be divided into two steps. First, we use this table to select neutron stars with magnetic fields B_eff of 10¹⁵ Gauss, 10¹⁴ Gauss, 10¹³ Gauss, spin periods P of 10, 5, and 1 s, and distances of 2, 1, and 0.5 kpc. We then use the information to calculate the radio signal conversed by GWs in the neutron star’s magnetic field along this path and the simple radiation mechanism. We have categorized the results of our calculations into a cross polarization and plus-polarization of conversion probability, as fully shown in Figures 17 and 18, respectively, in Appendix D for readability, and only shown in one diagram for comparison in Figure 2. Meanwhile, under identical parameter settings, we computed the conversion probabilities of GWs at various deflection angles at the radius r = r_occur of the neutron star. The calculation results for various frequencies are presented in Figure 3. Figure 3 illustrates that the conversion probability exhibits periodic variation in relation to the magnetic field distribution of neutron stars.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Because of the space-accumulation effect of the GZ effect and the change in the magnetar’s magnetic field’s direction, we can see from the two figures that the conversion probability first rises and then falls as the GW passes through the magnetar for a specific distance (F.-Y. Li et al. 2020). These two figures also show that the distance with the highest probability of GW–photon conversion is between ∼10⁶ and ∼10⁸ km when the frequency of the GW is between 10⁶ and 10¹¹ Hz. This distance exceeds the radius of a neutron star’s light cylinder, which is R_LC = cP/2π ≃ 5 × 10⁹ cmP₀ (Y.-P. Yang et al. 2020). Additionally, as the magnetar’s spin period and GW frequency increase, this distance gets smaller. It is important to note that some of the finer magnetic field structures in the neutron star magnetosphere cannot be reflected in the simple estimate because we consider the aligned GJ model in this paper. This is the main estimate we will do in the next paper and is highly worthy of further investigation.

As aforementioned, the determination of the specific intensity of the converted photon I_γ,ω at the frequency ω relies on the total conversion rate (G. Raffelt & L. Stodolsky 1988), denoting the integrated conversion probability per unit of time and solid angle, alongside the VHF GW’s differential energy fraction Ω_GW(ω). Considering the radio telescope’s frequency binning , the final specific intensity generated in the magnetic field can be expressed as

where H₀ is the Hubble constant (Planck Collaboration et al. 2020a), and represents the reduced Planck mass. Now, let us determine the magnitude of this solid angle Δθ. Considering a neutron star with a distance d from us, the region where the VHF GWs are converted into photons is spherical and has a radius of R_tot. This radius R_tot can be a local point of the neutron star’s magnetic field, a region, or the entire magnetic field. The size and choice of this radius will greatly affect the size of the final conversion probability. To determine the size of a R_tot, it is most convenient and optimistic to analyze an integral region of infinite size. This allows for the calculation of the maximum conversion probability and the optimal utilization of the telescope’s field of view. However, this is actually not a rational assumption. When the region is of considerable size, the magnetic field of the neutron star will weaken significantly, eventually being superseded by other magnetic fields in astrophysics. Consequently, we can no longer solely attribute the observational effect we are currently investigating to the neutron star. Hence, two factors influence the magnitude of a signal’s field of view: First, the strength of the magnetic field must be higher than the average magnetic field of the Universe B₀ ≲ 47 pGauss, which allows the magnetic field region to extend to a radius of up to 10⁹ km (A. Neronov & I. Vovk 2010; F. Tavecchio et al. 2010; R. Durrer & A. Neronov 2013; K. Takahashi et al. 2013; P. A. R. Ade et al. 2016; M. S. Pshirkov et al. 2016; K. Jedamzik & A. Saveliev 2019; V. Domcke & C. Garcia-Cely 2021). Second, it is evident that the signal’s effective field of vision cannot surpass the telescope’s field of view (see our discussion of this in Appendix B). Therefore, we can consider the following simple model: The telescope is centered on a neutron star, and the data at the back end of the observation are actually the integrated total voltage of all radio signals across the telescope’s field of view (A. R. Thompson et al. 2017). It is clear that R_tot is greater than r_occur for actual observation. Therefore, we can determine the total energy flux reaching Earth

Combining these points, we calculated the integral conversion probability with the change of size of the conversion region R_tot. The results are presented in Figure 4, utilizing the same parameters as in Figure 2. The complete results are available in Figure 19 in Appendix D. The absence of data points in the figure at lower frequencies is attributed to the radius at which conversion occurs, r_occur, being greater than the radius of the conversion region, R_tot, assumed in our calculation. The figure demonstrates that the decrease in the total conversion probability near the neutron star’s light cylinder, R_CL ≃ 2.5 × 10⁵ km, is due to changes in the topological configuration of the magnetic field in that region. Simultaneously, it is observed that, when the radius of this region exceeds 1 × 10⁸ km, the total conversion probability curve approaches a flat trend. This phenomenon occurs because the magnetic field of the neutron star decays with the third power of the radius, leading to the superposition of the conversion probability from the larger region being considered a higher-order term. Finally, we selected the radius of the conversion region as R_tot = 10⁹ km. The radius is less than the field-of-view radius of the telescope, and the magnetic field strength exceeds 47 pGauss. The selected conversion region represents an optimistic estimate, and actual observations are likely to differ significantly: (a) Due to the magnetic field crossing of the surrounding objects, the magnetic field of the neutron star does not have such a clear dividing line. (b) We assumed that the physical process of converting VHF GWs into electromagnetic waves occurs uniformly throughout the conversion region. However, VHF GWs may be localized to a limited area within this conversion region, necessitating observations to accurately identify the emission region of electromagnetic waves.

Zoom In Zoom Out Reset image size

Recent studies have revealed that, when nanohertz GWs pass through the Universe, they will produce weak interference modulation due to the diffraction effect of a large number of galactic disks, and the strain amplitude will change about (D. L. Jow & U.-L. Pen 2025). Similarly, when a VHF GW passes through a neutron star’s magnetic field, the GW is also diffracted by the neutron star, which we discuss in this subsection and give the amount of change in the modulated strain amplitude.

A solar mass neutron star with a radius of 10 km has a Schwarzschild radius of about 3 km < r₀, and the Fresnel length can be calculated to determine whether GWs will be diffracted as they pass by the neutron star. The Fresnel length, which relates to the GW frequency and the distance between the neutron star and radio telescope, can be described as (R. Takahashi & T. Nakamura 2003; J.-P. Macquart 2004; K. S. Thorne & R. D. Blandford 2017; H. G. Choi et al. 2021)

Here, d_eff = d_Ld_LS/d_S ≈ d_L denotes the effective angular diameter distance to the lens, where d_L is the distance from the neutron star to the Earth. This approximation holds under the assumption that the incident GWs can be treated as plane waves, implying that their sources either lie beyond the Milky Way or originate from the very early Universe. Then, the condition for diffraction is (R. Takahashi 2006; M. Oguri & R. Takahashi 2020; H. G. Choi et al. 2021). It is clear that GWs in the radio band 10⁶–10¹⁰ Hz combined with distant neutron stars satisfy the diffraction condition. The frequency dependence in the diffraction regime is measured by the shear γ(x) from the complex amplification factor with Kirchhoff–Fresnel integral F_amp(ω) =  (T. T. Nakamura & S. Deguchi 1999). Here, r denotes the physical displacement vector on the lens plane, measured from the center of the lens, and r_s represents the projected position of the source onto the same plane. The quantity T_d is the arrival-time difference between the deflected trajectory passing through r under the influence of the lens and the unperturbed, straight-line path that would occur in the absence of the lens. In an analogy with the energy density spectrum of GWs, reformulating the Kirchhoff–Fresnel integral in terms of a logarithmic frequency interval proves advantageous for the analysis that follows (H. G. Choi et al. 2021)

To obtain a fast estimate of the cumulative diffraction amplification factor, we consider only the single diffraction flux formed by the contribution of the isolated neutron star. The diffraction amplification factor, which is the energy amplification by the GW due to the diffraction of the neutron star, is the bridge between before and after GWs are diffracted

where h_diff(ω) denotes the lensed waveform, while h(ω) represents the unlensed waveform. And both are expressed in the frequency domain with respect to the frequency ω. Therefore, the GW sensitivity we observe needs to be divided by the diffraction amplification factor F_amp(ω) (D. L. Jow & U.-L. Pen 2025). Given that the spatial scale of the neutron star magnetosphere under consideration greatly exceeds the Fresnel scale, the characteristic transverse size of the lensing structure is much larger than the Fresnel radius. Under such conditions, the lens can be effectively treated as one-dimensional, and the resulting diffracted flux is well approximated by the diffraction pattern of a one-dimensional Gaussian lens

For example, if the distance parameter d_L = 2 kpc of PSR J0501+4516 is selected and the frequency ω is 10⁶ Hz, F_amp is equal to 0.176. It should be noted that GWs in the radio band are still tensor waves, and the diffraction formula we use is a reasonable approximation (S. R. Dolan 2017; S. Hou et al. 2019; J.-h. He 2020; D. L. Jow et al. 2020; Z. Li et al. 2022). In principle, a more precise tensor diffraction formula is needed to describe the properties of tensor waves, which we will discuss in a subsequent paper.

In this subsection, we will discuss the line-broadening mechanism of gravitons. The order of discussion is collision broadening, the Zeeman effect in a magnetic field, and natural broadening.

The Lorentz transformation allows us to represent the GW amplitudes h₊ and h_× under the GW’s transverse and traceless gauge

If we assume rotations around the axis, the combinations h_× ± ih₊ undergo a transformation

Another characteristic of GWs is that they adhere to the massless Klein–Gordon field equation

Therefore, we can examine the massless, one-dimensional representation of the Poincaré group, which is defined by the condition of four-momentum P_μP^μ = 0 and is characterized by a specific value of the helicity h. The helicity h is defined as the projection of total angular momentum onto the direction of motion, expressed as , where represents the propagation direction of the GW, and J = L + S denotes the total angular momentum, which is the sum of the orbital angular momentum L and spin angular momentum S. Upon rotating the direction of motion by an angle ψ, the helicity eigenstate ∣h〉 is transformed as follows:

By combining Equations (17) and (19), the helicity of GWs, as well as that of massless gravitons, can be determined to be h = ±2. In conclusion, gravitons, massless particles with a helicity ±2 that propagate along null geodesics, are compatible with the linearization of general relativity and contemporary quantum field theory. This indicates that their theoretical speed ought to equal the speed of light, denoted as c. Observations have verified this, indicating that the speed of gravitational waves is nearly equivalent to the speed of light (A. E. Romano & M. Sakellariadou 2025).

Because the massless gravitons propagate along null geodesics, there is no change in their velocity, and therefore no Maxwell velocity distribution law for them. Hence, the frequency distribution f_d(ω) at the temperature T broadened by the collision of spectral lines is invalid

where , where ω₀ represents the fundamental frequency, and the most probable speed is given by . Simultaneously, we observe that several studies suggest that gravitons can have mass (D. M. Eardley et al. 1973; B. P. Abbott et al. 2016b, 2017a, 2017b, 2019b; M. Isi et al. 2019; E. Payne et al. 2023; T. Flöss et al. 2024; Y. Hatta 2024). The distribution of collisional broadening can be derived in the presence of mass using Equation (20). Then, the full width at half-maximum (FWHM) points of massive gravitons can be determined.

where k_B signifies the Boltzmann constant, and T_e indicates the temperature of the neutron star magnetosphere. Figure 5 presents the FWHM for massive gravitons. The mass of gravitons is determined to be m_g = 10⁻⁹–10⁻⁴ eV, corresponding to frequencies of ω ≈ 10⁶–10¹¹ Hz. Furthermore, in Figure 6, we present the constraints imposed by current graviton observations and experiments on its mass range. The quality range limit of the method discussed in this paper is denoted by a brown slash. The quality limit range of our method is somewhat narrower than that of other experiments, featuring a marginally higher upper limit. The final observational results reveal distinct characteristics that enable the determination of various graviton masses and the analysis of gravitational theories.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Furthermore, the helicity of massless gravitons causes the Zeeman effect to inhibit the spectral lines from splitting into multiple transitions among their magnetic sublevels, as transitions with Δm = 0 and Δm = ±1 are permitted according to selection rules. The shifted transitional frequency Δω_Z between magnetic sublevels equals zero

where J₁ represents the angular momentum of the upper state, and J₂ denotes the angular momentum of the lower state. m₁ and m₂ represent the magnetic quantum numbers. Besides, μ_B denotes the Bohr magneton; ℏ represents the reduced Planck’s constant; g₁ and g₂ signify the g-factors of the two states; and B indicates the magnetic field. Next, we can calculate the natural line width using the uncertainty principle

where ΔE is represented as hΔω_nat, while Δt relates to the spontaneous emission rate of gravitons as per Bohr’s correspondence principle. Given the current lack of knowledge regarding the internal structure of gravitons, and considering that gravitons exhibit a helicity of ±2, it can be inferred that Δt approaches infinity. Consequently, the natural frequency, Δω_nat, tends toward zero, rendering the natural broadening negligible.

In conclusion, we believe that, in traditional general relativity, the broadening of massive gravitons is negligible, and massless gravitons exhibit no spectral line broadening. Consequently, the spectral line broadening of the photon signal transformed by the inverse GZ effect is minimal and can be effectively considered negligible. We can differentiate the radio signal from other typical astrophysical processes once we ascertain its shape and specific spectral line attributes. In the subsequent subsection, we will examine the frequency-dependent characteristics of this signal.

In contrast to gravitons, which exhibit minimal interaction with matter due to their exceedingly small cross section, photons traversing the Universe undergo scattering by matter or celestial entities, including electrons and protons. Therefore, it is essential to consider the energy depletion of photons in radio observations. To maintain generality, we concentrate exclusively on the scattering of photons converted by GWs through electrons during their propagation, as the energy variation of an electromagnetic wave is directly associated with its trajectory. In a medium of small optical thickness, the probability of photon scattering through the medium is given by . According to the law of large numbers, when a photon traverses the same medium M times (M → ∞), the total number of scattering events is Mτ_es. Consequently, the number of scattering events occurring through the medium only once is . In an optically thick medium, the number of times photons are scattered, denoted as N, can be calculated using random walk theory as , where τ_es = n_eσ_TR. Moreover, n_e represents the electron density. The Thomson scattering cross section is given by , where is the classical electron radius, and R signifies the medium scale. Consequently, for a medium with any optical thickness, the typical scattering number can be represented as . The Compton y-parameter serves as a metric for assessing the significance of scattering. The total energy change is determined by multiplying the average energy change per scattering event by the total number of scatterings. . Subsequent articles will analyze the interstellar scintillation of a pulsar to develop a more accurate model of energy loss in signal propagation. We choose representative electron densities n_e = 0.08 cm⁻³ and n_e = 10⁻⁴ cm⁻³ with electron temperatures for a diffuse ionized component of T_e = 10⁴ K and T_e = 10⁶ K. These data are utilized to plot Figure 7, illustrating the total frequency variation of converted photons from the magnetar to the telescope as a result of Compton scattering. This suggests that targeting neutron stars with elevated electron temperatures and densities along the path to Earth can enhance the distinguishability of the GW signal.

Zoom In Zoom Out Reset image size

Radio data processing typically entails the extraction of signals from radio observation data and the differentiation of various target sources through the utilization of polarization information from radiation sources. Four Stokes parameters can characterize the polarization features of a quasi-monochromatic electromagnetic wave

where and , assuming the scalar potentials are set to zero. The parameter I denotes total intensity, while Q and U represent linear polarization and its position angle, respectively. The symbol V signifies circular polarization. Utilizing the four Stokes parameters, the observed polarization angle is defined as , the linear polarization is given by , the circular polarization is expressed as , and the overall degree of polarization is represented by . Additionally, a four-vector can be defined using the Stokes parameters, which can be transformed by the generalized Faraday rotation tensor ρ_αβ (D. B. Melrose & R. C. McPhedran 1991; V. N. Sazonov 1969), incorporating Faraday rotation and conversion coefficients and the absorption tensor η_αβ in relation to the absorption coefficients of Stokes parameters (D. B. Melrose & R. C. McPhedran 1991; V. N. Sazonov 1969). This results in the general radiation transfer equation

where represent the vector of four Stokes parameters. The indices i = j = 1, 2, 3, 4 correspond to the components labeled as I, Q, U, V, while _i denotes the coefficients of spontaneous emission.

However, the polarization angle Φ is subject to the influence of Faraday rotation, which is characterized by the rotation measure (RM). Faraday rotation is an optical effect where the polarization plane undergoes a linear rotation with the square of the wavelength , expressed as . The RM is directly linked to the magnetic field aligned with the line of sight (LOS), considering the free electron density integrated along the path from the source to the observer. There are several methods for measuring the Faraday rotation, or RM of polarized astrophysical signals, such as RM synthesis, wavelet analysis, compressive sampling, and QU-fitting (X. H. Sun et al. 2015). In this paper, we employ RM synthesis along with Faraday synthesis. RM synthesis (B. J. Burn 1966; M. A. Brentjens & A. G. de Bruyn 2005) is a robust technique for quantifying Faraday rotation, akin to a Fourier transformation

where ϕ represents the Faraday depth; an extension of RM is used when the polarized signal is subject to varying degrees of Faraday rotation. Repeating this process for various ϕ values produces a Faraday dispersion function that shows the polarized intensity at different test levels. When RM synthesis is applied to emission spread over a large area of space, it often results in a complicated Faraday dispersion function with significant polarized emission at various Faraday depths (C. S. Anderson et al. 2016; J. M. Dickey et al. 2019). At times, the polarization effect from the detected source varies over time, requiring the use of the polarization position angle to analyze the results. The polarization position angle differs from the observed polarization angle Φ since it describes the orientation of the polarized signal before being affected by Faraday rotation. We can use the observed RM as a multiplicative phase factor to counteract the rotation of the spectrum and eliminate the impact of Faraday rotation

where [Q + iU]_obs is the observed spectrum, and [Q + iU]_source is the intrinsic polarization vector at the source, while RM and Φ₀ are fitted parameters. Φ₀ is the polarization position angle at a reference wavelength . In the case of calibrated polarized observations, Φ₀ is often referenced at infinite frequency where Faraday rotation is zero. In principle, any time dependence of Φ₀ can be determined by fitting the polarized signal through the burst duration. This time-resolved analysis is hard to do in practice because of S/N limitations. It is also not good for an automated pipeline that needs reliable ways to describe the polarized signal. An alternative method for characterizing time dependence in Φ₀ is to apply the observed polarization angle to the burst profiles of the derotated Stokes Q and U parameters , where Q_integ and U_integ are integrated over frequency to optimize the signal-to-noise of the Φ₀ measurement under the assumption that there is no frequency dependence in the intrinsic polarization angle at the source. Calculating the Φ₀(t) curve in this way makes it less sensitive to measurement errors associated with Stokes Q and U, yielding a more stable curve through the observational duration. In addition, it should be noted that the wind of magnetars can affect RM (M. Lyutikov 2022), but because it requires a relatively dense and slow wind, we will not consider this effect in this paper for the purpose of simplifying the model.

By employing the relationship between A_ω,λ(t, s) and h_ω,λ(t, s) in conjunction with RM synthesis, we can discern various sources of VHF GWs within the signal and eliminate some astrophysical signals resembling GW signals. This constitutes our polarization criterion. Since these radio signals are essentially conversed in resonance with GWs, the initial polarization of radio signals is exactly the same as that of GWs, resulting in linear polarization. By measuring the electron number density and magnetic field strength of the neutron star in the direction of the radio telescope’s LOS, we can obtain the RM of the signal RM = . We can derotate the spectrum by multiplying the phase factor, which helps to restore its shape prior to Faraday rotation. From this, we can find the ratio relationship between [Q + iU]_source and [Q + iU]_obs that is [Q + iU]_source/[Q + iU]_obs =, where is the frequency after the frequency shift. Taking into account the values of the parameters in Appendix B and the calculation results of the frequency drift shown in Figure 7, we can plot the result of the change of the real and imaginary parts of this ratio with the distance and frequency of the observing source in Figure 8, respectively. Figure 8 illustrates that, at lower observation frequencies, the vibrations of the real and imaginary components of [Q + iU]_source/[Q + iU]_obs are pronounced. Furthermore, increasing the distance does not mitigate this vibration, and the calculated results for varying distances are challenging to differentiate. At higher observation frequencies, an increase in distance results in more pronounced changes in the real and imaginary components of [Q + iU]_source/[Q + iU]_obs with respect to frequency. This analysis indicates that a higher frequency should be selected for the polarization observation of VHF GWs.

Zoom In Zoom Out Reset image size

The system-equivalent flux density (SEFD) is a figure of merit that characterizes the sensitivity of a radio telescope within a certain frequency band , where T_sys(ω) = T_AST(ω) + T_ATM(ω) + T_RT(ω) represents the frequency-dependent system temperature, and denotes the telescope gain, where A_eff = η_AA_phys is the effective collecting area. The effective collecting area is determined by multiplying the physical antenna aperture A_phys with an aperture efficiency η_A = η_sfη_blη_sη_tη_miscη_sloss. Here, η_sf denotes the reflection efficiency of the main reflector, which is governed by the Ruze equation , where ε is the rms surface error, and λ is the observational wavelength. The blockage efficiency, η_bl, accounts for the feed cabin obstruction. η_s represents the spillover efficiency, while η_t describes the illumination efficiency. η_misc includes various minor loss factors such as phase center offset and impedance mismatch. The factor most sensitive to the telescope pointing is η_sloss, which reflects the reduction in effective aperture area with zenith angle θ_ZA. Since the other efficiency terms remain nearly constant during observations, the variation in the total aperture efficiency η is primarily driven by η_sloss.

The main system temperature receives contributions from (a) astrophysical backgrounds , where τ(ω) is the atmospheric zenith opacity, and z_e is the elevation of the telescope, (b) the atmosphere , and (c) the telescope electronics T_rt(ω) = T_spill(ω) + T_radio(ω). In this paper, we consider the following radio telescopes: Five-hundred-meter Aperture Spherical Radio Telescope (FAST; P. Jiang et al. 2019, 2020; L. Qian et al. 2020; F. Staff 2021), TianMa Radio Telescope (TMRT; Z. Shen 2014, 2017; T. Staff 2021, 2022), Qitai Radio Telescope (QTT; Q. Staff 2014; N. Wang et al. 2023), Square Kilometre Array Phase 1 MID-197 (SKA1-MID; R. Braun et al. 2019; S. Staff 2019), Square Kilometre Array Phase 2 MID-2000 (SKA2-MID; S. Staff 2019), and Green Bank Telescope (GBT; G. Staff 2017; G. S. Staff 2017). We use several single-dish telescopes and radio arrays, which can be superimposed by considering the beamwidth 2θ_b and baseline lengths of a single-dish telescope. The beamwidth of a radio telescope is proportional to the wavelength of radiation under consideration, with the FWHM θ_b obeying the ratio , where D represents the telescope diameter. The individual beam elements, which are composed of dishes of the diameter D, have approximate gains . In our simplified treatment, the primary beam radius is simply . On the other hand, the synthesized beam diameter θ_synth = 2r_synth is determined by the maximum baseline of the array . We summarize these parameters of the telescope, as well as the number of antennas and the altitude of the telescope station site, in Table 2.

For astrophysical backgrounds , and T_ast =  T_cmb(ω) + T_rsb(ω) + ΔT_source(ω), it usually consists of the nearly isotropic CMB T_cmb ≈ 2.73 K, the average sky brightness temperature T_rsb(ω), which can be estimated from the Haslam map (C. G. T. Haslam et al. 1981, 1982), Wilkinson Microwave Anisotropy Probe (G. Hinshaw et al. 2013), Planck (N. Aghanim et al. 2020; Planck Collaboration et al. 2020b, 2020c, 2020d), and Beyond Planck data (T. L. Svalheim et al. 2023) assuming a power-law spectrum T_rsb(ω) = T₄₀₈(408 MHz/ω)^β with T₄₀₈ = 17.1, 25.2, 54.8 K for the 10th, 50th, and 90th percentiles of the all-sky distribution when gridded in an equal area projection at zero zenith angle, and ΔT_source(ω) is from the astronomical source being observed, which is related to the target source in the actual observation, and the value is so small that it can be ignored in theory. For example, in the v_RF ≈ 4.85 GHz sky survey made with the 300 foot telescope, the system noise was T_sys ≈ 60 K, but the faintest detected sources added only ΔT_source ≈ 0.01 K. Moreover, take the FAST telescope as an example. For a representative pulsar with an average flux density of 0.1 Jy, the resulting brightness temperature contribution is approximately 1.54 K. This is significantly smaller than the typical system temperature of several tens of kelvin, and even for exceptionally bright pulses, the time-averaged effect remains minor. Therefore, it is justified in most cases to assume ΔT_source ≪ T_sys when evaluating the system sensitivity in pulsar or continuum observations. For low-frequency class (LFC) telescopes like FAST, we can take into account the sky brightness temperature with the parameter β = 2.76 with the frequency from 500 MHz to 3 GHz (P. Platania et al. 1998) and the parameter β = 2.51 with the frequency from 1 to 500 MHz since the Galactic plane is a bright diffuse source at low frequencies (Y. Cong et al. 2021; M. O. Irfan et al. 2022). One can use the Ultralong-wavelength Sky Model, which can simulate the radio sky at frequencies below 10 MHz to obtain the low-frequency sky brightness temperature (Y. Cong et al. 2021). For medium frequency class (MFC) telescopes like the TMRT telescope, we can take into account the sky brightness temperature with the parameter β = 2.76 at the frequency ω ∈ [3, 10] GHz (P. Platania et al. 1998) and β = 3.00 at the frequency ω ∈ [10, 14] GHz (J. L. Weiland et al. 2022). For high-frequency class (HFC) telescopes like SKA1-MID, we can take into account the sky brightness temperature with the parameter β = 3.25 (J. L. Weiland et al. 2022). For very high-frequency class (VHFC) telescopes like GBT, we can take into account the sky brightness temperature with the parameter β = 3.25 at the frequency ω ∈ [30, 70] GHz (T. L. Svalheim et al. 2023) and β = 3.30 at the frequency ω ∈ [70, 120] GHz (T. L. Svalheim et al. 2023).

The is the surface air temperature in the telescope beam, which is related to the atmospheric environment of the real-time observation day. However, we can estimate it theoretically. We consider a layer of water vapor at a temperature T_atm and an atmospheric zenith opacity τ(ω) at a given observing frequency. The total optical depth is usually influenced by hydrosols τ_hy, oxygen τ_ox, dry-air τ_da, and water vapor τ_wv. And the oxygen and dry-air opacities are nearly constant, while the water vapor and hydrosol contributions vary significantly with weather. We use the atmospheric transmission at microwaves (ATM) model (J. R. Pardo et al. 2001) code as implemented within the Common Astronomy Software Applications (CASA; CASA Team et al. 2022) package to obtain atmospheric zenith opacity τ(ω). For all the telescopes, we set the water vapor scale height at 2 km without loss of generality. For the telescopes that have been built, we use local weather stations to obtain 10th, 50th, and 90th percentiles of data on the average annual precipitable water vapor (pwv), average annual temperature, average annual atmospheric pressure, and average annual humidity as our predictions. For the incompletely built SKA2-MID and QTT telescopes, we use the atmospheric prediction values consistent with SKA1-MID and the atmospheric data of Qitai Observatory, respectively. We present our simulations of the zenith atmospheric opacity of these telescopes as a function of observation frequency in Figure 9. And the three different colors show three different average annual pwvs. In the section results, we present the final results using the 50th average annual pwv.

Zoom In Zoom Out Reset image size

For the telescope electronics temperature T_rt(ω), we can estimate from the spillover radiation T_spill(ω) that the feed picks up in directions beyond the edge of the reflector, primarily from the ground. This contribution originates from the side lobe of the radio telescope. Unlike the other single dish with a fixed surface and horn, the illumination area of the horn varies for different zenith angles, leading to different values of T_spill(ω). And the radiometer noise temperature T_radio(ω) attributes to noise generated by the radiometer itself. For example, we can set T_spill(ω) = 3 K at all frequencies at zero zenith angle and from 0.35 to 1.05 GHz for SKA1-MID (R. Braun et al. 2019).

Finally, we present the simulation results of the telescope’s zenith atmosphere temperature with observation frequency in Figure 10. We use different colors to show the three different average annual pwvs, just like with zenith atmosphere opacity. For our final results, we use the simulation results of the system temperature at the 50th astrophysical background temperature and the average annual pwv. These telescope parameters are summarized in Table 4 in Appendix B based on the size of the telescope.

Zoom In Zoom Out Reset image size

As we will show in Section 5, we use two types of radio telescopes to discuss the signals of two different timescales, so we need to calculate four S/Ns. In this subsection, we discuss the four S/Ns separately.

For a single-dish telescope such as FAST, in the persistent event case, the S/N is equal to the flux F of the radio sources divided by the error of the measurement ΔF:

where n_pol denotes the number of polarization channels, and Δν signifies the frequency bandwidth associated with each radio telescope’s observations of neutron stars; and Δt_int indicates the integration time Δt_int = t_surτ, where t_sur represents the total time survey duration, and τ = λ_γ/2πD is the fractional time that an object on the celestial equator transits the field of view of the telescope; and SEFD is the SEFD mentioned earlier (H.-R. Yu et al. 2014). Therefore, it is evident that a longer observation time for the same target source will result in a higher S/N. However, during the actual observation, the observer will not have unlimited time. The location of the target source, its proper motion, and the radio telescope’s performance all influence the accessible observational time. We can select the target source based on the actual situation to extend the observational time, but the telescope cannot perform as expected; for example, for FAST, which is not a fully steerable telescope and has the largest target tracking time approximately 6 hr. Therefore, we estimate the observational time for all telescopes used to observe to be 6 hr (P. Jiang et al. 2020). For a transient event, the formula for calculating the energy flux becomes an integral, , and the integral time changes from the observation time to the duration of the signal (D. R. Lorimer et al. 2013). This limits the minimum sampling time of the telescope, and the approximate order of magnitude of the duration of the signal is mentioned in Section 5 and in the second subsection of Section 2. For example, this means that observations of transient signals cannot be received by using FAST’s SETI backend, which has a sampling time of 10 s (B.-L. Huang et al. 2023; X.-H. Luan et al. 2023; Z.-Z. Tao et al. 2023, 2022).

For a radio interferometer like Square Kilometre Array (SKA) in the persistent event case, we treat the signal as a collection of synthesized beams. Therefore, we start with the assumption that the observed signal is Gaussian, that each signal is independent of the other, and that there is only one kind of GW. And these assumptions are sufficient for discussion. Then, we can obtain the likelihood function of the observational energy flux F_obs,i, accounting for the measurement error σ_i and theoretical model energy flux F_model,i

where θ_set is an idiographic set of parameters such as the one-sided power spectral density of GWs S_h(ω), the amplitude of GWs Ω_GW, the characteristic strain of GWs h_c, and so on. And the χ² statistic is

where F_sys,i is the predicted background flux density for the single synthesized, which includes the SEFD and other known astrophysical processes of flux density. The single-dish error ΔF can be used to estimate the measurement error σ_i. Therefore, we can obtain the S/N from the χ² statistic (A. Hook et al. 2018)

where represents the array gain for that particular synthesized beam, based on the individual antennas and the number of array elements N. For each synthesized beam, there is a single signal flux F_model,i, a single integrated time Δt_int,i, a single beam’s frequency resolution Δν_i, and a single system temperature T_sys,i. We are aware that another possibility is the simultaneous observation of mixed multiple GW signals, a situation we will discuss in our future research, which will require significant effort. Because statisticians have discussed the construction of a likelihood as the product of sublikelihoods under the term composite likelihood inference, they discover that the statistical foundation is not solid. And there does not seem to exist a general strategy for constructing a combination of composite likelihoods that is both computationally convenient and statistically appealing. One can refer to these papers (B. Lindsay 1988; Y. Zhao & H. Joe 2005; F. Larribe & P. Fearnhead 2011; B. G. Lindsay et al. 2011; C. Varin et al. 2011; W. Hong et al. 2023). Finally, we modify the single signal flux F_model,i and the single integrate time Δt_int,i for radio interferometers in transient events, respectively, to the single peak density flux F_peak,i and the single signal duration Δt_sam,i.

With these preconditions, we can simulate the expected signal shapes for both timeframes. We use FAST as our telescope to generate simulation signals. The telescope’s and neutron star’s parameters for generating the simulation signal with FAST are as follows: The center frequency of the telescope is 1420 MHz, the bandwidth is 500 MHz, the number of sampling channels is 4000, and the data sampling time is 98.304 us. For a magnetar or pulsar, we set its magnetic field as B = 1.85 × 10¹⁴ Gauss, its spin period as , and the distance d = 2 kpc exactly as PSR J0501+4516. Since the energy of the final signal requires studying the SEFD of the telescope, we will discuss the parameters of the GW in Section 5. As discussed earlier, for the duration of the signal, the transient signal depends on the conversion probability with the distance the GW travels through the magnetic field, while the persistent signal is present throughout the entire observation process. We inject the 1σ Gaussian noise in the simulation process, and the simulated signal energy flux is at least 5σ relative to the FAST minimum detectable energy flux. The results of the simulation are shown in Section 5.2. After determining the shape of the radio signal that we have observed, we can proceed to analyze the specific spectral line properties of the source responsible for generating the signal.

Reviewing the discussion in Section 2.1, the distance at which this signal is generated is much greater than the radius of the optical column surface of the neutron star. And this distance decreases as the frequency of the GW increases, and the spin period of the magnetar increases, resulting in a signal that looks like a fast radio burst, but its signal duration is much longer than the fast radio burst (B. Zhang 2023). Therefore, we refer to this signal as a “novel rapid optical variation signal.”

For this “novel rapid optical variation signal,” from the conversion probability of GWs to photons, as long as the GW enters the magnetic field range of the neutron star, photons can be generated, so the overall signal can theoretically be an integral result. However, if the space-accumulation effect of the GZ effect is insufficient, it leads to a small number of produced photons, low-energy flux, and a faint signal. With the increase of GW frequency, the height of electromagnetic wave emission decreases, the magnetic field increases, and the conversion probability increases. Depending on the behavior of GWs, signals can be divided into transient signals and persistent signals, which are going to be discussed in Section 5. Therefore, without considering the dispersion in the propagation process, for a transient signal, after the observed frequency is fixed, the complete signal should present a likely bell-shaped signal on a two-dimensional graph where the horizontal axis is time, and the vertical axis is energy flux. The bell-shaped signal exhibits asymmetry in the distribution of peak energy flux between its left and right sides. The slope of the signal curve to the left of the peak is significantly greater than the slope of the signal curve to the right of the peak, which is estimated by the variational trend of the conversion probability. And the timeline shows that the signal to the left of the peak lasts much less than the signal to the right. In the whole observation frequency, on the two-dimensional heat map, the horizontal axis is the time, the vertical axis is the observation frequency, and the color bar of the heat map represents the signal intensity. The heat map should display a trapezoid signal, with the intensity increasing as the frequency increases. Similarly, when neglecting the effects of dispersion during signal propagation, for a persistent signal, after the observational frequency is fixed, the complete signal presents an overall integrated profile composed of transient signals on a two-dimensional graph where the horizontal axis is time, and the vertical axis is energy flux. After carefully deducting the effects of background noise, it should be an almost straight line, with only some slight intensity perturbation. In the whole observation frequency range, a trapezoidal signal should be presented on the two-dimensional heat map with time as the horizontal axis and observation frequency as the vertical axis, and the color bar of the heat map represents the signal intensity, and the intensity increases with the increase of frequency.

The primary difference among various radio telescopes regarding signal duration is in the transient signal. Due to the different sensitivity of the telescope, the time of the transient signal detected by the high-sensitivity telescope is longer. Since there is a signal whenever a GW enters the neutron star’s magnetic field, it can only be observed when the signal strength is strong, so the duration of the signal is calculated based on the conversion probability in Figures 17 and 18, the energy density of the GW, and the neutron star’s own properties such as rotation period, magnetic field strength, and other characteristics. Based on the current calculation results, a fixed characteristic amplitude h_c of the GW leads to an increase in the transient signal’s duration as the neutron star’s magnetic field intensity increases.

Furthermore, for the spectral index, we can use S ∝ αν^β to model the change of the energy flux S with the frequency ν, where β is the spectral index, approximately equal to 1, α corresponds to the GW energy density , and the neutron star conversion probability is a fitting constant. We can theoretically determine the conversion probability of the observational neutron star first, which enables us to derive the change in GW energy density with respect to frequency.

Combining the spatial-accumulation effects of the GZ effect, we endeavor to take full account of VHF GWs passing through magnetic fields around neutron stars. By doing so, we can increase the detectable field of view and then improve our detection capabilities along the trajectory along which GWs travel. The photons undergo an extensive and complex radiation process before being captured by radio telescopes. The process of VHF GWs and photons traversing the magnetic field, distinguishing the signals, and the motivation for detection are illustrated in Figure 11. The entire process is as follows: (A) The intense magnetic field of the neutron star creates an effective setting for the GZ effect to converse GWs into electromagnetic waves. The spatial-accumulation effect results in coherent superposition of the conversion probability of GWs as they traverse the magnetic field along the magnetospheric path. The impact of the residual magnetic fields along the observation path of neutron stars can be categorized into interstellar magnetic fields (V. Domcke & C. Garcia-Cely 2021) and planetary magnetic fields (T. Liu et al. 2024). In comparison, the intensity of electromagnetic wave signals converted from GWs in these magnetic fields is relatively low. The signal intensity produced by the interstellar magnetic field is 10² times lower than that conversed by the neutron star magnetic field, while the signal intensity from the planetary atmospheric magnetic field is also 10⁴ times lower than that of the neutron star magnetic field. The effect on our observations is minimal. The propagation path of the generated electromagnetic wave is influenced by the density distribution of the magnetospheric plasma. Free electrons induce a frequency shift in electromagnetic waves via Compton scattering, whereas the Faraday rotation effect alters the polarization angle, leading to variations in Stokes parameters associated with the magnetic field direction. The signals obtained from radio telescopes (B) and (C) require multilevel processing to delineate the characteristics of GWs. The performance variations among various radio telescopes directly influence the detection thresholds. This paper investigates the use of the equivalent flux density of the telescope system as a reference for background noise. It employs RM synthesis for polarization decomposition and a time-frequency heat map to identify the bell-shaped profile of transient signals and the broadband continuous spectrum of these signals. (D) By aligning the signal template with the observed data, the essential parameters of GWs can be determined: energy density and characteristic strain. The mass distribution of PBHs and the predicted parameters of quantum gravity theory can be constrained through multiband observations. (E) For scientific motivation, the exploration of VHF GWs, particularly in the MHz–GHz range, has opened a new avenue for probing physics beyond the standard cosmological model. Such frequency bands are largely inaccessible to current ground-based interferometers but may carry imprints of early Universe phenomena and fundamental interactions at energy scales close to the Planck scale. Here are some inspiring scientific goals. One compelling motivation arises from the inflationary paradigm. In canonical models of slow-roll inflation, the primordial tensor spectrum is predicted to be red tilted, leading to an exponential suppression of GW amplitudes at high frequencies. Consequently, any detectable relic signal in the GHz regime would challenge the standard inflationary framework. As discussed by S. Vagnozzi & A. Loeb (2022), the presence of a cosmic graviton background at such frequencies would lend support to alternative scenarios, including bouncing and emergent cosmologies, which generically predict enhanced high-frequency tensor power. Detection of such a signal would provide an empirical handle on trans-Planckian physics and the microphysical origin of the early Universe. VHF GWs also emerge naturally from the dynamics of PBHs. PBHs can form from enhanced scalar perturbations in the early Universe and may subsequently merge to generate a stochastic gravitational-wave background (SGWB) peaking in the MHz–GHz band. As shown by K. Kohri et al. (2025), the shape and amplitude of this SGWB are sensitive to the initial PBH mass distribution and formation history. Hence, constraints or detection of such a background would offer crucial insights into the physics of PBH formation and early Universe density fluctuations. In addition, topological defects such as cosmic strings, predicted by various high-energy symmetry-breaking mechanisms, provide another robust source of high-frequency gravitational radiation. Oscillating string loops can emit bursts of GWs through the formation of cusps and kinks, contributing to a distinct SGWB spectrum. The recent analysis by N. Aggarwal et al. (2025) emphasizes that GHz-band detectors could potentially identify the unique nonthermal spectral features associated with cosmic strings, thereby offering an observational probe of beyond-standard-model phase transitions. Finally, HFGWs provide a potential observational window into quantum gravity. In particular, loop quantum cosmology predicts a strongly blue-tilted tensor spectrum due to a superinflationary phase preceding standard expansion. E. J. Copeland et al. (2009) demonstrate that this leads to an enhancement of tensor modes at high frequencies while suppressing large-scale signals. Such a spectral signature, if observed, would offer direct evidence for quantum gravitational corrections to classical spacetime dynamics.

Zoom In Zoom Out Reset image size

Based on the behavior of VHF GWs and the width of the magnetic field, the radio signals can be categorized into two principal classes: (a) transient and (b) persistent signals. Transient radio signals we anticipated observing are more closely associated with GW transients generated by low-mass binary system mergers, whereas persistent radio signals align more closely with a primordial stochastic GW background from the early Universe. A photon signal is unequivocally transient when one of several GW signals traverses the magnetic fields of spheroidal celestial bodies, and the temporal gap between the current and the most recent GW signals exceeds the typical timescales of ∼100 s for magnetars or pulsars, 10 s for the Sun, and 10⁻² s for the Earth, characterized by the typical length scale t = L/c at 1 GHz. This characterization is elucidated in Section 2. The duration of the transient signal is theoretically determined by three major factors: First, the neutron star’s properties, including its magnetic field strength, rotation period, radius, and the density of surrounding charged matter, determine the conversion probability of GWs to electromagnetic waves. Second, the source of GWs determines the energy density of itself. Third, the temperature of the system, which includes the astrophysical background temperature, the atmosphere temperature, and the telescope electronics temperature, determines the sensitivity of the telescope. Therefore, a distinctly continuous signal emerges when a cohort of VHF GWs passes through the magnetic field at intervals shorter than the typical timescales or when several VHF GW signals sporadically traverse the interstellar magnetic fields.

In this research, we consider six radio telescopes: (a) single-dish telescopes FAST (P. Jiang et al. 2019, 2020; L. Qian et al. 2020; F. Staff 2021), TMRT (Z. Shen 2014, 2017; T. Staff 2021, 2022), QiTai Radio Telescope (QTT; Q. Staff 2014; N. Wang et al. 2023), and Green Bank Telescope (GBT; G. Staff 2017; G. S. Staff 2017); and (b) the radio interferometers SKA1-MID (R. Braun et al. 2019; S. Staff 2019) and SKA2-MID (S. Staff 2019), utilized for capturing radio signals transmuted by VHF GWs. By integrating the aforementioned three criteria, we simulate and depict the anticipated dual types of radio signals expected to be received by FAST in Figure 12. From the previous discussion, we can get that the system temperatures of FAST are 27.389 K in 1170 MHz and 26.829 K in 1670 MHz. Therefore, we can obtain the minimum detectable energy fluxes for persistent events from 3.093 × 10⁻⁵ Jy to 3.158 × 10⁻⁵ Jy that FAST can receive for 6 hr of observational time at a 500 MHz bandwidth with 1420 MHz as the center; these are the energy fluxes of 1σ. Obviously, the energy fluxes of 5σ are 1.547 × 10⁻⁴ Jy and 1.579 × 10⁻⁴ Jy. Correspondingly, outside 2 kpc, the energy fluxes of 5σ generated by the conversion of GWs into photons in the neutron star’s magnetic field are 1.875 × 10¹¹ Jy and 1.914 × 10¹¹ Jy. In combination with Figures 17, 18, 19, and Appendix B, we can calculate that must be at least 0.05 for the overall observed energy to have an S/N higher than 5. Similarly, the minimum detectable energy fluxes for transient events from 1170 to 1670 MHz are 1.695 × 10⁻³ Jy and 2.025 × 10⁻³ Jy. Clearly, the energy fluxes of 5σ are 8.073 × 10⁻³ Jy and 1.012 × 10⁻² Jy. Accordingly, outside 2 kpc, the energy fluxes of 5σ generated by the conversion of GWs into photons in the neutron star’s magnetic field are 1.027 × 10¹³ Jy and 1.227 × 10¹³ Jy. As with the persistent events, we can calculate that must be at least 8 for the overall observed energy to have an S/N higher than 5. To make this picture, we assume the energy of the GW stays the same over the bandwidth 500 MHz. And we also assume the characteristic strains of the GW from h_c ≈ 7.079 × 10⁻²³ to h_c ≈ 1.011 × 10⁻²² are varying with frequency for persistent events. And, we assume the energy of the GW stays the same over the bandwidth 500 MHz. And we also assume the characteristic strains of the GW from h_c ≈ 1.278 × 10⁻²¹ to h_c ≈ 8.955 × 10⁻²² are varying with frequency for transient events. It needs to be emphasized that the signal intensities shown in Figure 12 do not describe a comprehensive representation of our incoming observations. Instead, they solely illustrate the dynamic nature of the signal. Figure 12 (left panel) showcases the transient signal converted from PSR J0501+4516. The sudden, bright peak represents this signal, delineated by the red dashed line encompassing the entire signal. The strength of the signal depends on Ω_GW(ω), while the temporal width correlates with the magnetic field width of the magnetar and the frequency of the GW. In a rough estimation of a bell-like signal at the frequency of 1420 MHz, the duration of the signal to the left of the peak is , and the duration of the signal to the right of the peak is . Figure 12 (right panel) illustrates the persistent signal from PSR J0501+4516. The red dashed line also delineates the entire signal. Unlike the transient signal, this persistent signal is expected to persist throughout the entirety of the observation process and be detectable across the entire observation frequency spectrum.

Zoom In Zoom Out Reset image size

Integrating the aforementioned discourse on the two signal typologies, we contemplate the S/N for single-dish telescopes in persistent event scenarios, which equals the flux energy F of the radio sources divided by the measurement error ΔF, , where n_pol denotes the number of polarization channels, and Δν stands for the bandwidth, which we regard as the bandwidth for each radio telescope’s observation of magnetars. Δt_int = t_surτ represents the integration time, where t_sur signifies the total observation time, and τ = λ_γ/2πD is the fractional time an object on the celestial equator traverses the telescope’s field of view. The SEFD denotes the SEFD, , where T_sys(ω) = T_AST(ω) + T_ATM(ω) + T_RT(ω) is the frequency-dependent system temperature comprising (a) astrophysical backgrounds T_AST(ω), (b) atmospheric contributions T_ATM(ω), and (c) telescope electronics T_RT(ω). Subsequently, in transient event scenarios for single-dish telescopes, we should, respectively, adjust the energy density flux (D. R. Lorimer et al. 2013) F and integration time Δt_int to the peak density flux and the sampling time Δt_sam.

Additionally, in the case of persistent events with radio interferometers, we have the total autocorrelation (R. C. Walker 1989; A. R. Thompson et al. 2017; A. Hook et al. 2018; B. R. Safdi et al. 2019), where represents the array gain for a specific synthesized beam, and N denotes the number of array elements. Each term of the synthesized beam includes a single signal flux Fmodel, i, a single integration time, Δt_int,i, a single beam bandwidth, Δν, and a single system temperature, T_sys,i. Similarly to single-dish telescopes, in the case of transient events with radio interferometers, we, respectively, adjust the single signal flux F_i and the single integration time Δt_int,i to the single peak density flux F_peak,i and the single signal duration Δt_sam,i. We demonstrate the sensitivity of these two anticipated signals with PSR J0501+4516 in Figure 13, combined with some typical sources of VHF GWs (D. G. Figueroa et al. 2013; R. H. Cyburt et al. 2016; D. G. Figueroa & F. Torrenti 2017; P. Auclair et al. 2020; B. Barman et al. 2023). For persistent signals, we set the observation time to the maximum tracking time of a single FAST observation t_obs = 6 hr. Therefore, a detectable can be obtained by combining Equations (10) and (11) and the S/N

and the characteristic strain h_c(ω)

where F_i is the measured energy flux with i = 1 in 1σ of S/N and with i = 5 in 5σ of S/N. And we summarize these parameters of , h_c(ω), and F_i at 1 GHz under 6 hr of observation time in Table 3.

Zoom In Zoom Out Reset image size

Table 3. The List of Detection Capability of Radio Telescopes at 1 GHz with 6 hr of Observational Time

Telescope	Observation	1σ Detectable	5σ Detectable	Minimum Detectable Normalized	Minimum Detectable
Name	Bandwidth	Energy Flux	Energy Flux	Energy Density	Characteristic Strain hc(ω)
	(MHz)	(Jy)a	(Jy)b	(Dimensionless)c	(Dimensionless)d
FAST	500	2.98 × 10−5	1.49 × 10−4	(3.25 × 107, 1.62 × 108)	(7.20 × 10−24, 1.61 × 10−23)
TMRT	1500	1.12 × 10−4	5.64 × 10−4	(7.55 × 1010, 6.58 × 1012)	(3.47 × 10−22, 3.24 × 10−21)
SKA1-MID	300	1.17 × 10−9	5.85 × 10−9	(3.82 × 104, 1.92 × 105)	(2.47 × 10−25, 5.53 × 10−25)
SKA2-MID	300	1.15 × 10−10	5.75 × 10−10	(1.14 × 104, 5.72 × 104)	(1.35 × 10−25, 3.02 × 10−25)
GBT	800	1.66 × 10−4	8.30 × 10−4	(3.89 × 1010, 1.94 × 1011)	(2.49 × 10−22, 5.57 × 10−22)
QTT	1500	5.91 × 10−5	2.96 × 10−4	(3.21 × 109, 1.60 × 1010)	(7.16 × 10−23, 1.60 × 10−22)

Notes.

^aThe minimum detectable energy flux under the 1σ condition is derived from the ΔF estimation in Equation (28). ^bThe minimum detectable energy flux under the 5σ condition is derived from the ΔF estimation in Equation (28). ^cThe energy density in parentheses, from left to right, is the set of 1σ and 5σ, separated by commas. ^dThe characteristic strain in parentheses, from left to right, is the set of 1σ and 5σ, separated by commas.

Download table as:  ASCIITypeset image

From Figure 13, it is evident that the telescope exhibits greater sensitivity to persistent events than transient events, implying that persistent events can include a wider variety of GW signals, rendering the discovery of all such signals more challenging. Consequently, identifying a suspicious signal in the absence of a GW template is of paramount importance. We suggest employing signal cross correlation techniques (A. A. Penzias & R. W. Wilson 1965; C. J. Hogan & P. L. Bender 2001; M. Tinto et al. 2001; Y. Shao et al. 2023) to screen observational data at various observation frequencies.

To further investigate how the telescope’s detectability varies under different physical conditions, Figure 14 presents the upper bounds on VHF GW signals derived from six representative radio telescopes, for both transient (left panel) and persistent (right panel) scenarios. The solid and dashed colored lines correspond to S/Ns of 1 and 5, respectively, with the color scheme indicating the sensitivity range of each instrument. The evolution of the characteristic strain h_c for a range of GW sources is shown as a function of the magnetic field region’s radius, assuming the same neutron star magnetic parameters used in Figure 13 at 1 GHz.

Zoom In Zoom Out Reset image size

This analysis demonstrates that the predicted sensitivity is strongly influenced by the spatial extent of the coherent magnetic field. The optimistic cases correspond to magnetic regions bounded by the telescope’s full-beam field of view and the typical size of a pulsar’s magnetosphere, under the condition that the field strength remains above the interstellar background. In more conservative scenarios, the GW may only interact with a limited portion of the magnetosphere, or the effective signal may be reduced due to geometrical factors such as the finite source distance or beam coverage. These constraints reduce the conversion probability and lower the expected signal amplitude. For instance, one possibility is that the GW traverses only a small fraction of the pulsar’s magnetosphere, thereby significantly reducing the effective interaction volume over which GW–photon conversion can occur. In such cases, the interaction path length becomes comparable to or even shorter than the coherence scale of the background magnetic field, leading to a marked suppression in the conversion amplitude. Another scenario involves the regime in which only a single photon induced by the GZ effect may be generated and potentially detected, representing the quantum limit of detectability. This situation is especially relevant when the GW signal is weak, and the magnetic field strength or interaction region is limited. Moreover, some more geometrically constrained situations cannot be adequately modeled without reference to specific observational conditions. For example, the standard assumption that the GW source lies effectively at infinity may not hold in all cases. If the source is a compact object, such as a nearby stellar-mass binary, the angular extent of the incoming GW front may become nonnegligible. In such cases, only the component of the converted electromagnetic signal that is aligned with the observer’s LOS will be detectable, thereby decreasing the effective conversion probability. Similarly, for GW sources with limited spatial extent or anisotropic emission, the telescope may intercept only a subset of the generated electromagnetic signal—potentially through just one or a few beams in a multibeam system—which further suppresses the expected signal strength. Accordingly, Figure 14 provides a more comprehensive perspective on observational feasibility by bracketing both favorable and constrained conditions for a conversed photon in neutron star environments.