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

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−13 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.

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].
In this work, we investigate such resonantly converted monochromatic radio signal within the solar observation data collected by Low Frequency Array (LOFAR) telescope [87].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 −13 within the frequency range of 30 − 80 MHz.

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 where A ′ and γ represent dark photon and photon respectively, F µν and F ′µν represent the field strengths of photon and dark photon respectively, with the Greek letters µ, ν denoting the vector indices, m A ′ denotes the mass of dark photon, A ′ µ 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 nonrelativistic plasma, and can be represented as 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 ′ propagates in the plasma, it can resonantly convert into a SM photon when ω p ≈ m A ′ [86].In the solar corona, we have n e monotonically decreasing from 10 10 to 10 6 cm −3 with increasing height above the solar photosphere.Therefore, the corresponding plasma frequency scans from 4 × 10 range, the resonant conversion of DPDM into EM waves can occur at a specific radius r c satisfying ω p (r c ) = m A ′ .The frequency of the converted EM wave, m A ′ /(2π), 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 [87] and SKA [88].Since DM in the Galactic halo is nonrelativistic with the typical velocity v DM approximately 10 −3 times the speed of light, the converted EM wave is nearly monochromatic with a spread of about 10 −6 around its central value [86].The DM dispersion bandwidth B sig can be evaluated by Our analysis adopts the electron density profile for the quiet Sun provided by LOFAR observations [89], 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,91] where v rc is the radial velocity at the resonant layer.The prefactor 2/3 arises in (4) because the longitudinal mode of photons converted from the corresponding mode of dark photon cannot propagate out of the plasma.The conversion probability (4) accounts for two transverse modes, and we assume that dark photons polarize in three directions with equal probability.
Utilizing the conversion probability, the radiation power P per solid angle dΩ at the conversion layer can be derived as where the DM density is ρ DM = 0.3 GeV cm −3 [92,93] and v 0 , the initial DM velocity, follows a Maxwellian distribution f DM (v 0 ) with the most probable velocity of 235 km/s [94,95].v(r c ) = v 2 0 + 2G N M ⊙ /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 ′ → γ 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.[86].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 [96,97].We are using LOFAR data made in the tied-array beam mode.While this mode offers a nice angular resolution [87], 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 where d = 1 AU is the distance between Earth and the Sun.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 res = 97 kHz.In our case, B sig ≪ B res , so we have B = B 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 [97].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 [98] with the speed C s ≈ [T e (1 + 3T 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 −3 c in (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 raytracing simulation of radio photon propagation [97] considers angular dispersion due to inhomogeneities while ignoring energy dispersion.
In the context of ( 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 (6) parameterizes the scattering effect and is defined as where g(θ 1 , ϕ 1 ) 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 (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 (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 raytracing method developed in Ref. [97].We describe the scattering process of radio waves using the Fokker-Planck and Langevin equations based on the Hamilton equations for photons [97][98][99].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 Reference [96].Additionally, we consider the anisotropic density fluctuation magnitude as α anis = 0.1 [96].Here, α anis represents the anisotropy parameter, which is the ratio between the perpendicular and parallel correlation lengths [97].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.

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 [87], 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 ′ at 32 MHz [87].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 April 25, 2015, July 3, 2015, and September 3, 2015.The bandwidth is B 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 in-terval 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 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 Ōi and the standard deviation σ Ō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 σ sys i .The total uncertainty is in the quadrature form, We adopt the log-likelihood ratio test method [100] to set upper limits on the DPDM parameter space.We construct the likelihood function for a specific frequency bin i 0 in the Gaussian form [101] L(S, a) where B(a, f i ) is the polynomial function used for background fitting, the coefficients a = (a 1 , a 2 , a 3 ) are treated as nuisance parameters, S denotes the assumed DPDMinduced signal at bin i 0 , and δ ii0 is the Kronecker delta.
We then build the following test statistic [100,101] In the denominator, the likelihood L gets maximized at a = â and S = Ŝ; in the numerator, L gets maximized at a = ã for a specified S. The test statistic q S follows the half-χ 2 distribution, with the probability density function the cumulative distribution function of which is given by H(q S |S) = 1/2 1 + erf q S /2 , where erf(x) is the Gauss error function.Then, we can define the following criterion [100,101]: 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.6MHz 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.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 (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 −13 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 [102].
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.We also show the existing constraints (summarized in Ref. [52]) from the CMB distortion (95%) [6,53], the haloscope searches WISPDMX (95%) [103], and Dark E-field experiment (5σ) [102] 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 −3 , the same as our choice, and are scaled by the dark photon density (ρ A ′ /ρDM) 0.5 .

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 −13 , 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 LO-FAR 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.Ter-restrial 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 [104] satellite and the Parker Solar Probe [105], equipped with radio spectrometers, could offer an avenue for DPDM detection.

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 LO-FAR observations [89], which employed ray-tracing simulations to fit the solar intensity profile observed by LO-FAR 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. [90], are derived from the temperature (T ) and hydrogen density (n H ) profiles for the quiet Sun, based on the photospheric model from Ref. [106] and the coronal model from Refs.[107,108].The consistency and validation of these profiles have been confirmed by various research groups [109] using a chromosphere model from Ref. [110] and the coronal model from Refs.[107,108].These independent calculations consistently agree with each other and have been validated by observations of atomic lines in the soft X-ray range [111] and extreme ultraviolet range [106].
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 [106,111].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 [89] with the parameters N 0 and H 0 , and the latter is defined as where R ⊙ is the solar radius, k B is the Boltzmann constant, g ⊙ = 274 m s −2 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 [112].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 0 and the temperature T .To carry out our calculation, N 0 = 1.6 × 10 11 m −3 and T = 2 × 10 6 K are used from LOFAR observation fit [89].
In Fig. 1, we present a comparison between the profile we adopted from LOFAR observations [89], the profile from Ref. [90], the hydrostatic profile modeled by Eq. ( 12), and the r −2 profile.The r −2 density profile indicates that a constant solar wind speed has been attained [89].It can be seen that the solar profile from LO-FAR observations closely matches the hydrostatic profile at the high-frequency range, but exhibits a slower decline and transitions into the r −2 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 [89].
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 ′ → γ 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,91,113].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: µ , where A µ represents the vector field of photon.In this basis, the equations of motion become 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 T under a plane wave solution A(r, t) = Ã(r)exp(−iωt + ik r r) with frequency ω and wavenumber k = ω 2 − m 2 A ′ , 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 where This equation can be solved perturbatively by expanding the time-evolution operator in Dyson series [91].At the first-order, the conversion probability is given by [113] This formula can be further simplified to Eq. ( 4) by using saddle point approximation where f ′ (r 0 ) = 0.The thickness of resonant layer is on the order of 10 3 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 where z is the impact parameter for DPDM, b = r c v(r c )/v 0 is the maximum impact parameter for A ′ to reach the conversion layer at r = r c , v(r c ) = v 2 0 + 2G N M ⊙ /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 . 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 ′ → γ conversion probability by altering ⃗ ∇ω p ; (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 ′ →γ .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 ′ →γ .On the other hand, more resonant points where ω p (r) = m A ′ are introduced, which increases P A ′ →γ .It turns out that these two influences cancel each other out, resulting in P A ′ →γ 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. [96], 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 [96] ϵ e ≡ ∆n e /n e ≈ 10%, (20) and importantly, this fluctuation fraction remains constant as the radial distance changes [96,114].
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 [115,116] where C 2 N is the structural constant, q represents the spatial wavenumber, and α corresponds to the powerlaw 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 6 l i [96].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 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 ⟨∆n 2 e ⟩ ≡ (ϵ e n e ) 2 as, Using the Kolmogorov spectrum with α = 11/3, it can be determined that The fluctuations can be expressed in the Fourier modes, ∆n e (r) = qi 1 dq ∆ñ e (q)e iqr . ( We have rescaled the momentum q ≡ q/q o for convenience.Subsequently, the fluctuations averaged over the length scale Compared with Eq. ( 23), we have the fluctuations in the momentum space, The averaged derivative of n e (r) with respect to r, in the squared form, is then dq ′ ∆ñ e (q)∆ñ e (q ′ )qq ′ e i(q−q ′ )r 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 dq ′ ∆ñ e (q)∆ñ e (q ′ )q 2 q ′2 e i(q−q ′ )r 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 which turns out to be much larger than the dark photon wavelength, δl e k A ′ ≈ 30 ≫ 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, δl res = 2π/F ′′ (r c ).It is defined as the length along which the phase factor in Eq. ( 17), , changes by π.This is the length interval which dominantly contributes to P γ→A ′ .We have which obviously satisfies δl res ≪ δl e .Next, we evaluate the robustness of the saddle-point approximation.The crucial criterion is that the second derivative F ′′ (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 ′ (r) = 0).Then we calculate the following quantity with the help of Eqs. ( 29) and ( 30), 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. [117] 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 γ→A ′ remains unchanged.As stated before, there are two new effects that counteract each other in modifying the value of the probability P A ′ →γ : 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 dn to denote the ratio between P A ′ →γ with and without density fluctuations, and it can be calculated as We can provide a rough estimate of r 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 dn ≈ (∆r/δr) • (δr/δn e )/(∆r/δ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.
We then proceed to numerically compute the ratio r 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 dn ≃ 1 is indeed verified.For the convenience of numerical computation, we first need to discretize the density fluctuations, Eq. ( 23), as where ⟨∆n 2 e ⟩ (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 )+∆].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 )/∆.Next, the variance of density fluctuations for different momentum modes can be estimated as The density as a function of distance, n e (r), with the fluctuations taken into consideration, can be modeled as n e (r) = (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 ), 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 σ ne (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 −2 as shown in Eq. (38).
In this example, the corresponding resonant conversion layer is at r c ≃ 9.4 × 10 5 km from the solar center, which corresponds to q 0 r c ≃ 40.
To proceed, we take ∆ = 0.2 and thus we have N = 30 modes in total.However, computing r 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 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 dn with more modes gradually included in computations.We that the average value of r 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 ′ →γ , as the two effects of larger derivatives and more resonant points cancel out each other.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 LO-FAR stations Firstly, the Field of View (FOV) of LOFAR, or effectively, the Full Width Half Maximum (FWHM) of LOFAR, is determined by where λ is the observation wavelength, the coefficient η = 1.02 [118], and D ≃ 3.5 km is the station diameter according to Ref. [119].Therefore, the FWHM (for one beam) is approximately 10 −3 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 dP/dΩ × 4π.Therefore, the survived power at the last scattering sphere is given by Considering a virtual source point P 1 situated within a surface element dA 1 on the last scattering sphere (as depicted in the schematic diagram of Fig. 7), the power it radiates in the direction r is where the angular distribution function g(θ 1 , ϕ 1 ) accounts for the fact that after multiple random scattering events, the radiation from the surface element is not . Schematic diagram of the propagation of photons after the last scattering.RC denotes the conversion layer, and RS denotes the last scattering sphere.A surface element dA1, which containing a point P1, acts as the radiation source on the last scattering sphere.θ is the polar angle of P1.Another surface element dA2, encompassing P2, serves as the detection area on the Earth, which defines a solid angle dΩ1 about P1 in the direction of r. θi is the angle between the propagation vector r and the vector N i of dAi.The direction of N2 is aligned with the line connecting the centers of the Sun and the Earth.d is the distance from the Earth to the Sun.
simply in the radial direction.g(θ 1 , ϕ 1 ) is normalized as The relation dΩ 1 = dA 2 cos θ 2 /r 2 is useful where the cosine factor accounts for converting the receiving area dA 2 to the projected area normal to r.Then, Eq. ( 42) becomes where r is the distance from the surface element to the Earth.Meanwhile, since θ 2 is on the order of 10 −3 rad, it follows that cos θ 2 ≃ 1.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: As discussed in the main text, the angular distribution function g(θ 1 , ϕ 1 ) 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 where d = 1 AU is the distance from the Earth to the Sun.cos θ 2 = 1 − R 2 S sin θ 1 2 /d 2 is the geometric relation.The R S dependence in cos θ 2 is canceled out by the implicit R S dependence in g(θ 1 , ϕ 1 , R S ).For the simplest scenario without scattering, g(θ 1 , ϕ 1 ) = δ(θ 1 )/(2π sin θ 1 ), Eq. ( 46) becomes S sig = P sur • 1/d 2 • 1/B • dP/dΩ, 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.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 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 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 September 3, 2015, respectively.After undergoing the data cleaning process, each bin contains 920 and 1040 time bins, respectively.The top and bottom panels of 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.
Fig. 9 display the histograms for the flux at f 200 and f 400 , respectively, and these plots align well with the 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 where F µν ≡ ε µναβ F αβ /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 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 The Sun possesses a dipole-like magnetic field but suffers from large fluctuations [120,121].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.35R ⊙ [122].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 |B T | to be 1 Gauss at 1.05R ⊙ and extrapolate this value to obtain |B T | ≈ 0.11-0.82Gauss for our frequency range, following the attenuation relation ∝ R −3 .
The upper limit for the dark photon case can be directly translated into that for the axion case using the relation (49).We adopt |B T | 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 .The constraints on the parameter space of axion-like particle dark matter.95% C.L. upper limit on axion-photon coupling gaγγ from 17 minutes observation of LOFAR data is shown in the cyan shaded region.We also show the existing constraints (summarized in Ref. [52]) from various experiments and astrophysical observations in gray shaded regions, including Light-Shining-through-a-Wall experiments: CROWS [123] (95%), ALPS [124] (95%), and OS-QAR (95%) [125]; helioscope: CAST (95%) [126]; haloscope: ADMX SLIC (90%) [127]; astrophysical bounds: magnetic white dwarf polarization (MWDP) (95%) [70], Globular Clusters (95%) [82,83], pulsars (95%) [62], as well as quasars and blazars (QBs, shown in dashed gray) (95%) [75][76][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 −3 , and we have rescaled it to be 0.3 GeV cm −3 in the plot for comparison.
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 [123], ALPS [124] and OSQAR [125], it is not as competitive as the direct search experiments such as CAST [126] or ADMX SLIC [127] (in very narrow bands), and the astrophysical bounds from observations of magnetic white dwarf polarization [70], Globular Clusters [82,83], pulsars [62], as well as quasars and blazars [75][76][77].

−6 to 4 ×Figure 1 .
Figure1.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[89].In comparison, the solid orange line represents the profile from Ref.[90].The dashed cyan line represents a simple hydrostatic model and the dashed red line represents an r −2 profile[89].The gray shaded region denotes the frequency region around 30−80 MHz which our study focuses on.

Figure 2 .
Figure 2. The propagation coefficients as functions of frequency.The survival probability Psur and the smearing factor β as functions of the photon frequency f are depicted as the solid black line and the solid purple line, respectively.

2 Figure 3 .
Figure 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 April 25, 2015, July 3, 2015, and September 3, 2015 are represented by the blue, orange, and green curves, respectively.

Figure 4 .
Figure 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 ′ from 17 minutes observation of LOFAR data is shown in the cyan shaded region.We also show the existing constraints (summarized in Ref.[52]) from the CMB distortion (95%)[6,53], the haloscope searches WISPDMX (95%)[103], and Dark E-field experiment (5σ)[102] 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 −3 , the same as our choice, and are scaled by the dark photon density (ρ A ′ /ρDM) 0.5 .

Figure 5 .
Figure 5.The electron density profile exhibiting fluctuations.The electron density ne 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 ne(rc) 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.

Figure 6 .
Figure 6.The ratio between the conversion probabilities with and without density fluctuations.The ratio r 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.

Figure 8 .
Figure 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.

FluxFluxFigure 9 .
Figure 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 f200 = 49.21MHz is shown in the orange shaded region, while the Gaussian distribution with mean value F = 1.589, standard deviation σ = 0.0039 is shown in the solid blue line.b The distribution for the 400th bin with f400 = 68.74MHz is shown in the orange shaded region, while the Gaussian distribution with mean value F = 2.994, standard deviation σ = 0.0058 is shown in the solid blue line.These distributions exhibit a good fit to a Gaussian distribution.