Searching for ultralight dark matter conversion in solar corona using Low Frequency Array data

Abstract

Ultralight dark photons and axions are well-motivated hypothetical dark matter candidates. Both dark photon dark matter and axion dark matter can resonantly convert into electromagnetic waves in the solar corona when their mass is equal to the solar plasma frequency. The resultant electromagnetic waves appear as monochromatic signals within the radio-frequency range with an energy equal to the dark matter mass, which can be detected via radio telescopes for solar observations. Here we show our search for converted monochromatic signals in the observational data collected by the high-sensitivity Low Frequency Array (LOFAR) telescope and establish an upper limit on the kinetic mixing coupling between dark photon dark matter and photon, which can reach values as low as 10⁻¹³ within the frequency range of 30 − 80 MHz. This limit represents an improvement of approximately one order of magnitude better than the existing constraint from the cosmic microwave background observation. Additionally, we derive an upper limit on the axion-photon coupling within the same frequency range, which is better than the constraints from Light-Shining-through-a-Wall experiments while not exceeding the CERN Axion Solar Telescope (CAST) experiment or other astrophysical bounds.

Introduction

Due to the absence of significant results in the search for weakly interacting massive particles (WIMPs)^1,2,3, increasing attention has shifted towards the ultralight dark matter (DM) candidates, including dark photons, quantum chromodynamic (QCD) axions and axion-like particles. Dark photon is a hypothetical vector ultralight DM candidate^4,5,6,7, constituting one of the simplest extensions of the Standard Model (SM) by incorporating a massive vector field coupled to the photon field through the kinetic mixing marginal operator^{8,9,10,11,12,13}. There are several ways to produce the right amount of dark photon dark matter (DPDM) during the early Universe, including the misalignment mechanism with a non-minimal coupling to the Ricci scalar^5,6,14,15,16, inflationary fluctuations^{7,17,18,19,20,21,22,23,24,25,26}, parametric resonances^{27,28,29,30,31,32}, or the decay of the cosmic strings³³. The QCD axion, initially introduced to address the strong CP problem as a hypothetical particle^34,35,36,37, where “CP” stands for the combination of charge conjugation symmetry and parity symmetry, has been shown to be a good DM candidate³⁸. Axion-like particles arising in, e.g., string-theory models³⁹, coupled to SM particles in a similar way, also stand as promising DM candidates. Axions or axion-like particles can be generated by the misalignment mechanism^40,41,42, or the decay of topological objects^43,44 during the early Universe.

The couplings between dark photon or axions and SM particles provide important tools in searching for these ultralight particles. Various types of experiments are looking for the signals associated with photons, including haloscopes for Galactic halo DM^45,46, helioscopes for ultralight particles emitted from the Sun^45,46, and the “Light Shining through the wall” (LSW) methods^47,48. Dark photons and axions can also be detected via WIMP detectors^49,50. Moreover, many experimental results initially intended for axion DM can be reinterpreted for dark photons. A comprehensive summary of experimental constraints (including projected ones) for dark photons and axions can be found in ref. ^51,52.

The other meaningful way to look for axions or dark photons is to investigate anomalous signals in various astrophysical environments, such as the cosmic microwave background (CMB) spectral distortion constraints on dark photons^6,53, gamma-ray constraints on axion DM^54,55, neutron stars^{56,57,58,59,60,61,62,63,64,65,66,67}, white dwarfs^66,68,69,70, supernovae^71,72, quasars and blazars^{73,74,75,76,77,78}, the Sun, red giants and horizontal branch stars^79,80,81, and globular clusters^82,83. These searches assume that the ultralight particles are either DM or sourced inside the astrophysical objects. Remarkably, the Sun, being our closest star, offers a good laboratory for probing ultralight particles. Previous works have set constraints on ultralight particles generated inside the Sun via stellar cooling^79,80,84 and axion decay⁸⁵. On the other hand, ref. ⁸⁶ proposed that DPDM can resonantly convert into monochromatic radio-frequency electromagnetic (EM) waves in the solar corona. This phenomenon occurs at a radius where the plasma frequency equals the DPDM mass. Furthermore, with the presence of the solar magnetic field, axion DM can also resonantly convert into radio waves in the solar corona.

In this work, we investigate such resonantly converted monochromatic radio signal within the solar observation data collected by Low Frequency Array (LOFAR) telescope⁸⁷. To calculate the signal, we carry out simulations of EM wave propagation inside solar corona of the quiet Sun. Subsequently, we compare the signal with the LOFAR data to deduce the upper limits for both the DPDM model and axion DM model. However, due to the relatively weak nature of the solar coronal magnetic field, our constraint on axion parameters does not exceed many existing constraints. Therefore, we focus on the dark photon case in the main text while leaving the detailed discussion of the axion case in the methods, subsection Constraint on axion-like particle dark matter. We set the 95% C.L. upper limit on the kinetic mixing coupling between DPDM and photon to about 10⁻¹³ within the frequency range of 30−80 MHz.

Results

Resonant conversion of ultralight DM into photons in solar plasma

For the DPDM model, dark photons interact with SM particles through kinetic mixing, and the corresponding Lagrangian can be written as

$${{{{{{{{\mathcal{L}}}}}}}}}_{{A}^{{\prime} }\gamma }=-\frac{1}{4}{F}_{\mu \nu }^{{\prime} }{F}^{{\prime} \mu \nu }+\frac{1}{2}{m}_{{A}^{{\prime} }}^{2}{A}_{\mu }^{{\prime} }{A}^{{\prime} \mu }-\frac{1}{2}\epsilon {F}_{\mu \nu }{F}^{{\prime} \mu \nu },$$

(1)

where ${A}^{{\prime} }$ and γ represent dark photon and photon respectively, F_μν and ${F}^{{\prime} \mu \nu }$ represent the field strengths of photon and dark photon respectively, with the Greek letters μ, ν denoting the vector indices, ${m}_{{A}^{{\prime} }}$ denotes the mass of dark photon, ${A}_{\mu }^{{\prime} }$ is the vector field of dark photon, and ϵ stands for the kinetic mixing parameter.

In the solar corona, the presence of free electrons gives rise to a plasma frequency denoted by ω_p, serving as the effective mass for the EM wave. This quantity is determined by the free electron density n_e in the non-relativistic plasma, and can be represented as

$$\begin{array}{r}{\omega }_{p}={\left(\frac{4\pi {\alpha }_{{{{{{{{\rm{EM}}}}}}}}}{n}_{e}}{{m}_{e}}\right)}^{\frac{1}{2}}={\left(\frac{{n}_{e}}{7.3\times 1{0}^{8}{{{{{{{{\rm{cm}}}}}}}}}^{-3}}\right)}^{\frac{1}{2}}\mu {{{{{{{\rm{eV}}}}}}}},\end{array}$$

(2)

where α_EM represents the fine structure constant and m_e is the electron mass. It is noteworthy that we employ natural units throughout our paper, thereby setting ℏ and c to unity: ℏ = c = 1. When a dark photon ${A}^{{\prime} }$ propagates in the plasma, it can resonantly convert into a SM photon when ${\omega }_{p} \, \approx \, {m}_{{A}^{{\prime} }}$⁸⁶. In the solar corona, we have n_e monotonically decreasing from 10¹⁰ to 10⁶ cm⁻³ with increasing height above the solar photosphere. Therefore, the corresponding plasma frequency scans from 4 × 10⁻⁶ to 4 × 10⁻⁸ eV. If the DM mass ${m}_{{{{{{{{{\rm{A}}}}}}}}}^{{\prime} }}$ falls within this range, the resonant conversion of DPDM into EM waves can occur at a specific radius r_c satisfying ${\omega }_{p}({r}_{c})={m}_{{{{{{{{{\rm{A}}}}}}}}}^{{\prime} }}$. The frequency of the converted EM wave, ${m}_{{{{{{{{{\rm{A}}}}}}}}}^{{\prime} }}/(2\pi )$, lies within the radio-frequency range of about 10−1000 MHz. Therefore, it can be tested by various radio telescopes engaged in solar physics programs, such as LOFAR⁸⁷ and SKA⁸⁸. Since DM in the Galactic halo is non-relativistic with the typical velocity v_DM approximately 10⁻³ times the speed of light, the converted EM wave is nearly monochromatic with a spread of about 10⁻⁶ around its central value⁸⁶. The DM dispersion bandwidth B_sig can be evaluated by

$${B}_{{{{{{{{\rm{sig}}}}}}}}} \, \approx \, \frac{{m}_{{A}^{{\prime} }}{v}_{{{{{{{{\rm{DM}}}}}}}}}^{2}}{2\pi }\, \approx \, 130\,{{{{{{{\rm{Hz}}}}}}}}\left(\frac{{m}_{{A}^{{\prime} }}}{\mu {{{{{{{\rm{eV}}}}}}}}}\right).$$

(3)

Our analysis adopts the electron density profile for the quiet Sun provided by LOFAR observations⁸⁹, shown as the solid blue line in Fig. 1, and the DM wind constantly passes through the solar atmosphere. For specific details regarding different solar density profiles in the context of the quiet Sun, refer to the methods, subsection The solar model. The probability of DPDM resonantly converting into photons is^86,90

