Improved axion haloscope search analysis

The most significant and practical figure of merit in axion haloscope searches is the scanning rate, because of the unknown axion mass. Even under the best experimental parameters, the only way to improve the figure of merit is to increase the experimentally designed signal to noise ratio in the axion haloscope search analysis procedure. In this paper, we report an improved axion haloscope search analysis, developed using CAPP-8TB. By fully incorporating the known negative correlations induced by the background parametrizations, the figure of merit of the CAPP-8TB haloscope search was increased by about 17% for the scanning rate, with about an 8% improvement in the signal to noise ratio. In addition, the physics results can be retrieved directly from the standard Gaussian statistics of the data, at any confidence level.


I. INTRODUCTION
The axion [1] is an elementary particle considered to result from a breakdown in a new symmetry first proposed by Peccei and Quinn (PQ symmetry) [2] to solve the strong CP problem in the standard model of particle physics [3]. This particle is massive, stable, and was born cold by the PQ symmetry breaking, which makes the axion one of the most promising candidates for cold dark matter (CDM), which constitutes about 85% of the matter in the Universe according to cosmological measurements and the standard model of big bang cosmology [4].
The method of searching for the axion proposed by Sikivie [5], also known as the axion haloscope search, involves a microwave resonant cavity with a strong static magnetic field that induces axions to convert to microwave photons. Using the resonant cavity, the axion signal power can be enhanced when the axion mass m a matches the resonant frequency of the cavity mode ν, m a = hν/c 2 . Because the axion mass is unknown, however, the resonant cavity has to be tunable, to allow axion haloscope searches to scan all frequencies corresponding to possible axion masses. Because of this frequency scanning procedure, the most significant and practical figure of merit in axion haloscope searches is the scanning rate [6] for a target signal to noise ratio SNR target . P aγγ a is the axion signal power proportional to B 2 V CQ L for an axion-photon coupling strength [5,6]. The details of B, V , C, and Q L can be also found there [5,6], P n is the noise power proportional to the noise temperature T n , and b a is the axion signal window. * Now at Korea Astronomy and Space Science Institute, Daejeon 34055, Republic of Korea † Corresponding author : brko@ibs.re.kr From the radiometer equation [7], the SNR is where σ Pn is the noise power fluctuation. The subscript "designed" stands for experimentally designed, thus the SNR designed must be designed to be the same as the SNR target or the SNR target must be set to be the SNR designed in order to have an axion haloscope search in an experimentally designed time. In general, the experimentally achieved signal to noise ratio SNR achieved , however, is smaller than the designed one by holding the relation SNR achieved = ǫ SNR SNR designed , where ǫ SNR is the reconstruction efficiency of the SNR designed in an axion haloscope search analysis procedure whose values vary from about 50% to 80% depending on the analysis strategy [8][9][10].
Having said that, the scanning rate guides the following two cases; (I) in order for axion haloscope searches to achieve the target sensitivity or SNR target in an experimentally designed time, the experimental parameters, B, V , C, Q L , and T n have to be designed to meet the condition SNR achieved = SNR target , (II) in order for axion haloscope searches to achieve the target sensitivity or SNR target , axion haloscopes have to take data until the condition SNR achieved = SNR target is satisfied, under the experimental parameters, B, V , C, Q L , and T n at their best.
In both cases, the figure of merit of the experiments can be enhanced by improving the SNR achieved or ǫ SNR , which results in more sensitive results for (I) and shorter data acquisition periods for (II).
In this paper, we report an improved axion haloscope search analysis adopted in the recent axion dark matter search with the CAPP-8TB haloscope [11]. With the improved axion haloscope search analysis, the figure of merit of the CAPP-8TB haloscope search was effectively increased by about 17%.

