Improved axion haloscope search analysis

One of the most significant and practical figures of merit in axion haloscope searches is the scanning rate, because of the unknown axion mass. 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 using the data taken by the CAPP-8TB haloscope. By correcting for the background biased by the background parametrizations in the presence of axion signals, we realized a signal to noise ratio efficiency of about 100%. Given the axion haloscope search analyses to date, the scanning rate can be improved by 21%, with about a 10% improvement in the signal to noise ratio. This improvement is another low cost innovation in axion haloscope searches, where all the experimental parameters are currently at their best.


Introduction
The axion [1,2] is an elementary particle considered to result from a breakdown in a new symmetry first proposed by Peccei and Quinn (PQ symmetry) [3] to solve the strong CP problem in the Standard Model of particle physics [4][5][6][7][8]. This particle is massive, stable, and is born cold by the PQ symmetry breaking. This makes the axion one of the most promising candidates for cold dark matter, which constitutes about 85% of the matter in the Universe according to cosmological measurements and the standard model of Big Bang cosmology [9].
The method of searching for the axion proposed by Sikivie [10,11], also known as the axion haloscope search, involves a microwave resonant cavity with a strong static magnetic field that induces axions to convert into 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 [12] for a target signal to noise ratio SNR target . P aγγ a is the axion signal power proportional to g 2 aγγ B 2 V CQ L [10][11][12], where g aγγ is the axion-photon coupling strength, B is a static magnetic field provided by magnets in the axion haloscopes, V is the cavity volume, C is a form factor representing the overlap between the electric field of the cavity mode and JHEP04(2021)297 the static magnetic field whose general definition can be found in ref. [13], and Q L is the loaded quality factor of the cavity mode. P n is the noise power proportional to the noise temperature T n and the axion signal window b a .
From the radiometer equation [14], 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. However, generally the experimentally achieved signal to noise ratio SNR achieved 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 90% depending on the analysis strategy [15][16][17][18][19].
Having said that, the scanning rate guides the following two cases.
(I) 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) 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, even 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 using the data taken by the CAPP-8TB haloscope [20]. Using the improved axion haloscope search analysis in this paper, the figure of merit of axion haloscope searches can be effectively increased by about 21%.