$${P}_{{A}^{{\prime} }\to \gamma }({v}_{rc})=\frac{2}{3}\times \pi {\epsilon }^{2}{m}_{{A}^{{\prime} }}{v}_{rc}^{-1}{\left| \frac{\partial \ln {\omega }_{p}^{2}(r)}{\partial r}\right| }_{r={r}_{c}}^{-1},$$

(4)

where v_rc is the radial velocity at the resonant layer. The prefactor 2/3 arises in Eq. (4) because the longitudinal mode of photons converted from the corresponding mode of dark photon cannot propagate out of the plasma. The conversion probability Eq. (4) accounts for two transverse modes, and we assume that dark photons polarize in three directions with equal probability.

Fig. 1: Comparison between different electron density profiles.

Various density profiles are depicted using different lines. The solid blue line represents the profile derived from LOFAR observations⁸⁹. In comparison, the solid orange line represents the profile from ref. ¹⁰⁴. The dashed cyan line represents a simple hydrostatic model and the dashed red line represents an r⁻² profile⁸⁹. The gray shaded region denotes the frequency region ~30−80 MHz which our study focuses on.

Utilizing the conversion probability, the radiation power ${{{{{{{\mathcal{P}}}}}}}}$ per solid angle dΩ at the conversion layer can be derived as

$$\frac{d{{{{{{{\mathcal{P}}}}}}}}}{d\Omega }=\int\,d{{{{{{{{\bf{v}}}}}}}}}_{0}\,{f}_{{{{{{{{\rm{DM}}}}}}}}}({{{{{{{{\bf{v}}}}}}}}}_{0}){P}_{{A}^{{\prime} }\to \gamma }({v}_{0}){\rho }_{{{{{{{{\rm{DM}}}}}}}}}v({r}_{c}){r}_{c}^{2},$$

(5)

where the DM density is ρ_DM = 0.3 GeV cm^−3 91,92 and v₀, the initial DM velocity, follows a Maxwellian distribution f_DM(v₀) with the most probable velocity of 235 km/s^93,94. $v({r}_{c})=\scriptstyle\sqrt{{v}_{0}^{2}+2{G}_{N}{M}_{\odot }/{r}_{c}}$ corresponds to the DM velocity at the conversion layer including the gravitational effect of the Sun, with G_N standing for the gravitational constant and M_⊙ representing the solar mass. Detailed derivations for the conversion probability and the radiation power can be found in the methods, subsection The conversion probability of ${A}^{{\prime} }\to \gamma$ and the radiation power. We have furthermore demonstrated that electron density fluctuations do not change the result of the conversion probability in the methods, subsection Impact of small-scale fluctuations on conversion probability.

Propagation of converted photons in solar plasma

The converted EM waves propagating through the corona will experience interactions with the plasma, including both absorption and scattering processes. The absorption of these converted photons is mainly through the inverse bremsstrahlung process. Due to refraction, the converted EM wave would propagate radially outward once it exits of the resonant region if scatterings between the EM wave and the plasma were absent⁸⁶. However, the presence of scattering in the inhomogeneous plasma will randomize the direction of EM waves, leading to a broadened angular distribution of the outgoing EM waves^95,96. We are using LOFAR data made in the tied-array beam mode. While this mode offers a nice angular resolution⁸⁷, the field of view (FOV) of each LOFAR beam is significantly smaller than the total angular span of the Sun. Consequently, we expect the scattering effect to suppress the signal observed by the LOFAR detector.

When accounting for both absorption and scattering effects, the spectral flux density received by LOFAR can be expressed as

$${S}_{{{{{{{{\rm{sig}}}}}}}}}=\frac{1}{{{{{{{{\mathcal{B}}}}}}}}}\frac{1}{{d}^{2}}\frac{d{{{{{{{\mathcal{P}}}}}}}}}{d\Omega }{P}_{{{{{{{{\rm{sur}}}}}}}}}(f)\beta (f),$$

(6)