II. AXION HALOSCOPE SEARCH ANALYSIS STRATEGY
The simplest axion haloscope search analysis strategy is the one-bin search that was employed in Ref. [10], where all the axion signal power belongs to a single frequency bin width corresponding to the axion signal window b a , if axions are there. The price for the simplicity in the one-bin search is a very low reconstruction efficiency of the axion signal power ǫ P aγγ a . As described in Refs. [8,10], we lose about 20% of the signal power by choosing a signal window to get an optimized SNR, and an additional 20% from the frequency binning choice.
The other strategy is the multi-bin co-adding search developed by the Axion Dark Matter eXperiment (ADMX) [8], where all the axion signal power is distributed over the multi frequency bins obeying an axion signal shape [12] (also shown in Fig. 1). The multi frequency bin width also corresponds to the axion signal window b a . The multi-bin co-adding search overcomes the inefficiency caused by the selection of frequency binning. We can retain about 99.9% of the total signal power by taking the signal window of 5000 Hz or, equivalently, 10 frequency bins with a resolution bandwidth of 500 Hz.
The multi-bin co-adding search with signal likelihood weighting developed by the Haloscope at Yale Sensitive to Axion CDM (HAYSTAC) [9] overcomes the inefficiency resulting from the choice of signal window.
As reported in [8][9][10], however, the background subtraction with a five-parameter fit [8,10] or a Savitzky-Golay filter [9] also generates an additional inefficiency of about 20% in the SNR reconstruction in the axion haloscope search analysis procedure. Improving the inefficiency induced by the background subtraction might be another low-cost innovation in axion haloscope searches, and that is the main contribution of this paper.

III. AXION HALOSCOPE SEARCH ANALYSIS PROCEDURE
The axion haloscope search analysis procedure can be divided into the following three steps; Step-1 ; background parametrization for the background subtraction 1 , Step-2 ; combining all the power spectra as a single power spectrum, taking into account the overlaps among the power spectra, Step-3 ; constructing a "grand power spectrum" by coadding multi-bins with the signal likelihood weighting.
After the background subtraction in Step-1, the normalized power excess at each stage is likely to follow a Gaussian distribution in the absence of axion signals. The mean of the Gaussian depends on the expected SNR, while its width should be unity with the correct error estimation. All of the effort in this work was made with the goal of obtaining such a Gaussian distribution in each step, from Step-1 to Step-3.

IV. DATA AND PARAMETERS
Here we describe the data used for this work; (i) experimental data from the CAPP-8TB experiment [11], (ii) 5000 simulated axion haloscope search experiments with a flat background and axion signals at a particular frequency, referred to as "flat background simulation data", (iii) 5000 simulated axion haloscope search experiments with the CAPP-8TB backgrounds and axion signals at a particular frequency, referred to as "CAPP-8TB simulation data", (iv) 5000 simulated axion haloscope search experiments with the CAPP-8TB backgrounds only, referred to as "background only simulation data".
Here we also include some assumptions and parameters adopted by the CAPP-8TB axion haloscope search [11]. We assume the axions have an isothermal distribution, thus distribute over a boosted Maxwellian shape as shown in Fig. 1 with the parameters, an axion rms velocity of about 270 km/s and the Earth rms velocity of 230 km/s with respect to the galaxy frame [12]. We took 5000 Hz as the signal window of the axions considered in CAPP-8TB, then can retain about 99.9% of the total signal power as shown in Fig. 1. In the Step-2 and Step-3 procedures, the resolution bandwidth (RBW) was set to be 500 Hz [11]. With our 10 co-adding bins in Step-3, therefore, P aγγ a is retained almost 100%.