Axion haloscope search analysis strategy
The simplest axion haloscope search analysis strategy is the one-bin search that was employed in ref. [19], 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 of this one-bin search is a very low reconstruction efficiency of the axion signal power P aγγ a .
As described in refs. [15,16,19], 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) [15,16], where all the axion signal power is distributed over the multi frequency bins obeying an axion signal shape [21] (also shown in figure 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. The multi-bin co-adding search with signal likelihood weighting developed by the Haloscope At Yale Sensitive To Axion Cold dark matter (HAYSTAC) [17,18] overcomes the inefficiency resulting from the choice of signal window.
As reported in [15][16][17][18][19], however, the background subtraction with a 5-parameter fit [15,16,19] or a Savitzky-Golay filter [17,18] also generates an additional inefficiency of about 20% in the SNR reconstruction in the axion haloscope search analysis procedure. The HAYSTAC analysis procedure recovered this inefficiency using a single frequencyindependent numerical factor, resulting in an improvement of about 8% in SNR [17,18]. Another low cost innovation in axion haloscope searches might be to improve the inefficiency induced by the background subtraction, and that is the main contribution of this paper.

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 co-adding multi-bins with 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, because the background in axion haloscope searches is very like stationary Poisson noise. The mean of

JHEP04(2021)297
the Gaussian depends on the expected SNR, while its width should be unity if the power excess errors are estimated correctly. All of the effort in this work is focused on obtaining such a Gaussian distribution in each step, from Step-1 to Step-3.

Data and parameters
Here we describe the data used for this work. It consists of: (i) experimental data from the CAPP-8TB experiment [20], (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", (v) other 4 sets of 5000 simulated axion haloscope search experiments with the CAPP-8TB backgrounds, where each set has axion signals at different frequencies, thus has different SNRs, referred to as "other CAPP-8TB simulation data".
Here we also employ some assumptions and parameters adopted by the CAPP-8TB axion haloscope search [20]. We assume the axions have an isothermal distribution, thus distribute over a boosted Maxwellian shape, as shown in figure 1, with the following 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 [21]. We took 5000 Hz as the signal window of the axions considered in CAPP-8TB. We can then retain about 99.9% of the total signal power as shown in figure 1. After Step-1, the five nonoverlapping frequency bins in each backgroundsubtracted power spectrum were merged so that the resolution bandwidth (RBW) became 500 Hz, which we refer to as "Step-1.5". An RBW of 500 Hz was chosen not only to retain the axion signal shape shown in figure 1, but also to avoid unnecessarily inordinate calculation time. The Step-2 and Step-3 procedures then went through with the RBW of 500 Hz [20]. With our 10 co-adding bins in Step-3, therefore, P aγγ a is almost 100% retained.

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 figure 1.

JHEP04(2021)297
number of co-adding bins  The flat background simulation data was fed through the procedure described in section 3. 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 power spectral line in the single power spectrum was weighted accordingly [17,18] using the axion signal likelihood shown in figure 1. It was further normalized by the corresponding noise fluctuation, which was also weighted according to the axion signal shape [17,18]. With a signal window of 5000 Hz or, equivalently, 10 co-adding bins, figure 3 shows two distributions of the normalized power excess from a particular frequency where we put simulated axion signals, with (circles) and without (rectangles) the signal weighting, respectively. Figure 3 also shows the Gaussian fit results, including the means (µ) and widths (σ) of the distributions. Both means follow the predicted SNRs shown in figure 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 ex-JHEP04(2021)297  , and from a single simulated axion haloscope search experiment with CAPP-8TB backgrounds and axion signals at a particular frequency (inverted triangles in the 2nd and 3rd columns). Distributions in the 2nd column were obtained after subtracting the background with perfect fit, while those in the 1st and 3rd columns were with a 5-parameter fit. From top to bottom row, they are the distributions after the Step-1, Step-2, and Step-3 procedures, respectively. Lines are a Gaussian fit resulting in µ and σ with negligible fit errors of O(10 −3 ). 2 perimental data, but with a 5-parameter fit [15,16] for the background parametrization and subtraction. The same procedure except for the background subtraction using the simulation input background functions, i.e., a perfect fit, was also applied just to the simulation data to determine the inefficiency resulting from the background subtraction [15][16][17][18][19]. The normalized power excess distributions from the CAPP-8TB experimental data are shown as triangles in the 1st column of figure 4. The distributions from a single simulated axion haloscope search experiment with CAPP-8TB backgrounds and axion signals at a particular frequency are shown in the 2nd and 3rd columns of figure 4 as inverted triangles. The distributions in the 1st and 3rd columns were obtained after subtracting the background with a 5-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 5-parameter background parametrization. This has been previously observed in axion haloscope search experiments [15][16][17][18][19]. Without shrinking the distribution after the Step-3 procedure with perfect fit (bottom center of figure 4), the narrower width of the normalized power excess distribution must be induced from the background parametrization and subtraction thereafter [15][16][17][18][19].

JHEP04(2021)297
Our validations shown in figures 3 and 4 demonstrate not only our solid understanding of the analysis procedure but also the validity of our simulation data.

Incorporating the full correlations
Equation (6.1) shows the power excess in one of the frequency bins of the grand power spectrum with the signal weighting [17,18] and eq. (6.2) shows the full error propagation of the power excess shown in eq. (6.1) where P i and σ P i are the 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 the P i and σ P i weightings. With the signal window b a of 5000 Hz as shown in figure 1, L i meets the condition co-adding bins i=1 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 could mean that the noise fluctuations in the grand power spectrum σ Pg are overestimated with the terms for i = j only in eq. (6.2) or, equivalently, with co-adding bins i=1 σ 2 P i L 2 i only. The correlation coefficients ρ ij in eq. (6.2) are very likely to be negative due to overestimation of the noise fluctuations in the grand power spectrum [15][16][17][18][19].
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 [20], 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 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 can explain the narrower widths of the normalized power excess distributions shown in the bottom left (experimental data) and right (simulation data) of figure 4. 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 [20], which corresponds to the bottom left data in figure 4. After the full incorporation of the correlation matrices the triangles now follow a standard Gaussian distribution.
The results from the CAPP-8TB simulation data are also shown in (a), (b), and (c) in figure 7, respectively [20]. was also from 5000 simulated CAPP-8TB experiments. Both the inverted triangles and 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 circle distributions in figures 7(a), 7(b), and 7(c), the SNR efficiencies with a 5-parameter fit are 90.7% and 84.0% with and without the correlation matrices, respectively. The SNR of 90.7% with the full frequency-dependent correlations shows good agreement with what the HAYSTAC achieved with a single frequency-independent numerical factor [17,18] and was applied to our previous publication [20]. One also can expect such an SNR agreement between the HAYSTAC and CAPP-8TB experiments based on figure 5, where each correlation coefficient map does not depend significantly on the frequency, or is almost unaffected by the presence of the axion signals.

Improvement
Following the HAYSTAC method, the normalized power excess in the Step-3 can be expressed as Normalized grand power excess = P grand σ P grand ξ Step-3 , Step-3 is a frequency-independent scale factor to remedy the bias induced from the background parametrizations. Its value for the CAPP-8TB experiment is 0.92, as shown in figure 7(b). Another frequency-independent scale factor ξ Step-1.5 of 0.98 was used in Step-1.5 for the distributions shown in the middle left (experimental data) and right (simulation data) of figure 4, hence for the CAPP-8TB results [20], which is also induced by the background parametrizations. The result using a ξ Step-3 of 0.92 is in good agreement with that using the full correlations. Having said that, through the rest of this work we will employ frequency-independent scale factors. Equation (7.1) itself is valid whether the scale factor affects P grand or σ grand , though it was applied to the latter. First, we looked into which are affected by a 5-parameter JHEP04(2021)297  fit using the CAPP-8TB simulation data. Table 1 shows the category of the normalized grand power excess with different P grand and σ grand depending mainly on the background parametrizations. Figure 8 shows the normalized grand power excess distributions at the simulated axion signal frequency, where they are categorized according to table 1. Since we found the 5-parameter fit actually affects only P grand from the simulation study shown in table 1 and figure 8, we will correct for P grand only to improve SNR .
Second, we scrutinized the signal region and found that the 5-parameter fit little affects the axion signal power, but mainly distorts the background in the presence of the axion signals. The dashed lines in figure 9 are the P grand difference between (C) and (A) in table 1 and show no axion signal excess, which implies the axion signal power identified using the 5-parameter fit and perfect fit are the same. The solid line in figure 9 is the P grand difference between (F) and (A) in table 1 and the axion signal excess is not due to the 5-parameter fit, but due to the two frequency-independent scale factors ξ Step-1.5 and ξ Step-3 , which can result in overestimated SNRs. The difference between the dashed and solid lines in the axion signal sideband regions are also due to ξ Step-1.5 and ξ Step-3 .  We corrected for the background biased by the background parametrizations in the presence of axion signals using the CAPP-8TB simulation data and background only simulation data. The CAPP-8TB simulation data was used for the background parametrizations. The parametrizations were then applied to the background only simulation data to extract the background corrections shown as rectangles (blue) in figure 10, which are actually the power excess for the correlation coefficient calculations in section 6, after further going through the Step-3 procedure. However, the background correction using the rectangles in figure 10 cannot suppress the overestimated axion signal excess, which is shown in figures 9 and 10 by the solid line. This can result in overestimated SNRs as mentioned earlier. Hence we refer to the correction as "underestimated background correction". To avoid the overestimated axion signal excess, we applied another scale factor ζ of 0.44 to the underestimated background correction, which is shown in figure 10 as circles (red) which are referred to as "scaled background correction". The scale factor of 0.44 was obtained at the axion signal frequency by equating the difference between the two background corrections (rectangles and circles in figure 10) and the overestimated axion signal excess (solid line in figure 9), where the scaled background correction is the product of the scale factor and the underestimated background correction. The left and right plots in figure 10 show the background corrections in the frequency regions without and with the simulated axion signals, respectively. The narrower band with circles (red) in the left plot resulted from the scale factor ζ. Note that the noise power with the co-adding and signal weighting is about 6.5 zW in the signal region, thus the scaled background correction is at most 0.03% of the noise power. The underestimated background correction is at most 0.06% of the noise power, which can explain why the σ grand from the 5-parameter fit has no problems, as shown in figure 8(B). Figure 11 shows the SNRs after applying the two background cor-JHEP04(2021)297 rections, where rectangles (blue) are SNRs with the underestimated background correction and circles (red) are those with the scaled background correction. The SNR using the scaled background correction shows good agreement with that using perfect fit (solid line) in the axion signal frequency, while the SNR with the underestimated background correction is overestimated by a factor of 1/(ξ Step-1.5 ξ Step-3 ), as shown in the bottom plot of figure 11, as expected earlier. Finally, we confirmed that the scale factor ζ is a frequency-independent numerical factor from the other CAPP-8TB simulation data, which has to be. Otherwise, the direction of the background correction in this paper is undesirable for axion haloscope searches because of the unknown axion mass.

Summary
We report an approach to improve axion haloscope search analyses using the data obtained with the CAPP-8TB haloscope. By correcting for the background biased by the background parametrizations in the presence of axion signals, we realized an SNR of about 100%. Given the axion haloscope search analyses to date, the scanning rate can be improved by 21%, with about a 10% improvement in the SNR. For the CAPP-8TB results [20], the limits of g aγγ can be improved by about 5% with an improved SNR , e.g., from about 1.00×10 −14 to 0.95×10 −14 GeV −1 . This improvement is another low cost innovation in axion haloscope searches, where all the experimental parameters are currently at their best.