where d = 1 AU is the distance between Earth and the Sun. ${{{{{{{\mathcal{B}}}}}}}}$ represents the bandwidth, which is the larger one between the DM dispersion bandwidth B_sig which is about 130 Hz and the spectral resolution of the telescope ${B}_{{{{{{{{\rm{res}}}}}}}}}=97$ kHz. In our case, ${B}_{{{{{{{{\rm{sig}}}}}}}}}\, \ll \,{B}_{{{{{{{{\rm{res}}}}}}}}}$, so we have ${{{{{{{\mathcal{B}}}}}}}}={B}_{{{{{{{{\rm{res}}}}}}}}}$. The survival probability P_sur(f) and the factor β(f) are defined later. It is noteworthy that the energy dispersion could be enlarged by scatterings with the plasma inhomogeneities. However, this impact is negligible because the inhomogeneities can be treated as effectively static, given their velocities are much lower than the speed of light, and only elastic scatterings need to be considered⁹⁶. The speed of inhomogeneities may become important for photons with the smallest velocities just after conversion. The typical density fluctuation is the ion-sound waves⁹⁷ with the speed ${C}_{s} \, \approx \, \sqrt{\left.\right[{T}_{e}(1+3{T}_{i}/{T}_{e})/{m}_{i}}$ which is about 100 km/s (T_e, T_i, and m_i are respectively the electron temperature, ion temperature, and ion mass), which is comparable with the DM velocity v_DM which is approximately 10⁻³c in Eq. (3). This similarity implies that the line width cannot be broaden significantly. As a result, the signal line still safely locates within a single LOFAR frequency bin. Furthermore, the effect of inhomogeneities on energy dispersion diminishes quickly as the converted photons rapidly become relativistic after leaving the conversion layer and as the electron density n_e decreases. Therefore, the ray-tracing simulation of radio photon propagation⁹⁶ considers angular dispersion due to inhomogeneities while ignoring energy dispersion.

In the context of Eq. (6), the term P_sur corresponds to the survival probability of the converted photon. It is important to note that for each converted photon, P_sur also depends on the path it travels. Therefore, numerical simulations are essential for accurately calculating P_sur. The β factor in Eq. (6) parameterizes the scattering effect and is defined as

$$\beta (f)=\frac{{d}^{2}}{{R}_{S}^{2}}{\int_{{{{{{{{\rm{beam}}}}}}}}}}\frac{g({\theta }_{1},{\phi }_{1})}{{r}^{2}}dS,$$

(7)

where g(θ₁, ϕ₁) is the angular distribution function of scattered photons at the last scattering radius R_S, beyond which the scattering process can be neglected. The value of R_S(f) is determined by numerical simulation and typically ranges from about 5 to 7R_⊙, with a slight dependence on photon frequency. The integration in Eq. (7) is over the last scattering surface, and r signifies the distance from the integrated surface element dS to LOFAR. The detailed derivation and computation of Eq. (6) involve intricate but fundamental geometric analyses, and are presented in the methods, subsection The effective spectral flux density received by LOFAR stations.

For simulating the propagation of converted photons within the corona plasma, considering both absorption and scattering effects, we employ the Monte Carlo ray-tracing method developed in ref. ⁹⁶. We describe the scattering process of radio waves using the Fokker-Planck and Langevin equations based on the Hamilton equations for photons^96,97,98. In our simulation, we utilize the Kolmogorov spectrum to describe electron density fluctuations in the quiet Sun, with δn_e/n_e = 0.1, following the work of ref. ⁹⁵. Additionally, we consider the anisotropic density fluctuation magnitude as α_anis = 0.1⁹⁵. Here, α_anis represents the anisotropy parameter, which is the ratio between the perpendicular and parallel correlation lengths⁹⁶.

Then for each frequency, we calculate P_sur(f) and β(f), and simulation results are presented in Fig. 2. It is noticeable that the absorption effect becomes more prominent as the frequency increases. Similarly, the smearing effect exhibits a similar trend, primarily due to the diminishing FOV of LOFAR with increasing frequency. This reduction in FOV at higher frequencies amplifies the impact of the smearing effect on the observations.

Fig. 2: The propagation coefficients as functions of frequency.

The survival probability P_sur and the smearing factor β as functions of the photon frequency f are depicted as the solid black line and the solid purple line, respectively.

LOFAR data analysis and setting constraints on the ultralight DM couplings

LOFAR is an advanced radio interferometer with high resolution and sensitivity. The observation data of the beam-formed mode⁸⁷, where 24 LOFAR core stations in the Netherlands are combined to form 127 tied-array beams. This mode offers significantly increased frequency resolution while reducing spatial resolution. However, the 3.5 km baseline of the LOFAR core limits the FOV to only about ${5}^{{\prime} }$ at 32 MHz⁸⁷. The observation data we use is the spectral flux density calibrated in solar flux unit (sfu) within the frequency range of 30-80 MHz. Since some beams are outside the solar surface, only beams with fluxes greater than half of the maximum beam flux are selected. We have data from three different observation periods, all with an observation duration of 17 minutes, which were carried out on 25 April 2015, 3 July 2015, and 3 September 2015. The bandwidth is ${B}_{{{{{{{{\rm{res}}}}}}}}}=97$ kHz.

The data from the selected beams is averaged. The resulting averaged data is distributed across 516 frequency bins, with each bin containing 6000 time bins. To eliminate burst-like noises, a data-cleaning process is employed. Firstly, the 6000-time series is divided into 150 intervals, each encompassing 40 time bins, which is sufficient for capturing statistical behavior. The interval with the lowest mean is selected as the reference interval. Subsequently, the mean μ_t and standard deviation σ_t of each interval are compared to those of the reference interval. Intervals meeting the conditions μ_t[test] < μ_t[ref] + 2σ_t[ref] and σ_t[test] < 2σ_t[ref] are retained. This data-cleaning process only removes transient noises while preserving the time-independent ultralight DM signal.

After data cleaning, for each frequency bin, i, we can get the average value ${\bar{O}}_{i}$ and the standard deviation ${\sigma }_{{\bar{O}}_{i}}$ as the statistical uncertainty of the time series. We parameterize the background locally by fitting each frequency bin and its adjacent k bins with a polynomial function of degree n. In practice, we choose k = 10 and n = 3. Then, we use the least square method to evaluate the deviation of data to the background fit. The fitting deviation is taken to be the systematic uncertainty ${\sigma }_{i}^{sys}$. The total uncertainty is in the quadrature form, ${\sigma }_{i}^{2}={\sigma }_{{\bar{O}}_{i}}^{2}+{({\sigma }_{i}^{sys})}^{2}$. It turns out that ${\sigma }_{i}^{sys}$ always dominates in σⁱ.

We adopt the log-likelihood ratio test method⁹⁹ to set upper limits on the DPDM parameter space. We construct the likelihood function for a specific frequency bin i₀ in the Gaussian form¹⁰⁰

$$L(S,a)=\mathop{\prod }\limits_{i={i}_{0}-5}^{{i}_{0}+5}\frac{1}{{\sigma }_{i}\sqrt{2\pi }}\exp \left[-\frac{1}{2}{\left(\frac{B(a,{f}_{i})+S{\delta }_{i{i}_{0}}-{\bar{O}}_{i}}{{\sigma }_{i}}\right)}^{2}\right],$$

(8)

where B(a, f_i) is the polynomial function used for background fitting, the coefficients a = (a₁, a₂, a₃) are treated as nuisance parameters, S denotes the assumed DPDM-induced signal at bin i₀, and ${\delta }_{i{i}_{0}}$ is the Kronecker delta. We then build the following test statistic^99,100

$${q}_{S}=\left\{\begin{array}{ll}-2\ln \left[\frac{L(S,\tilde{a})}{L(\hat{S},\hat{a})}\right],\quad &\hat{S} \,\le \,S \\ 0,\hfill \quad &\hat{S} \, > \, S \end{array}\right..$$

(9)

In the denominator, the likelihood L gets maximized at $a=\hat{a}$ and $S=\hat{S}$; in the numerator, L gets maximized at $a=\tilde{a}$ for a specified S. The test statistic q_S follows the half-χ² distribution, with the probability density function

$$h({q}_{S}| S)=\frac{1}{2}\delta ({q}_{S})+\frac{1}{2}\frac{1}{\sqrt{2\pi }}\frac{1}{\sqrt{{q}_{S}}}{{{{{{{{\rm{e}}}}}}}}}^{-{q}_{S}/2},$$

(10)

the cumulative distribution function of which is given by $H({q}_{S}| S)=1/2\left(1+{{{{{{{\rm{erf}}}}}}}}\left(\sqrt{{q}_{S}/2}\right)\right)$, where ${{{{{{{\rm{erf}}}}}}}}(x)$ is the Gauss error function. Then, we can define the following criterion^99,100:

$${p}_{S}=\frac{1-{{{{{{{\rm{erf}}}}}}}}(\sqrt{{q}_{S}/2})}{1-{{{{{{{\rm{erf}}}}}}}}(\sqrt{{q}_{0}/2})},$$

(11)

which measures how far the assumed signal is away from the null result S = 0. To obtain the 95% confidence level (C.L.) upper limit ${S}_{\lim }$, we set p_S = 0.05. The results of ${S}_{\lim }$ as functions of frequency are shown in Fig. 3, with the datasets used stemming from three observation periods represented by different colors. Among the constraints from the three datasets, we select the strongest constraint at each frequency bin to determine the final upper limit. In the 112th frequency bin with f = 40.6 MHz for all three periods of observations, an increasing intensity was observed. However, this bin was identified as a bad channel and subsequently excluded from our analysis. Moreover, similar issues were identified in the 25th, 26th, and 27th bins (f = 32.2 MHz), the 34th bin (f = 33.0 MHz), and the 101st bin (f = 39.5 MHz) of the observations on 25 April 2015, in the 46th bin (f = 34.2 MHz) and the 208th bin (f = 50.0 MHz) of the observations on 3 July 2015. These bad channels were also removed from our analysis.

Fig. 3: The model-independent constraints on the monochromatic signal.

Model-independent 95% C.L. upper limits ${S}_{\lim }$ regarding photon frequency f are derived from LOFAR data on a constant monochromatic signal. The limits obtained from the observation data on 25 April 2015, 3 July 2015, and 3 September 2015 are represented by the blue, orange, and green curves, respectively.

We calculate the 95% C.L. upper limits on the kinetic mixing parameter ϵ for the DPDM model by requiring ${S}_{\lim }$ equal to S_sig in Eq. (6). The upper limit on ϵ derived from LOFAR data for DPDM is depicted in Fig. 4, which shows that the upper limit on ϵ is about 10⁻¹³ within the frequency range 30−80 MHz. It is about one order of magnitude better than the existing CMB constraint^6,53, and is complementary to other searches for DPDM at higher frequency, such as the Dark E-field experiment¹⁰¹.

Fig. 4: The constraints on the parameter space of DPDM.

95% C.L. upper limit on the kinetic mixing parameter ϵ for DPDM regarding the DPDM mass ${m}_{A{\prime} }$ from 17 minutes observation of LOFAR data is shown in the cyan shaded region. We also show the existing constraints (summarized in ref. ⁵²) from the CMB distortion (95%)^6,53, the haloscope searches WISPDMX (95%)¹²⁷, and Dark E-field experiment (5σ)¹⁰¹ in gray shaded regions. Different constraints may choose different confidence levels, and we keep their original choice unchanged as labeled in the parentheses following each experiment. The existing constraints also assume the dark matter density ρ_DM = 0.3 GeV cm⁻³, the same as our choice, and are scaled by the dark photon density ${({\rho }_{A{\prime} }/{\rho }_{{{{{{{{\rm{DM}}}}}}}}})}^{0.5}$.

Based on the same data analysis method, we can set upper limits on the axion-photon coupling for the case of axion DM. However, due to the relatively weak solar coronal magnetic field, our resulting constraint for the axion case is not as strong as many existing constraints. This portion of our analysis is detailed in the methods, subsection Constraint on axion-like particle dark matter.

Discussion

When DPDM or axion DM pass across the Sun, they can resonantly convert into EM waves in the solar corona. To explore this phenomenon, we conducted numerical simulations of the converted photons propagating in the plasma, including the effects of absorption and scattering. Radio telescopes for solar observations are capable of detecting the monochromatic converted EM waves. We used three datasets of 17-minute observation data from LOFAR to search for such signals. We found that this method sets a stringent limit on the kinetic mixing parameter for dark photons, specifically ϵ at approximately 10⁻¹³, within the frequency range 30−80 MHz. This limit is about one order of magnitude stronger than the constraint derived from CMB observations. Similarly, we obtain an upper limit on the axion-photon coupling g_aγγ for the axion DM model in the same frequency range. The constraint on g_aγγ is better than that from Light-Shining-through-a-Wall experiments but is not comparable with the CAST and astrophysical bounds. The LOFAR data analysis in this work shows great potential in searching for ultralight DM with radio telescopes. With greater sensitivity, we expect future radio programs such as the SKA telescope are expected to yield even greater sensitivity in the search for DPDM and axion DM. Terrestrial radio telescopes cannot search for DPDM with frequencies lower than 10 MHz due to the screening effect from the ionosphere. In these cases, the use of solar probes, such as the STEREO¹⁰² satellite and the Parker Solar Probe¹⁰³, equipped with radio spectrometers, could offer an avenue for DPDM detection.

Methods

The solar model

In our study, we centered our attention on the quiet Sun due to its reduced occurrence of active events like turbulence and flares. To conduct our calculations, we utilized the electron number density (n_e) profile derived from LOFAR observations⁸⁹, which employed ray-tracing simulations to fit the solar intensity profile observed by LOFAR in the frequency range of 30−80 MHz.

There have been other density profiles for the quiet Sun, but their differences are within factor of a few. For example, the density profiles based on the work of V. De La Luz, et al.¹⁰⁴, are derived from the temperature (T) and hydrogen density (n_H) profiles for the quiet Sun, based on the photospheric model from ref. ¹⁰⁵ and the coronal model from refs. ^106,107. The consistency and validation of these profiles have been confirmed by various research groups¹⁰⁸ using a chromosphere model from ref. ¹⁰⁹ and the coronal model from refs. ^106,107. These independent calculations consistently agree with each other and have been validated by observations of atomic lines in the soft X-ray range¹¹⁰ and extreme ultraviolet range¹⁰⁵.

Furthermore, one can adopt a spherically symmetric and hydrostatic model for the quiet Sun, where gas pressure and gravitational force remain in equilibrium, resulting in a static configuration over time. The hydrostatic equilibrium in the quiet Sun region has been confirmed in previous studies^105,110. Here we provide a simple analytical expression to parameterize the hydrostatic density model, which can be expressed in an exponential form. In this simplified form, the electron density is modeled as⁸⁹

$$\begin{array}{r}{n}_{e}={N}_{0}\exp (1/({H}_{0}r)),\end{array}$$

(12)

with the parameters N₀ and H₀, and the latter is defined as

$${H}_{0}=\frac{{k}_{B}T}{0.6{m}_{p}{g}_{\odot }}\frac{1}{{R}_{\odot }^{2}},$$

(13)

where R_⊙ is the solar radius, k_B is the Boltzmann constant, g_⊙ = 274 m s⁻² signifies the gravitational acceleration at the coronal base, and 0.6 times the proton mass m_p gives the average particle mass in the corona¹¹¹. The temperature T corresponds to a scale height temperature determined by both electron and ion temperatures. There are two parameters to be determined: the density N₀ and the temperature T. To carry out our calculation, N₀ = 1.6 × 10¹¹ m⁻³ and T = 2 × 10⁶ K are used from LOFAR observation fit⁸⁹.

In Fig. 1, we present a comparison between the profile we adopted from LOFAR observations⁸⁹, the profile from ref. ¹⁰⁴, the hydrostatic profile modeled by Eq. (12), and the r⁻² profile. The r⁻² density profile indicates that a constant solar wind speed has been attained⁸⁹. It can be seen that the solar profile from LOFAR observations closely matches the hydrostatic profile at the high-frequency range, but exhibits a slower decline and transitions into the r⁻² profile at the low-frequency range. This behavior is expected, as it signifies the shift from subsonic plasma flow in the corona to the supersonic solar wind⁸⁹.

As a result, the variation in the electron number density (n_e) profile for the quiet Sun across different observations remains within a factor of a few. Although this variation does affect the plasma frequency, it is proportional to the square root of the electron number density, and it only shifts the location of the resonant region. Additionally, the derivative of n_e with respect to radius plays a role in determining the conversion probability, yet its effect is also relatively minor. These uncertainties are small and have a negligible effect on the resulting photon signal.

The conversion probability of ${A}^{{\prime} }\to \gamma$ and the radiation power

The conversion probability for DPDM to photon can be calculated either by quantum field theory (QFT) as a 1 → 1 process, or by solving linearized wave equation^86,90,112. Here we take linearized wave method as an example, providing formulas for conversion probability and radiation power. It is important to note that in this subsection, our formulas are derived from the solar profile without accounting for small-scale fluctuations. In the following subsection, we will provide estimations of the impact of these small-scale fluctuations.

We can eliminate the kinetic mixing term by performing a rotation of the vector fields in Eq. (1) into interaction basis: ${A}_{\mu }\to 1/\sqrt{1-{\epsilon }^{2}}{A}_{\mu },{A}_{\mu }^{{\prime} }\to -\epsilon /\sqrt{1-{\epsilon }^{2}}{A}_{\mu }+{A}_{\mu }^{{\prime} }$, where A_μ represents the vector field of photon. In this basis, the equations of motion become

$$\left[\frac{{\partial }^{2}}{\partial {t}^{2}}-\frac{{\partial }^{2}}{\partial {r}^{2}}+\left(\begin{array}{cc}{\omega }_{p}^{2}&-\epsilon {m}_{{A}^{{\prime} }}^{2}\\ -\epsilon {m}_{{A}^{{\prime} }}^{2}&{m}_{{A}^{{\prime} }}^{2}\end{array}\right)\right]\left(\begin{array}{c}A(r,t)\\ {A}^{{\prime} }(r,t)\end{array}\right)=0,$$

(14)

which are coupled wave equations.

These second-order coupled equation can be approximated to first-order linearized wave equations using the WKB approximation, as the spatial variation of the plasma frequency occurs on a much larger scale than the wavelength of DPDM. Consequently, we have ${\partial }_{t}^{2}-{\partial }_{r}^{2}=-{\omega }^{2}-{\partial }_{r}^{2} \, \approx -2{k}_{r}({k}_{r}+i{\partial }_{r})-{m}_{{A}^{{\prime} }}^{2}-{k}_{T}^{2}$ under a plane wave solution $A(r,t)=\tilde{A}(r)\exp (-i\omega t+i{k}_{r}r)$ with frequency ω and wavenumber $k=\sqrt{{\omega }^{2}-{m}_{{A}^{{\prime} }}^{2}}$, where k_r and k_T is the longitudinal and transverse components of momentum k. The resulting first-order linearized wave equation can be expressed as

$$\begin{array}{r}(i{\partial }_{r}-{H}_{0}-{H}_{I})\left(\begin{array}{c}\tilde{A}(r)\\ \tilde{{A}^{{\prime} }}(r)\end{array}\right)=0,\end{array}$$

(15)

where

$${H}_{0}= \frac{1}{2{k}_{r}}\left(\begin{array}{cc}{\omega }_{p}^{2}-{m}_{{A}^{{\prime} }}^{2}-{k}_{T}^{2}&0\\ 0&-{k}_{T}^{2}\end{array}\right),\\ {H}_{I}= \frac{1}{2{k}_{r}}\left(\begin{array}{cc}0&-\epsilon {m}_{{A}^{{\prime} }}^{2}\\ -\epsilon {m}_{{A}^{{\prime} }}^{2}&0.\end{array}\right).$$

(16)

This equation can be solved perturbatively by expanding the time-evolution operator in Dyson series⁹⁰. At the first-order, the conversion probability is given by¹¹²

$${P}_{{A}^{{\prime} }\to \gamma }={\left| \int\nolimits_{{r}_{0}}^{r}d{r}^{{\prime} }\frac{-\epsilon {m}_{{A}^{{\prime} }}^{2}}{2{k}_{r}}{e}^{i\int\nolimits_{{r}_{0}}^{{r}^{{\prime} }}d{r}^{{\prime\prime} }\frac{1}{2{k}_{r}}\left[{\omega }_{p}{({r}^{{\prime\prime} })}^{2}-{m}_{{A}^{{\prime} }}^{2}\right]}\right| }^{2}.$$

(17)

This formula can be further simplified to Eq. (4) by using saddle point approximation

$$\int\nolimits_{-\infty }^{\infty }dr{e}^{-f(r)} \, \approx \, {e}^{-f({r}_{0})}\sqrt{\frac{2\pi }{{f}^{{\prime\prime} }({r}_{0})}},$$

(18)

where ${f}^{{\prime} }({r}_{0})=0$. The thickness of resonant layer is on the order of 10³ km for the frequency range under consideration. The WKB approximation and the saddle point approximation can be tested even in the presence with small-scale fluctuations, as we will demonstrate in the next subsection. Additionally, we will numerically show that the value of the conversion probability remains unaffected by small-scale fluctuations in the upcoming subsection.

The radiation power can be obtained from the conversion probability. Taking into account the gravitational focus effect and considering incoming DPDM at infinity moving under the influence of the gravitational potential, the radiation power per solid angle is

$$\frac{d{{{{{{{\mathcal{P}}}}}}}}}{d\Omega }\, \approx \, 2\frac{1}{4\pi }{\rho }_{{{{{{{{\rm{DM}}}}}}}}}\int\,d{{{{{{{{\bf{v}}}}}}}}}_{0}\,{f}_{{{{{{{{\rm{DM}}}}}}}}}{v}_{0}\int\nolimits_{0}^{b}dz2\pi z{P}_{{A}^{{\prime} }\to \gamma }({v}_{r})\\= {\rho }_{{{{{{{{\rm{DM}}}}}}}}}\int\,d{{{{{{{{\bf{v}}}}}}}}}_{0}\,{f}_{{{{{{{{\rm{DM}}}}}}}}}{P}_{{A}^{{\prime} }\to \gamma }({v}_{0})v({r}_{c}){r}_{c}^{2},$$

(19)

where z is the impact parameter for DPDM, b = r_cv(r_c)/v₀ is the maximum impact parameter for ${A}^{{\prime} }$ to reach the conversion layer at $r={r}_{c},v({r}_{c})=\sqrt{{v}_{0}^{2}+2{G}_{{{{{{{{\rm{N}}}}}}}}}{M}_{\odot }/{r}_{c}}$ is the velocity of DPDM at r_c, and the radial direction velocity of DPDM at r_c with different impact parameter is ${v}_{r}(z)=\sqrt{2{G}_{{{{{{{{\rm{N}}}}}}}}}{M}_{\odot }/{r}_{c}+{v}_{0}^{2}-{v}_{0}^{2}{z}^{2}/{r}_{c}^{2}}$. The factor of 2 accounts for both incoming and outgoing DPDM as incoming DPDM will be totally reflected.

Impact of small-scale fluctuations on conversion probability

In this subsection, we will estimate the influence arising from small-scale inhomogeneities in the plasma by incorporating density fluctuations.

Density fluctuation can lead to three main effects: (1) modifying the magnitude of the ${A}^{{\prime} }\to \gamma$ conversion probability by altering $\left\vert \overrightarrow{\nabla }{\omega }_{p}\right\vert$; (2) introducing non-spherical modifications to the conversion surface; and (3) introducing scattering and absorption of the converted photons, resulting in smearing of their velocity directions and a reduction in photon flux. The third effect has been addressed in the Results section when accounting for the propagation effects, utilizing the Monte Carlo ray-tracing simulations.

Regarding the first effect, the inclusion of density fluctuations introduces two opposing influences that modify the conversion probability, ${P}_{{A}^{{\prime} }\to \gamma }$. On one hand, the derivative of electron density with respective to distance becomes larger due to the fluctuations, leading to a decrease in ${P}_{{A}^{{\prime} }\to \gamma }$. On the other hand, more resonant points where ${\omega }_{p}(r)={m}_{{A}^{{\prime} }}$ are introduced, which increases ${P}_{{A}^{{\prime} }\to \gamma }$. It turns out that these two influences cancel each other out, resulting in ${P}_{{A}^{{\prime} }\to \gamma }$ with density fluctuations remaining the same as the original value. In addition, the non-spherical effect in the second effect is insignificant. In the following, we will quantify the first and second effects.

In ref. ⁹⁵, an advanced Monte Carlo simulation technique was employed to address density fluctuations and their impact on photon refraction and scattering during plasma propagation. Their findings suggested that refraction and scattering might be the primary factors contributing to the observed lower brightness temperatures in quiet-Sun radio emissions across different frequencies, deviating from expected values. Hence, we adopt their mathematical framework for describing density fluctuations and incorporate it into our own research.

First and foremost, we emphasize that the density fluctuations in the plasma density are relatively small, with an approximate magnitude of⁹⁵

$${\epsilon }_{e}\equiv \Delta {n}_{e}/{n}_{e} \, \approx \, 10\%,$$

(20)

and importantly, this fluctuation fraction remains constant as the radial distance changes^95,113.

The probability distribution of plasma density fluctuations is described by the spatial power spectrum. For the solar corona of the quiet Sun, the spatial power spectrum of density fluctuations can be expressed as^114,115

$$P(q)={C}_{N}^{2}{q}^{-\alpha },\quad {q}_{o} \, < \, q \, < \, {q}_{i},$$

(21)

where ${C}_{N}^{2}$ is the structural constant, q represents the spatial wavenumber, and α corresponds to the power-law exponent, which is chosen as α = 11/3 to reflect the Kolmogorov spectrum.

The scale of density turbulence is defined as l ≡ 2π/q, with l_i = 2π/q_i and l_o = 2π/q_o denoting the inner and outer scales of the density turbulence, respectively. It is reasonable to assume that l_o ≈ 10⁶l_i⁹⁵. Consequently, the steep shape of the Kolmogorov-type spectrum for P(q) indicates that density fluctuations predominantly occur on larger scales.

The inner scale, denoted as l_i, can be associated with the ion inertial scale, given by the expression

$${l}_{i}=\frac{684}{\sqrt{{n}_{e}/{{{{{{{{\rm{cm}}}}}}}}}^{-3}}}\,{{{{{{{\rm{km}}}}}}}}.$$

(22)

Consequently, for the plasma layers corresponding to frequencies of 30, 40, 60, and 80 MHz, the respective inner scales l_i are estimated to be 0.2 km, 0.15 km, 0.1 km, and 0.075 km.

The spatial power spectrum can be normalized to the variance of the density fluctuations $\langle \Delta {n}_{e}^{2}\rangle \equiv {({\epsilon }_{e}{n}_{e})}^{2}$ as,

$$\int\nolimits_{{q}_{o}}^{{q}_{i}}P(q)4\pi {q}^{2}dq=\langle \Delta {n}_{e}^{2}\rangle={({\epsilon }_{e}{n}_{e})}^{2}.$$

(23)

Using the Kolmogorov spectrum with α = 11/3, it can be determined that

$${C}_{N}^{2}=\frac{{q}_{o}^{\alpha -3}}{6\pi }\langle \Delta {n}_{e}^{2}\rangle,$$

(24)

$$P(q)=\frac{{q}_{o}^{\alpha -3}}{6\pi }{q}^{-\alpha }\langle \Delta {n}_{e}^{2}\rangle .$$

(25)

The fluctuations can be expressed in the Fourier modes,

$$\Delta {n}_{e}(r)=\int\nolimits_{1}^{{\tilde{q}}_{i}}d\tilde{q}\,\Delta {\tilde{n}}_{e}(q){{{{{{{{\rm{e}}}}}}}}}^{iqr}.$$

(26)

We have rescaled the momentum $\tilde{q}\equiv q/{q}_{o}$ for convenience. Subsequently, the fluctuations averaged over the length scale l_o = 2π/q_o are

$$\left\langle \Delta {n}_{e}^{2}\right\rangle= \frac{1}{{l}_{o}}\int\,dr\int\nolimits_{1}^{{\tilde{q}}_{i}}d\tilde{q}\int\nolimits_{1}^{{\tilde{q}}_{i}}d{\tilde{q}}^{{\prime} }\,\Delta {\tilde{n}}_{e}(q)\Delta {\tilde{n}}_{e}({q}^{{\prime} }){{{{{{{{\rm{e}}}}}}}}}^{i(q-{q}^{{\prime} })r}\\= \int\nolimits_{1}^{{\tilde{q}}_{i}}d\tilde{q}\,\langle \Delta {\tilde{n}}_{e}^{2}(q)\rangle .$$

(27)

Compared with Eq. (23), we have the fluctuations in the momentum space,

$$\langle \Delta {\tilde{n}}_{e}^{2}(q)\rangle=P(q)4\pi {q}^{2}{q}_{o}=\frac{2}{3}{\tilde{q}}^{2-\alpha }\langle \Delta {n}_{e}^{2}\rangle .$$

(28)

The averaged derivative of n_e(r) with respect to r, in the squared form, is then

$$\langle {({n}_{e}^{{\prime} })}^{2}\rangle \simeq \, \langle {(\Delta {n}_{e}^{{\prime} })}^{2}\rangle \\= \frac{1}{{l}_{o}}\int\,dr\int\nolimits_{1}^{{\tilde{q}}_{i}}d\tilde{q}\int\nolimits_{1}^{{\tilde{q}}_{i}}d{\tilde{q}}^{{\prime} }\,\Delta {\tilde{n}}_{e}(q)\Delta {\tilde{n}}_{e}({q}^{{\prime} })q{q}^{{\prime} }{{{{{{{{\rm{e}}}}}}}}}^{i(q-{q}^{{\prime} })r}\\= \int\nolimits_{1}^{{\tilde{q}}_{i}}d\tilde{q}\,\langle \Delta {\tilde{n}}_{e}^{2}(q)\rangle \cdot {q}^{2}\simeq \frac{2}{3}\langle \Delta {n}_{e}^{2}\rangle \frac{1}{5-\alpha }{q}_{o}^{2}{\tilde{q}}_{i}^{5-\alpha }.$$

(29)

In the first step, the derivative of the background electron density, n_e,bkg(r), has been omitted due to its relatively small magnitude compared to that of the fluctuations. Similarly, the averaged second derivative of n_e(r) with respect to r, in the squared form, is

$$\langle {({n}_{e}^{{\prime\prime} })}^{2}\rangle \simeq \, \langle {(\Delta {n}_{e}^{{\prime\prime} })}^{2}\rangle \\= \frac{1}{{l}_{o}}\int\,dr\int\nolimits_{1}^{{\tilde{q}}_{i}}d\tilde{q}\int\nolimits_{1}^{{\tilde{q}}_{i}}d{\tilde{q}}^{{\prime} }\,\Delta {\tilde{n}}_{e}(q)\Delta {\tilde{n}}_{e}({q}^{{\prime} }){q}^{2}{q}^{{\prime} 2}{{{{{{{{\rm{e}}}}}}}}}^{i(q-{q}^{{\prime} })r}\\= \int\nolimits_{1}^{{\tilde{q}}_{i}}d\tilde{q}\,\langle \Delta {\tilde{n}}_{e}^{2}(q)\rangle \cdot {q}^{4}\simeq \frac{2}{3}\langle \Delta {n}_{e}^{2}\rangle \frac{1}{7-\alpha }{q}_{o}^{4}{\tilde{q}}_{i}^{7-\alpha }.$$

(30)

Next, we are going to examine whether the WKB approximation and saddle-point approximation are threatened by the inclusion of density fluctuations.

Using Eq. (29), we can estimate the typical length scale of density variations as

$$\delta {l}_{e}={\left| \frac{{n}_{e}^{{\prime} }}{{n}_{e}}\right| }^{-1} \simeq {\left[\sqrt{\frac{2}{3}}{\left(\frac{1}{5-\alpha }\right)}^{\frac{1}{2}}{\epsilon }_{e}{q}_{o}{\tilde{q}}_{i}^{\frac{5-\alpha }{2}}\right]}^{-1} \simeq \, 1{0}^{-3}{q}_{o}^{-1},$$

(31)

which turns out to be much larger than the dark photon wavelength, $\delta {l}_{e}{k}_{{A}^{{\prime} }} \,\approx \, 30 \, \gg \, 1$. Therefore, the WKB approximation applied in deriving Eq. (4) remains justified even with the density fluctuations included.

Another important length scale is the resonant conversion length, $\delta {l}_{{{{{{{{\rm{res}}}}}}}}}=\sqrt{2\pi /{F}^{{\prime\prime} }({r}_{c})}$. It is defined as the length along which the phase factor in Eq. (17), $F(r)\equiv \int\,dr[{\omega }_{p}^{2}(r)-{m}_{{A}^{{\prime} }}^{2}]/(2{k}_{{A}^{{\prime} }})$, changes by π. This is the length interval which dominantly contributes to ${P}_{\gamma \to {A}^{{\prime} }}$. We have

$$\delta {l}_{{{{{{{{\rm{res}}}}}}}}}=\sqrt{2\pi \frac{{k}_{{A}^{{\prime} }}}{{\omega }_{p}}\frac{1}{{\omega }_{p}^{{\prime} }}}\simeq \sqrt{2\pi }{v}_{{{{{{{{\rm{DM}}}}}}}}}{\left(\frac{\delta {l}_{e}}{{k}_{{A}^{{\prime} }}}\right)}^{1/2},$$

(32)

which obviously satisfies $\delta {l}_{{{{{{{{\rm{res}}}}}}}}}\ll\, \delta {l}_{e}$.

Next, we evaluate the robustness of the saddle-point approximation. The crucial criterion is that the second derivative $F^{\prime} (r)$ plays a dominant role in the Taylor series of F(r) compared with the higher derivative terms (note that at the resonant point, $F^{\prime} (r)=0$). Then we calculate the following quantity with the help of Eqs. (29) and (30),

$$\frac{\frac{1}{2!}{F}^{{\prime\prime} }(r)}{\frac{1}{3!}\delta {l}_{{{{{{{{{{\rm{res}}}}}}}}}}}{F}^{{\prime\prime} {\prime} }(r)} \simeq \frac{3}{\sqrt{2\pi }}\frac{{[{F}^{{\prime\prime} }(r)]}^{3/2}}{{F}^{{\prime\prime} {\prime} }(r)}\\ \simeq \frac{3}{\sqrt{2\pi }}{\left(\frac{1}{2{k}_{A^{\prime} }}\frac{{\omega }_{p}^{2}}{{n}_{e}}\right)}^{1/2}\frac{{[n^{\prime} (r)]}^{3/2}}{{n}_{e}^{{\prime\prime} }(r)}\simeq 5.$$

(33)

We see that the second derivative indeed plays a dominant role, indicating that the saddle-point approximation still holds true with an acceptable accuracy. Also, we notice that ref. ¹¹⁶ numerically shows that the saddle-point approximation works well when the ratio in Eq. (33) is larger than unity.

The above arguments show that the form of conversion probability, Eq. (4), is still correct when considering the density fluctuations. However, its numerical value may be altered by the inclusion of density fluctuations. But as we will demonstrate below, the value of ${P}_{\gamma \to A{\prime} }$ remains unchanged. As stated before, there are two new effects that counteract each other in modifying the value of the probability ${P}_{A{\prime} \to \gamma }$: the larger derivative of n_e with respect to distance and more resonant points (>1). An example of an n_e(r) profile with density fluctuations is shown in Fig. 5 where the effects of a larger derivative and more intersections can be seen. We use ${r}_{{{{{{{{\rm{dn}}}}}}}}}$ to denote the ratio between ${P}_{{A^{\prime}} \to \gamma }$ with and without density fluctuations, and it can be calculated as

$${r}_{{{{{{{{\rm{dn}}}}}}}}}= \frac{{\sum }_{{n}_{e}(r^{\prime} )={n}_{e}({r}_{c})}{\left| \frac{1}{{n}_{e}(r)}\frac{d{n}_{e}(r)}{dr}\right| }_{r=r^{\prime} }^{-1}}{{\left| \frac{1}{{n}_{e,{{{{{{{\rm{bkg}}}}}}}}}(r)}\frac{d{n}_{e,{{{{{{{\rm{bkg}}}}}}}}}(r)}{dr}\right| }_{r={r}_{c}}^{-1}} \\= \frac{{\sum }_{{n}_{e}(r^{\prime} )={n}_{e}({r}_{c})}{\left| \frac{dr}{d{n}_{e}(r)}\right| }_{r=r^{\prime} }}{{\left| \frac{dr}{d{n}_{e,{{{{{{{\rm{bkg}}}}}}}}}(r)}\right| }_{r={r}_{c}}}.$$

(34)

We can provide a rough estimate of ${r}_{{{{{{{{\rm{dn}}}}}}}}}$. Suppose the fluctuation amplitude is δn_e in a length scale δr. The intersections can only occur within the interval Δr along which the background density n_e,bkg(r) changes by δn_e. The number of intersections can be estimated as Δr/δr. Then, we have ${r}_{{{{{{{{\rm{dn}}}}}}}}} \, \approx \, (\Delta r/\delta r)\cdot (\delta r/\delta {n}_{e})/(\Delta r/\delta {n}_{e})$ which is about 1. This suggests that the two effects cancel each other out. Thus, we anticipate that the average conversion probability, when density fluctuations are considered, will remain the same as our original value.

Fig. 5: The electron density profile exhibiting fluctuations.

The electron density n_e with fluctuations is shown in blue solid line. This is plotted with the first 12 modes included. This profile is centered around the resonant layer corresponding to 40 MHz. The electron density n_e(r_c) for the 40 MHz frequency is illustrated by the orange line. The two vertical green solid line denotes the interval Δr where the intersections can occur.

We then proceed to numerically compute the ratio ${r}_{{{{{{{{\rm{dn}}}}}}}}}$ in Eq. (34) using a large sample of n_e profile generated by Monte Carlo method. As we will see below, the concise result ${r}_{{{{{{{{\rm{dn}}}}}}}}}\simeq 1$ is indeed verified.

For the convenience of numerical computation, we first need to discretize the density fluctuations, Eq. (23), as

$$\langle \Delta {n}_{e}^{2}\rangle= \frac{2}{3}\langle \Delta {n}_{e}^{2}\rangle \int\nolimits_{0}^{{\log }_{10}{\tilde{q}}_{i}}{\tilde{q}}^{3-\alpha }\ln (10)\cdot d{\log }_{10}\tilde{q}\\= \frac{2}{3}\langle \Delta {n}_{e}^{2}\rangle \mathop{\sum }\limits_{n=0}^{N}{\tilde{q}}^{3-\alpha }\cdot \Delta \cdot \ln (10)\\= \frac{2\ln (10)}{3}\mathop{\sum }\limits_{n=0}^{N}{\langle \Delta {n}_{e}^{2}\rangle }^{({{{{{{{\rm{disc}}}}}}}}.)}({q}_{n})\cdot \Delta,$$

(35)

where

$$\begin{array}{rcl}&&{\langle \Delta {n}_{e}^{2}\rangle }^{({{{{{{{\rm{disc}}}}}}}}.)}({q}_{n})\equiv \langle \Delta {n}_{e}^{2}\rangle {\tilde{q}}_{n}^{3-\alpha },\\ &&{q}_{n}=1{0}^{n\Delta }{q}_{o},\quad {\tilde{q}}_{n}=1{0}^{n\Delta }.\end{array}$$

(36)

${\langle \Delta {n}_{e}^{2}\rangle }^{({{{{{{{\rm{disc}}}}}}}}.)}({q}_{n})$ is the variance (squared) of density fluctuations of the q_n mode in the interval [${\log }_{10}({q}_{n}/{q}_{o}),{\log }_{10}({q}_{n}/{q}_{o})+\Delta$]. The discretization is carried out in the logarithmic scale, as the momentum span is broad, spanning 6 orders of magnitude from q_o to q_i. The total number of modes is $N={\log }_{10}({q}_{i}/{q}_{o})/\Delta$. Next, the variance of density fluctuations for different momentum modes can be estimated as

$$\begin{array}{r}{\sigma }_{{n}_{e}}({q}_{n})\simeq \sqrt{{\langle \Delta {n}_{e}^{2}\rangle }^{({{{{{{{\rm{disc}}}}}}}}.)}({q}_{n})}\simeq {\epsilon }_{e}{n}_{e}{\left(\frac{{q}_{n}}{{q}_{o}}\right)}^{\frac{3-\alpha }{2}}.\end{array}$$

(37)

The density as a function of distance, n_e(r), with the fluctuations taken into consideration, can be modeled as

$${n}_{e}(r)=\frac{{r}_{c}^{2}}{{r}^{2}}\left({n}_{e,{{{{{{{\rm{bkg}}}}}}}}}({r}_{c})+\mathop{\sum }\limits_{n=0}^{N}\delta {n}_{e}({q}_{n})\Delta \cdot \sin \left[{q}_{n}(r-{r}_{c})+{\phi }_{n}\right]\right).$$

(38)

Note that we have used sine functions for simplicity in numerical evaluation. We also add random phases ϕ_n for each mode. Based on (37), we have the variance of δn_e(k_n)/n_e(r_c),

$$\frac{{\sigma }_{{n}_{e}}({k}_{n})}{{n}_{e}({r}_{c})}\simeq {\epsilon }_{e}\cdot 1{0}^{n\Delta \cdot (\frac{3-\alpha }{2})}.$$

(39)

Then, we employ Monte Carlo method to generate the values for δn_e(k_n) following a Gaussian distribution with a mean value of zero and a variance of ${\sigma }_{{n}_{e}}({k}_{n})$ for each k mode. Additionally, the phase ϕ_n randomly picks up a value between 0 and 2π for each k mode. To check the effect of the fluctuations on the conversion probability, we take the frequency 40 MHz as an example. At this frequency, the solar wind dominates so that the background density profile can be taken as r⁻² as shown in Eq. (38). In this example, the corresponding resonant conversion layer is at r_c ≃ 9.4 × 10⁵ km from the solar center, which corresponds to q₀r_c ≃ 40.

To proceed, we take Δ = 0.2 and thus we have N = 30 modes in total. However, computing ${r}_{{{{{{{{\rm{dn}}}}}}}}}$ in the presence of the full N = 30 modes turns out to be challenging due to numerical limitations. Consequently, we perform the calculations for subsets of modes, specifically considering the first 2, 4, 6, . . . , 20 modes individually. For each chosen number of modes, we iteratively evaluate ${r}_{{{{{{{{\rm{dn}}}}}}}}}$ 1000 times using different Monte Carlo n_e profiles and then take the average to ensure statistical stability. In Fig. 6, we present the result of ${r}_{{{{{{{{\rm{dn}}}}}}}}}$ with more modes gradually included in computations. We see that the average value of ${r}_{{{{{{{{\rm{dn}}}}}}}}}$ converges to approximately 1, insensitive to the number of modes. Therefore, we conclude that including density fluctuations does not significantly change the value of ${P}_{A^{\prime} \to \gamma }$, as the two effects of larger derivatives and more resonant points cancel out each other.

Fig. 6: The ratio between the conversion probabilities with and without density fluctuations.

The ratio ${r}_{{{{{{{{\rm{dn}}}}}}}}}$ is computed numerically for various numbers of modes k at the 40 MHz frequency, and is shown as the blue dots, while the orange line marks the position of unity as a reference. These calculations are performed over a total of 1000 samples and the resulting values are averaged.

Next, we check the second effect of density fluctuations, which concerns the modification to the shape of the conversion surface. The deformation of the conversion surface is within the length interval Δr around r_c. As illustrated in Fig. 6, Δr/r_c ≪ 1. Therefore, the deformation effect is negligible compared with the orignal conversion sphere located at r = r_c without density fluctuations included.

In summary, we have provided a quantitative demonstration that inhomogeneities have a minimal impact on our calculations. This conclusion holds true for the condition of deriving the conversion probability, the magnitude of the conversion probability, and the deformation of the conversion sphere. There are two key factors contribute to the result: Firstly, the fraction of density fluctuation remains small, at approximately 10%. Secondly, the density fluctuation predominantly occurs at larger scales, indicating that small-scale turbulence has a limited effect.

The effective spectral flux density received by LOFAR stations

Firstly, the Field of View (FOV) of LOFAR, or effectively, the Full Width Half Maximum (FWHM) of LOFAR, is determined by

$${{{{{{{\rm{FWHM}}}}}}}}=\eta \times \frac{\lambda }{D},$$

(40)

where λ is the observation wavelength, the coefficient η = 1.02¹¹⁷, and D ≃ 3.5 km is the station diameter according to ref. ¹¹⁸. Therefore, the FWHM (for one beam) is approximately 10⁻³ rad.

We can effectively define the last scattering sphere of radius R_S, beyond which the scattering effect can be ignored, allowing the radio waves to propagate in straight lines for r > R_S. The total radiation power for dark photon signal at frequency f after conversion is $d{{{{{{{\mathcal{P}}}}}}}}/d\Omega \times 4\pi$. Therefore, the survived power at the last scattering sphere is given by

$${{{{{{{\mathcal{P}}}}}}}}={P}_{{{{{{{{\rm{sur}}}}}}}}}(f)4\pi \frac{d{{{{{{{\mathcal{P}}}}}}}}}{d\Omega }.$$

(41)

Considering a virtual source point P₁ situated within a surface element dA₁ on the last scattering sphere (as depicted in the schematic diagram of Fig. 7), the power it radiates in the direction r is

$$\begin{array}{r}d{{{{{{{{\mathcal{P}}}}}}}}}^{{\prime} }={{{{{{{\mathcal{P}}}}}}}}\frac{d{A}_{1}}{4\pi {R}_{S}^{2}}g({\theta }_{1},{\phi }_{1})d{\Omega }_{1},\end{array}$$

(42)

where the angular distribution function g(θ₁, ϕ₁) accounts for the fact that after multiple random scattering events, the radiation from the surface element is not simply in the radial direction. g(θ₁, ϕ₁) is normalized as

$$1=\int\,g({\theta }_{1},{\phi }_{1})d{\Omega }_{1}.$$

(43)

The relation $d{\Omega }_{1}=d{A}_{2}\cos {\theta }_{2}/{r}^{2}$ is useful where the cosine factor accounts for converting the receiving area dA₂ to the projected area normal to r. Then, Eq. (42) becomes

$$\begin{array}{r}d{{{{{{{{\mathcal{P}}}}}}}}}^{{\prime} }={{{{{{{\mathcal{P}}}}}}}}\frac{d{A}_{1}}{4\pi {R}_{S}^{2}}g({\theta }_{1},{\phi }_{1})\frac{d{A}_{2}\cos {\theta }_{2}}{{r}^{2}},\end{array}$$

(44)

where r is the distance from the surface element to the Earth. Meanwhile, since θ₂ is on the order of 10⁻³ rad, it follows that $\cos {\theta }_{2}\simeq 1$.

Fig. 7: Schematic diagram of the propagation of photons after the last scattering.

R_C denotes the conversion layer, and R_S denotes the last scattering sphere. A surface element dA₁, which containing a point P₁, acts as the radiation source on the last scattering sphere. θ is the polar angle of P₁. Another surface element dA₂, encompassing P₂, serves as the detection area on the Earth, which defines a solid angle dΩ₁ about P₁ in the direction of r. θ_i is the angle between the propagation vector r and the normal vector N_i of dA_i. The direction of N₂ is aligned with the line connecting the centers of the Sun and the Earth. d is the distance from the Earth to the Sun.

By substituting Eq. (41) into Eq. (44) and integrating over the area on the last scattering sphere covered by the beams, the effective spectral flux density (power per unit area and unit frequency) received by LOFAR is derived as:

$${S}_{{{{{{{{\rm{sig}}}}}}}}}={P}_{{{{{{{{\rm{sur}}}}}}}}}\frac{1}{{{{{{{{\mathcal{B}}}}}}}}}\frac{1}{{R}_{S}^{2}}\frac{d{{{{{{{\mathcal{P}}}}}}}}}{d\Omega }{\int}_{{{{{{{{\rm{beam}}}}}}}}}\frac{g({\theta }_{1},{\phi }_{1})}{{r}^{2}}d{A}_{1}.$$

(45)

As discussed in the main text, the angular distribution function g(θ₁, ϕ₁) can be determined by numerical simulations. The integration is performed in the spherical coordinates (θ, ϕ) with the Solar center as the origin. Consequently, it can be transformed into

$${S}_{{{{{{{{\rm{sig}}}}}}}}}={P}_{{{{{{{{\rm{sur}}}}}}}}}\frac{1}{{d}^{2}}\frac{1}{{{{{{{{\mathcal{B}}}}}}}}}\frac{d{{{{{{{\mathcal{P}}}}}}}}}{d\Omega }{\int}_{{{{{{{{\rm{beam}}}}}}}}}g({\theta }_{1},{\phi }_{1})\frac{\sin {\theta }_{1}}{\cos {\theta }_{2}}d{\theta }_{1}d{\phi }_{1}.$$

(46)

where d = 1 AU is the distance from the Earth to the Sun. $\cos {\theta }_{2}=\sqrt{1-{R}_{S}^{2}\sin {{\theta }_{1}}^{2}/{d}^{2}}$ is the geometric relation. The R_S dependence in $\cos {\theta }_{2}$ is canceled out by the implicit R_S dependence in g(θ₁, ϕ₁, R_S). For the simplest scenario without scattering, $g({\theta }_{1},{\phi }_{1})=\delta ({\theta }_{1})/(2\pi \sin {\theta }_{1})$, Eq. (46) becomes ${S}_{{{{{{{{\rm{sig}}}}}}}}}={P}_{{{{{{{{\rm{sur}}}}}}}}}\cdot 1/{d}^{2}\cdot 1/{{{{{{{\mathcal{B}}}}}}}}\cdot d{{{{{{{\mathcal{P}}}}}}}}/d\Omega$, as expected. It is worth noting that since the data is averaged over the beams with flux larger than 50% of the maximum beam flux, the spherical surface integral is over the area covered by these selected beams, and then divided by the number of selected beams.

Statistics of robustness of background fitting parameter choosing

The upper limits on mono-chromatic signals, determined through the log-likelihood ratio test, exhibit robustness against variations in the parameters used for fitting the background. These parameters, denoted as n for the degree of the polynomial function and k for the number of bins included in the calculation, do not significantly affect the results. Typically, quadratic and trilinear polynomial forms are employed, yet even with these different choices, the outcomes remain largely unchanged. To demonstrate this resilience, we conducted a comprehensive analysis on LOFAR data collected on September 3, 2015, using varying degrees of polynomials and numbers of adjacent bins. Specifically, we examined three cases: 10 adjacent bins with a 3rd-degree polynomial, 10 adjacent bins with a 2nd-degree polynomial, and 8 adjacent bins with a 3rd-degree polynomial. The results of these analyses are presented in Fig. 8.

Fig. 8: The constraints on the monochromatic signal with different background fitting parameters.

The 95% C.L. upper limits from LOFAR data on September 3, 2015 with a constant mono-chromatic signal using different background fitting parameters. The orange, cyan and blue limits represent using 10 adjacent bins with a 3rd-degree polynomial, 10 adjacent bins with a 2nd-degree polynomial and 8 adjacent bins with a 3rd-degree polynomial, respectively, with n representing the degree of polynomial and k representing the number of adjacent bins.

Our investigation demonstrates that the derived signal limits exhibit a remarkably stability and are impervious to the specific choices of n and k. This robustness serves to reinforce the reliability and consistency of our method in establishing upper limits on the mono-chromatic signal from the LOFAR data.

Statistics of the Gaussian feature of LOFAR data

In our fitting process, the flux F(t_i, f_j) is characterized by its time index t_i and frequency index f_j. To analyze each frequency bin f_j, we calculate the average flux over time, denoted as $\bar{F}({f}_{j})$, and assume that it varies smoothly in frequency, which is fitted by using 3rd-degree polynomials. Within a fixed frequency bin, we consider the fluxes of different time bins to follow a Gaussian distribution, with $\bar{F}({f}_{j})$ serving as the mean of the Gaussian function.

To validate the assumption of Gaussian distribution, we specifically examine two frequency bins (j = 200, 400), corresponding to 49.21 MHz and 68.74 MHz on 3 September 2015, respectively. After undergoing the data cleaning process, each bin contains 920 and 1040 time bins, respectively. The top and bottom panels of Fig. 9 display the histograms for the flux at f₂₀₀ and f₄₀₀, respectively, and these plots align well with the Gaussian distribution.

Fig. 9: The distributions of flux in the LOFAR data.

They are obtained from time-bins after data cleaning process on September 3, 2015. a The distribution for the 200th bin with f₂₀₀ = 49.21 MHz is shown in the orange shaded region, while the Gaussian distribution with mean value $\bar{F}=1.589$, standard deviation σ = 0.0039 is shown in the solid blue line. b The distribution for the 400th bin with f₄₀₀ = 68.74 MHz is shown in the orange shaded region, while the Gaussian distribution with mean value $\bar{F}=2.994$, standard deviation σ = 0.0058 is shown in the solid blue line. These distributions exhibit a good fit to a Gaussian distribution.

Constraint on axion-like particle dark matter

In the axion DM model, the axion a, as a pseudo-scalar particle, interacts with the SM photon via

$${{{{{{{{\mathcal{L}}}}}}}}}_{a\gamma }=\frac{1}{2}{\partial }_{\mu }a{\partial }^{\mu }a-\frac{1}{2}{m}_{a}^{2}{a}^{2}+\frac{1}{4}{g}_{a\gamma \gamma }a{F}_{\mu \nu }{\tilde{F}}^{\mu \nu },$$

(47)

where ${\tilde{F}}^{\mu \nu }\equiv {\varepsilon }^{\mu \nu \alpha \beta }{F}_{\alpha \beta }/2$ is the dual EM field strength, m_a is the axion mass, a is the axion field and g_aγγ is the coupling strength between the axion and EM field. The last term in (47) can be simplified as − g_aγγaE ⋅ B.

Similar to the dark photon scenario, the probability of axion DM converting into photons is given by

$${P}_{a\to \gamma }({v}_{rc})=\pi \frac{{g}_{a\gamma \gamma }^{2}{\left\vert {{{{{{{{\bf{B}}}}}}}}}_{T}\right\vert }^{2}}{{m}_{a}}{v}_{rc}^{-1}{\left| \frac{\partial \ln {\omega }_{p}^{2}(r)}{\partial r}\right| }_{r={r}_{c}}^{-1}.$$

(48)

B_T is the magnetic field transverse to the direction of the axion propagation. The key difference from the dark photon case is that, the conversion of axions into photons requires the presence of a magnetic field. The probabilities (4) and (48) in the two cases are related via the expression

$$\sqrt{\frac{2}{3}}\epsilon {m}_{A{\prime} }^{2}\iff {g}_{a\gamma \gamma }\left\vert {{{{{{{{\bf{B}}}}}}}}}_{T}\right\vert {m}_{a}.$$

(49)

The Sun possesses a dipole-like magnetic field but suffers from large fluctuations^119,120. The global map of the magnetic field in solar corona obtained using the technique of the Coronal Multi-channel Polarimeter shows that the magnetic field strength is about 1-4 Gauss in the corona at the distance of 1.05-1.35 R_⊙¹²¹. In our case, the resonant conversion happens at the range of about 2.18-1.12R_⊙ (corresponding to frequencies in the range of 30–80 MHz; see Fig. 1). To proceed conservatively, we estimate $\left\vert {{{{{{{{\bf{B}}}}}}}}}_{T}\right\vert$ to be 1 Gauss at 1.05R_⊙ and extrapolate this value to obtain $\left\vert {{{{{{{{\bf{B}}}}}}}}}_{T}\right\vert \, \approx \, 0.11$-0.82 Gauss for our frequency range, following the attenuation relation ∝ R⁻³.

The upper limit for the dark photon case can be directly translated into that for the axion case using the relation (49). We adopt $\left\vert {{{{{{{{\bf{B}}}}}}}}}_{T}\right\vert$ as a function of distance using the extrapolation above. Consequently, we plot the constraint on g_aγγ in Fig. 10. However, there is a large uncertainty in our estimation of the magnetic field, which overshadows other statistical and systematic uncertainties. Therefore, the zig-zag features shown in Fig. 4 become less meaningful in the axion DM case. As a result, in Fig. 10, we average the upper limits over every 20 frequency bins to indicate the sensitivity of the LOFAR data on axion DM model. The resulting graph shows that while our limit exceeds the existing constraints from Light-Shining-through-a-Wall experiments, including CROWS¹²², ALPS¹²³, and OSQAR¹²⁴, it is not as competitive as the direct search experiments such as CAST¹²⁵ or ADMX SLIC¹²⁶ (in very narrow bands), and the astrophysical bounds from observations of magnetic white dwarf polarization⁷⁰, Globular Clusters^82,83, pulsars⁶², as well as quasars and blazars^75,76,77.

Fig. 10: The constraints on the parameter space of axion-like particle dark matter.

95% C.L. upper limit on axion-photon coupling g_aγγ from 17 minutes observation of LOFAR data is shown in the cyan shaded region. We also show the existing constraints (summarized in ref. ⁵²) from various experiments and astrophysical observations in gray shaded regions, including Light-Shining-through-a-Wall experiments: CROWS¹²² (95%), ALPS¹²³ (95%), and OSQAR (95%)¹²⁴; helioscope: CAST (95%)¹²⁵; haloscope: ADMX SLIC (90%)¹²⁶; astrophysical bounds: magnetic white dwarf polarization (MWDP) (95%)⁷⁰, Globular Clusters (95%)^82,83, pulsars (95%)⁶², as well as quasars and blazars (QBs, shown in dashed gray) (95%)^{75 -- 77}. Different constraints may choose different confidence levels, and we keep their original choice unchanged as labeled in the parentheses following each experiment. The ADMX SLIC constraint assumes axions to be dark matter, ρ = 0.45 GeV cm⁻³, and we have rescaled it to be 0.3 GeV cm⁻³ in the plot for comparison.