V. VALIDATIONS
Using the flat background simulation data, first, we validated our understanding of the axion haloscope search analysis procedure, especially the signal weighting in the Step-3 procedure as well as our simulation data, particularly the signal injections that must reflect our cavity responses. With a flat background, one can easily predict solid SNRs with little background dependence in the Step-2 procedure. Figure 2 shows the predicted SNR as a function of the number of co-adding bins with and without the signal weighting. The weighting factors were obtained by integrating over the relevant regions of the axion signal likelihood shown in Fig. 1. The flat background simulation data was fed through the procedure described in Sec. III. The backgrounds were subtracted without any fit parametrizations. All the background-subtracted power spectra were then combined as a single power spectrum. In constructing a grand power spectrum, each of the backgroundsubtracted and combined power spectral line was weighted accordingly [9] using the axion signal likelihood shown in Fig. 1. It was further normalized by the corresponding noise fluctuation, which was also weighted according to the axion signal shape [9]. With a signal window of 5000 Hz or, equivalently, 10 co-adding bins, Fig. 3 shows two distributions of the normalized power excess from a particular frequency where we put simulated axion signals, with (solid circles) and without (rectangles) the signal weighting, respectively. Figure 3 also shows a Gaussian fit results, including the means (µ) and widths (σ) of the distributions. Both means follow the predicted SNR shown in Fig. 2 and the widths are unity as they should be. Having demonstrated our solid understanding of the simulation data and axion haloscope search analysis procedure from the flat background simulation data, we applied the same analysis procedure to the CAPP-8TB simulation data as well as the CAPP-8TB experimental data, but with a five-parameter fit [8] for the background parametrization and subtraction. The same procedure except for the background subtraction using the simulation input background functions, i.e., perfect fit was also applied just to the simulation data to determine the inefficiency induced from the background subtraction [8][9][10]. The normalized power excess distributions from the CAPP-8TB experimental data are shown as triangles in the 1st column of Fig. 4. The distributions from a single simulated axion haloscope search experiment with CAPP-8TB backgrounds and axion signals at a particular frequency are in the 2nd and 3rd columns of Fig. 4 as inverted triangles. The distributions in the 1st and 3rd columns were obtained after subtracting the background with a five-parameter fit, while in the 2nd they were obtained with a perfect fit. From top to bottom, they are the distributions after the Step-1, Step-2, and Step-3 procedures, respectively. The distributions after the Step-3 procedure are narrower than unity for both the experimental and simulation data with a fiveparameter background parametrization. This has been observed in axion haloscope search experiments [8][9][10]. Without shrinking the distribution after the Step-3 procedure with perfect fit (bottom center of Fig. 4), the narrower width of the normalized power excess distribution must be induced from the background parametrization and subtraction thereafter [8][9][10].
Our validations shown in Figs. 3 and 4 demonstrate not only our solid understanding of the analysis procedure but also the validity of our simulation data.

VI. IMPROVEMENT
Equation (3) shows the power excess in one of the frequency bins of the grand power spectrum with the signal weighting [9] P grand = co-adding bins and Eq. (4) shows the full error propagation of the power excess shown in Eq. (3) where P i and σ Pi are power excess and associated error in the i th frequency bin of the combined power spectrum after the Step-2 procedure, respectively, and L i is an axion signal likelihood for both P i and σ Pi weightings. With the signal window b a of 5000 Hz as shown in L i ≃ 1 with 10 co-adding frequency bins and an RBW of 500 Hz for the CAPP-8TB axion haloscope search. ρ ij is the correlation coefficient between the frequency bins participating in the co-adding procedure, thus is unity for i = j. The narrower width of the normalized power excess distribution from the normalized grand power spectrum means that the noise fluctuations in the grand power spectrum σ P grand are overestimated with the terms for i = j only in Eq. (4) or, equivalently, with Pi L 2 i only. The correlation coefficients ρ ij in Eq. (4) are very likely to be negative due to overestimation of the noise fluctuations in the grand power spectrum [8][9][10].
With 10 co-adding bins, each frequency bin in the grand power spectrum has a full error contribution from 10 × 10 combinations or, equivalently, from a relevant 10 × 10 correlation matrix. The total number of frequencies in our grand power spectrum is 100001 with a frequency band of 50 MHz (1600 to 1650 MHz) and an RBW of 500 Hz [11], hence we need to construct 100001 10 × 10 correlation matrices. These 100001 correlation matrices were constructed from the CAPP-8TB simulation data and background only simulation data. The elements in In the top 2-D plots, ρij for i = j are unity with 100% correlation, but are set to zero to distinguish the differences among ρij for i = j. ρij for i = j are also left out in the bottom histograms and that is the reason the entries in them are not 100, but 90 from a 10×10 matrix. Plots in the 1st column were obtained from a frequency that has axion signals, and those in the next three columns from other frequencies with no axion signals.
each correlation matrix were calculated as the standard Pearson correlation coefficient. The CAPP-8TB simulation data was used for the background parametrizations and the parametrizations were then applied to the background only simulation data to extract the power excess for the correlation coefficient calculations. Figure 5 shows examples of the correlations obtained from the simulation data. The 1st column shows the correlation coefficient map and the distribution of the correlation coefficients from a frequency that has input axion signals, and the other columns show them from frequencies that have no axion signals. As predicted earlier, most of the correlation coefficients were constructed as negative values, which induces the σ P grand overestimation or, equivalently, the narrower widths of the normalized power excess distributions shown in the bottom left (experimental data) and right (simulation data) of Fig. 4. Figure 5 also shows that the negative correlations are stronger in the signal region, which also makes sense given the following. In the presence of a signal on top of the background, the background level in the signal region, in general, can be overestimated by any background only fit parametrizations. The background subtraction with the overestimated background level then can result in a more highly underestimated power excess overall, i.e., more negatively biased power excess, in the signal region. This will be confirmed with later simulations (see Fig. 8). The 100001 correlation matrices constructed from the simulation data were then fed through the CAPP-8TB experimental data and CAPP-8TB simulation data in the Step-3 procedure. Figure 6 was obtained from the CAPP-8TB experimental data, which corresponds to the bottom left data in Fig. 4. After the full incorporation of the correlation matrices the triangles now follow a standard Gaussian distribution. With the standard Gaussian distribution from the data, one can retrieve the physics results at any confidence level directly from the statistics  Figure 7(c) was also from 5000 simulated CAPP-8TB experiments. Both the inverted triangles and solid circles follow a Gaussian distribution with a width exhibiting unity, thanks to the full incorporation of the correlation matrices. From the mean values of the solid circle distributions in Figs. 7(a), 7(b), and 7(c), the SNR efficiencies with a five-parameter fit are 90.7% and 84.0% with and without the correlation matrices, respectively. The former shows a relative improvement of about 8% in ǫ SNR with respect to the latter.
We also obtained SNR as a function of frequency from the CAPP-8TB simulation data, as shown in Fig. 8 that shows around the signal region only. The blue circles and red rectangles in Fig. 8 were obtained with perfect fit and a five-parameter fit, respectively. The red rectangles show more negatively biased power excess around the signal region as pointed out earlier.
According to the radiometry equation shown in Eq. (2), the inefficiency in ǫ SNR can be from the reconstruction efficiency of the axion signal power ǫ P aγγ a or an overestimation of the noise fluctuation σ Pn , or both. With the σ Pn improvement from this paper, however, the reduced SNR shown as red rectangles in Fig. 8 must be from the reconstruction efficiency of the axion signal power ǫ P aγγ a that was deducted from a five-parameter fit parametrization and subtraction thereafter.

VII. SUMMARY
We report an improvement in axion haloscope search analysis developed using CAPP-8TB. The improvement results from the full incorporation of the correlation matrices between frequency bins participating in the coadding procedure.
This improvement not only effectively increased the figure of merit of the CAPP-8TB haloscope search by about 17%, with an ǫ SNR improvement of about 8%, but also allows physics results at any confidence level to be retrieved directly from the standard Gaussian statistics of the data, without having to rely on Monte Carlo simulations.

VIII. ACKNOWLEDGMENTS
This work is supported by the Institute for Basic Science (IBS) under Project Code No. IBS-R017-D1-2020-a00.