Exploration of the polarization angle variability of the Crab Nebula with P OLARBEAR and its application to the search for axion-like particles

The Crab Nebula, also known as Tau A

The Crab Nebula, also known as Tau A, is a polarized astronomical source at millimeter wavelengths.It has been used as a stable light source for polarization angle calibration in millimeter-wave astronomy.However, it is known that its intensity and polarization vary as a function of time at a variety of wavelengths.Thus, it is of interest to verify the stability of the millimeter-wave polarization.If detected, polarization variability may be used to better understand the dynamics of Tau A, and for understanding the validity of Tau A as a calibrator.One intriguing application of such observation is to use it for the search of axion-like particles (ALPs).Ultralight ALPs couple to photons through a Chern-Simons term, and induce a temporal oscillation in the polarization angle of linearly polarized sources.After assessing a number of systematic errors and testing for internal consistency, we evaluate the variability of the polarization angle of the Crab Nebula using 2015 and 2016 observations with the 150 GHz Polarbear instrument.We place a median 95% upper bound of polarization oscillation amplitude A < 0.065 • over the oscillation frequencies from 0.75 year −1 to 0.66 hour −1 .Assuming that no sources other than ALP are causing Tau A's polarization angle variation, that the ALP constitutes all the dark matter, and that the ALP field is a stochastic Gaussian field, this bound translates into a median 95% upper bound of ALP-photon coupling gaγγ < 2.16 × 10 −12 GeV −1 × (ma/10 −21 eV) in the mass range from 9.9 × 10 −23 eV to 7.7 × 10 −19 eV.This demonstrates that this type of analysis using bright polarized sources is as competitive as those using the polarization of cosmic microwave background in constraining ALPs.

I. INTRODUCTION
Tau A 1 is a polarized astronomical source at millimeter wavelengths.It has been used as a stable light source for polarization angle calibration in millimeter-wave astronomy [1][2][3].It is also known that its intensity varies as a function of time at a wide variety of timescales and wavelengths [4][5][6][7][8][9][10].The precise time-resolved polarimetry of Tau A is also becoming an active research topic in astrophysics [9][10][11].Thus, it is of interest to verify the stability of the polarization angle of Tau A in millimeter wavelengths.If detected, variability must be considered as the systematic error in the use of Tau A as a polarization angle calibrator, and may also be important for understanding the dynamics of Tau A.
One intriguing application of the time-resolved polarimetry of Tau A is to use it to search for axion-like particles (ALPs).Axions or ALPs, which are pseudoscalar fields coupled to photons through a Chern-Simons term, are potential candidates for the dark matter.The axion was first proposed as a pseudo-Nambu-Goldstone boson arising from a solution to the strong CP problem in quantum chromodynamics [12][13][14][15], while similar particles called ALPs are predicted by string theories [16][17][18].Although ALPs do not solve the strong CP problems, the ALPs can be bosonic ultralight dark matter and fuzzy dark matter, because their coupling to photons does not have a fixed relationship with their mass [19,20].Ultralight dark matter with masses around 10 −22 eV is especially interesting as it could also resolve small scale tensions that exist in the standard cold dark matter model [21,22].Beyond their gravitational effects, ALPs introduce birefringence which presents a unique detection pathway for the presence of ALPs [23,24].If ALPs are present, one will observe a net rotation of the polarization angle based on the change in the ALP field between where a photon was emitted and detected: where g aγγ is the ALP-photon coupling constant and ϕ is the ALP field.The search for polarization angle oscillation due to ALPs has been conducted using the light from astrophysical objects, radio galaxies [25], jets in active galaxies [26], protoplanetary disks [27], pulsars [28][29][30][31], and the cosmic microwave background (CMB) [32][33][34][35] and [36,hereafter PB23].
In this paper, we report our analysis using the observations by the Polarbear experiment, a ground-based experiment installed on the 2.5 m aperture off-axis Gregorian-Dragone type Huan Tran telescope observing the polarization of the CMB from the Atacama Desert in Chile [37,38], at 150 GHz from 2015 to 2016.This paper is organized as follows.Section II gives the brief review of the phenomenology of the ultralight dark matter and how the polarization oscillation signal appears for the case of Tau A. Section III describes the procedure of data analysis to calibrate and evaluate polarization angle of Tau A and its variation.Section IV describes the data validation methods.Section V describes the results.Section VI describes dedicated studies for the evaluation of systematic errors.Discussion and conclusion are described in Sec.VII and Sec.VIII respectively.

II. MODELING OF ALP-INDUCED SIGNAL
This section discusses the statistics of ultralight dark matter in general and its phenomenology in the case of ultralight ALPs through the weak interaction with the photons, and models the expected ALP-induced polarization oscillation signal of Tau A. We model the ALPinduced polarization signal assuming that the ALP with a single mass constitutes all the dark matter.

A. Statistics of the ultralight dark matter
During the formation of the Milky Way galaxy, the dark matter relaxes into gravitational potential wells and forms dark matter halos with a Maxwellian velocity distribution of characteristic dispersion velocity, v vir ∼10 −3 c, described by the standard halo model (SHM) [39,40].
We consider the statistical properties of such an ultralight dark matter field, called a virialized ultralight field (VULF) [41][42][43].Neglecting the velocity of the dark matter field, the ultralight dark matter field oscillates at its Compton frequency (= m a c 2 ℏ −1 ), where m a is the mass of dark matter.The velocity distribution of the SHM leads to the spectral broadening of the oscillation frequency characterized by the coherence time τ coh (= m a v 2 vir ℏ −1 ), and also leads to the spatial dispersion characterized by the coherence length λ coh (= ℏm −1 a v −1 vir ).In the limit that the observation time is much shorter than the coherence time of the dark matter, the time variation of the dark matter field is described as a single mode sinusoid as where ϕ c and ϕ s are real-valued amplitudes.Since the field modes of different frequencies and the random phase where ϕ 0 = ϕ 2 c + ϕ 2 s follows a Rayleigh distribution with scale factor of ϕ DM / √ 2, and θ = arctan(ϕ s /ϕ c ) follows uniform distribution from 0 to 2π.The ϕ DM is determined by the average local dark matter density ρ DM ∼0.3 GeV/cm 3 [44], and thus

B. Axion-like polarization oscillation of Tau A
We discuss an axion-like polarization oscillation signal for the case of Tau A. Tau A is 2000 ± 500 pc away from the Earth [45], and its diameter is about 1.7 pc [46].As shown in Fig. 1, the distance to Tau A is much longer than the coherence length of an ALP across the mass scales from 9.9 × 10 −23 eV to 7.7 × 10 −19 eV, which is investigated in this study.Thus the oscillation amplitudes of the local field and the field at Tau A are independent.Since the size of Tau A is much larger than the Compton wavelength (= ℏm −1 a c −1 ) of an ALP, the oscillation at Tau A averages out, because various phases of oscillation at Tau A contributes to the signal [See Eq. (2.18) of 47].Therefore, the search of an axion-like polarization oscillation by Tau A is only sensitive to the local ALP field, and the signal (Eq.( 1)) may be modeled as where ψ 0 is the polarization angle of Tau A, and ϕ local (t) is the local stochastic ALP field, which oscillates according to Eq. ( 2).

III. ANALYSIS
The data processing and map-making pipeline follows Polarbear Collaboration [48, hereafter PB20] with the improved of half-wave plate (HWP) angle estimation introduced in the improved result [49].The data set of this study was analyzed for the calibration of PB20.In this study, we re-analyze the same data set adopting a blind-analysis framework [50].The polarization angle of Tau A in each observation is not seen (blinded) until the validation of the data analysis is completed (Sec.IV), i.e. the criteria for the internal consistency test (Sec.IV B) are met, and the possible systematic errors are evaluated (Secs.IV C and VI).The polarization maps and the absolute angle of Tau A (ψ 0 ) averaged over the entire observation period are not blinded.

A. Observations
In the 4th and 5th seasons, over 486 days from 2015 to 2016, Polarbear observed Tau A several times per week, for a total of 220 observations.The Polarbear receiver employs a total of 1274 transition-edge sensor (TES) bolometers cooled to 0.3 K observing the sky through lenslet-coupled dipole antennas [51].While Polarbear has a beam with a 3.6 ′ FWHM, Tau A has an angular size of 7 ′ × 5 ′ [52].Therefore, we do not examine the detailed structure of Tau A. Figure 2 shows the map of Tau A observed by Polarbear.Throughout this study, the polarization angle of Tau A is evaluated by integrating the map around Tau A within 12 ′ diameter (Sec.III E).
Tau A is at declination of 22  .In our observing schedule, Polarbear always observes Tau A when it sets, from 39 • to 30 • in elevation, using a same raster scan pattern.The telescope scans 5.5 • back-and-forth in azimuth at 0.2 • /s, while the elevation continuously tracks Tau A. Between each sweep in azimuth, called a subscan, the elevation changes by 2 ′ , about half of the beam width.The typical observation time for each observation is one hour.Before and after every observation, a chopped thermal source calibration is carried out to characterize the relative detector gain, gain drift, and detector time constant.These values can be different for each observation depending on how each TES was tuned which depends on optical conditions of the atmosphere.

B. Data selection
We apply the same data selection criteria as PB20 except for the following optimization of the data selection conditions.The data volume is summarized in Table I.Observation efficiency is low because the primary science target of Polarbear is the CMB, not Tau A. Since Polarbear has a 2.4 • field of view, the Tau A signal is contained in less than 1% of the entire data.The remaining data outside of 12 ′ diameter around Tau A are used for evaluating the data quality.The data selection efficiency is summarized in Table II.After the data selection, the number of observations is reduced from 220 to 95.
The data selection criteria for Tau A observations differ from PB20 in three ways.First, the threshold for weather conditions is tightened to average precipitable water vapor (PWV) < 2.5 mm, and newly introduced criterion allowing data only when max PWV − min PWV < 0.6 mm.These criteria are motivated by the estimated intensity to polarization leakage (Sec.III D 2).Second, the threshold for common mode power spectral density (PSD) is relaxed, because the common mode atmospheric signal is larger in Tau A observations than in CMB observations.During the observation, the telescope tracks Tau A and changes elevation, CMB observations are per-formed with the constant elevation scans (CES).Third, the Sun-Moon avoidance cut is relaxed to allow Tau A observations more than 20 • away from either source, because the observation patch for Tau A is smaller than the CMB patch in PB20.

C. Time domain processing and demodulation
This section provides a brief review of the processing of the timestreams, measured by the TESs in the Polarbear receiver.Each TES bolometer is coupled to a dipole-slot antenna and measures a single polarization of incident light.Polarbear employs a half-wave plate (HWP) at the prime focus which continuously rotates at 2 Hz.Thus, the idealistic timesteams of the detector are modulated as follows [53, hereafter T17]: where I(t), Q(t) and U (t) are the incident Stokes parameters defined in telescope coordinate, χ(t) is the rotation angle of the HWP, and is the modulation function.We obtain the intensity signal d 0 (t) = I(t) by applying a low-pass filter, and also we demodulate the timestreams and extract the polarization signal as d ideal In the presence of the static instrumental polarization (A 4 ), the mis-estimation of the HWP angle (∆χ), the polarization angle of detector (θ det ), the intensity to polarization (I2P) leakage coefficient (λ 4 ), the time constant of detector (τ ), and the timing offset (∆t) the modulated timestreams is modified as (8) where ∆I(t) = I(t) − ⟨I(t)⟩ t is the drift of the intensity signal, and the modified modulation function m ′ (χ) is Therefore, the demodulated timestreams are modified as  where ∆τ , ∆θ det and ∆λ 4 are the mis-estimation of τ , θ det and λ 4 , respectively.The argument of complex phase of Eq. ( 11) corresponds to the mis-estimation of polarization angle in telescope coordinate, therefore the corresponding mis-estimation of the polarization angle of Tau A is expressed as ∆ψ = 2 (∆χ + χ∆t + χ∆τ ) + ∆θ det .
Also, the mis-estimation of the I2P leakage coefficient ∆λ 4 biases the Tau A polarization angle as where Q Tau A , U Tau A and I Tau A are the Stokes parameters of the incident light from Tau A. This bias is approximately expressed as where p frac is a complex polarization fraction of Tau A defined as p frac I Tau A = Q Tau A + iU Tau A .Each of these calibrations and imperfections are characterized as we describe in the following section.

D. Calibration
This section highlights the improved polarization angle calibration methods in this study.The calibration of the pointing model of the telescope, pointing offsets of detectors, effective beam function, relative gain variations, detector time constants, and polarization efficiencies are the same as those used in PB20.

Relative polarization angle between detectors
The relative polarization angle between detectors is calibrated from the average of Tau A polarization angle (ψ 0 ) over the entire observation period for each detector in the same way as PB20, with all the improvements of the I2P evaluation and the calibration of relative polarization angle between observations described in the following sections.The calibration of relative polarization angle between detectors is performed iteratively until it converges.Typically two iterations are required.The uncertainty of the relative polarization angle between detectors is estimated to be 0.1 • for each detector, which is negligibly small as we discuss in Sec VI H.

Intensity to polarization leakage
The estimation of the I2P leakage is updated from PB20.To validate the I2P leakage estimates, two different I2P leakage estimates are performed on Jupiter observations and the results are compared.There are two major sources of I2P leakage, one is a leakage due to the imperfection of optical components and the other is due to the nonlinearity of the detectors (T17).The I2P leakage estimated by two different methods includes both effects.
The first method uses the unpolarized atomospheric 1/f signal, which closely follows the method used in PB20 developed by T17.We take the correlation between the one-hour intensity timestream and one-hour polarization timestream and obtain the I2P leakage coefficient λ 4 .This assumes that the intensity timestream is dominated by the signal of unpolarized atmospheric fluctuation, and the fluctuation of the polarization timestream is dominated by the intensity timestream leaked into polarization by the I2P leakage effect.Tau A, or Jupiter for validation, is masked out from the timestream so that the Tau A (Jupiter) signal does not bias the estimation of λ 4 .
The second method uses Jupiter signal.The I2P can be decomposed into monopole, dipole, and quadrupole terms [54].We measure the monopole term by integrat- ing the map around Jupiter within 12 ′ diameter as where I p , Q p and U p are the p-th I, Q and U map pixel.The contamination from higher order multipoles is found to be small [T17].Here we assume that the polarization of Jupiter is purely from I2P from our instruments.
Jupiter is slightly polarized due to synchrotron radiation [55,56], and its polarization at 150 GHz is expected to be ≪ 1%.
Figure 3 shows the comparison of the observation-byobservation λ 4 by the two methods.The λ 4 of method 1 shows larger PWV dependence than that of method 2, and the λ 4 at larger PWV show closer values to that of method 2. We attribute the PWV dependence of method 1 to bias, based on the verification of the systematic errors discussed in Sec IV C.
The I2P estimates by method 1 may be biased, because at lower PWV, the 1/f fluctuation from non-optical effects, such as the fluctuation of the focal plane temperature or the temperature of the primary or secondary mirror, the warm readout electronics may not be negligible [57] 2 .On the other hand, method 2 is robust to any 1/f fluctuations because it is estimated by the point source, Jupiter.However, method 2 also has its drawbacks, such as the fact that they are observed under different observing conditions and that Jupiter is about 10 times brighter (T Jupiter ∼2 K) than Tau A, so it may be more affected by the detector non-linearities and have constant bias.
In our Tau A analysis, we construct the I2P leakage model by making a modification to method 1 and applying it to the Tau A map-making process.We model the PWV dependence of method 1 as a linear function for each detector and substitute the median PWV (=0.7 mm) of Tau A observations, to make representative λ 4 for each detector.We call this method 1 ′ .Figure 4 shows the comparison of detector-by-detector I2P coefficient λ 4 by method 1 ′ and method 2. They show reasonable agreement.We discuss the uncertainty of I2P leakage estimates by the difference of method 1 ′ and 2 in Sec.VI.

Relative polarization angle between observations
As we have shown in Eq. ( 12), a mis-calibration of the relative polarization angle between observations may be caused by the mis-estimation of the HWP angle (∆χ), the time constant of detectors (∆τ ), and the timing offset (∆t).In PB20, we calibrated the time constant of each detector in each observation by the average of the chopped thermal source calibration right before and after the Tau A observations, and the achieved polarization angle error was sufficiently small.
In the search for polarization oscillation, the impor-tance of the calibration of relative polarization angle between observations increases compared to the analysis which takes an average of all observations, e.g.CMB power spectrum analysis in PB20.In this study, we further calibrate the relative polarization angle between observations assuming that the telescope and the polarization angle of the instrumental polarization due to the primary mirror ( 1 2 arg(A 4 )) are stable. 3Because our HWP modulator is the second optical element in the path of light from the sky to the focal plane [T17], the reflection at the primary mirror produces the instrumental polarization [58].This becomes a good calibrator for the relative polarization angle between observations, because it illuminates all detectors stably and uniformly (|A 4 |∼0.1 K, T17).The instrumental polarization is estimated as a signal synchronous to the 4th harmonic of the rotation angle of the HWP, by averaging the modulated timestreams (Eq.( 8)) over the HWP angle for each observation [T17], masking Tau A signal and glitches in the time domain.
The amount of the corrected polarization angle variation between observations is 0.16 • root mean square (RMS).This variation is measured within the statistical uncertainty of the instrumental polarization angle 0.02 • , which is a median of the statistical uncertainty per observation in terms of the standard deviation (STD).

E. Estimation of polarization angle of Tau A and its statistical uncertainty
The polarization angle of Tau A (ψ) is estimated in the map-domain for each observation, with the same filteredbinned map-making as PB20.We estimate ψ by integrating over the map around Tau A within 12 ′ diameter as where d t is the measured ψ, and t is the average observation time over one hour.The diameter of the integration region is chosen so that the region covers the size of Tau A convolved by our beam size (FHWM 3.6 ′ ) including the measured boresight pointing drift of our telescope (Sec.VI F).This angle estimation method is not the most efficient method, but it is immune to the complexity of the structure of Tau A, and robust to systematic errors.The statistical uncertainty of each observation σ t was estimated by the bootstrap method, also called "signflip noise estimation," which is widely used for noise estimation in the CMB analysis.First, we generate the signflip noise only map for each observation by randomly assigning a +1 or −1 factor to each detector and coadding individual detector maps 4 .The random sign is assigned so that the total data weight is even for both signs.The signflip map will cancel the Tau A signal as well as the 1/f noise which is common among all the detectors.Then, we sample the noise realizations of Stokes parameters of Q and U from signflip map over independent 48 regions with the same size as the one used to estimate ψ as δQ = noise region p One of the statistical noise realizations of Tau A polarization angle δψ is obtained by adding δQ and δU to the Tau A signal as The statistical uncertainty σ t is obtained from the STD of the bootstrapped noise realizations δψ.The typical value of σ t is 0.16 • from the median of all observations.

F. Estimation of Tau A polarization angle spectrum
As discussed in Sec.III A, the Tau A observations span over 486 days, and are performed daily.Therefore, the Nyquist frequency (f NYQ ) is well-defined in our data set.Because of the aliasing, the frequency bins above f NYQ are strongly correlated with the frequency bins below f NYQ , and almost all degrees of freedom in our data are contained in the following set of 242 bins: where ∆T ≃ 1 − 1/366.24day is the observation interval, one local sidereal day.The harmonics of f NYQ (= 1/(2∆T )) are not included, because the harmonics of f NYQ correspond to the constant mode and we do not have sensitivity to them.Over the 486 days of the observation period, we have 95 observations of Tau A. Therefore, the number of frequency bins (242) is larger than the effective number of degrees of freedom (95).This induces a bin-by-bin correlations among the frequency bins below f NYQ .We also accept the bin-by-bin correlations due to aliasing, and the scientific result is extended to The maximum frequency to which we have sensitivity is determined by the duration of a single Tau A observation, one hour (Sec.III A).
We estimate the frequency spectrum by least-square spectral analysis [59], similar to Fourier analysis.Hereafter we call the amplitude of least-square spectral analysis "LSSA."The LSSA of the measured Tau A angle d t with the error bar σ t is df = argmax where argmax A∈C [f (A)] represents a complex amplitude that maximizes a function f (A).The solution of Eq. ( 21) is where c t = cos(2πf t), s t = sin(2πf t).Hereafter, we express the complex amplitude of the LSSA as where Similarly, in the following we express the LSSA of the ALP oscillation and noise as φf = t F f t ϕ t and ñf = t F f t n t , respectively.Note that φf and ñf are defined using the same noise weight, F f t .
The LSSA estimate of the signal amplitude at each frequency is biased for multiple reasons.We quantify the bias on the oscillation amplitude averaging over the phase of the oscillation by the transfer function F f , the detailed derivation of which is described in Appendix A. Figure 5 represents the individual transfer functions and the overall transfer function.The overall transfer function is evaluated by simulating all the transfer functions at once.

G. Estimation of polarization oscillation amplitude
We search for a single mode polarization oscillation without knowing its frequency and phase.Therefore, we take the convolution for all possible phases, whose probability density function is flat, and estimate the amplitude A obs f assuming the presence of a signal at each frequency.We maximize the following likelihood function: where A f and θ f the real-valued true amplitude and phase of the oscillation signal at f , and P (θ f ) is the uniform probability density.
is the probability density for obtaining df for given signal parameters, A f and θ f .
is the difference of the LSSA of data minus signal, and Ñf is a noise covariance matrix with elements where ⟨•⟩ represents average over signflip noise simulations.Note that when the noise covariance is diagonal and 25) is simplified as where I 0 is the modified Bessel function of the first kind.The equivalent likelihood is obtained by considering that 2| df To evaluate the detection significance, we form the following test statistic, which quantifies the global significance of the preference for the signal over the null hypothesis: where is the test statistic for the local significance at each frequency.A best f and θ best are the parameters that maximize Eq. ( 26).

H. Estimation of ALP-photon coupling
We assume the stochastic ALP field (Eq.( 2)) to estimate the ALP-photon coupling, g aγγ .We consider the estimator for the square of the ALP photon coupling g 2 aγγ as and translate it into g aγγ .We treat ĝ2 aγγ as a single parameter and allow it to be negative.The probability density for ĝ2 aγγ is constructed by Monte Carlo simulations, using the probability density for obtaining df where Cf is a covariance matrix of the noise and signal The signal covariance matrix is Note that Eq. ( 33) is obtained by the multivariate Gaussian convolution of the probability density for the local stochastic ALP amplitude and the probability density for obtaining df for given signal parameters (Eq.( 26)), where we relate the stochastic ALP amplitude φf with A f exp(iθ f ) of Eq. ( 27).
It is immediately shown that ĝ2 aγγ is unbiased because the variance of df is As discussed in Appendix B, the efficiency (variance) of this estimator is comparable to the estimator used in PB23, while this estimator is computationally faster.

IV. DATA VALIDATION
To validate the data analysis, we test the internal consistency of the data using a blind analysis framework.Hereafter we call this test as a null test.We follow the similar formalism used in PB23.We split the data in half and form null statistics by taking a difference to cancel the oscillation signal but amplify the systematic errors, and analyze the consistency of the null statistics with the null hypothesis.The data-splits are largely categorized into two types of splits to probe for different type of systematics.
Bolometer-splits: These splits probe the miscalibration between detectors or the systematic error in each observation such as the time drift of the detector characteristics over one hour.
Observation-splits: These splits probe the systematic variations between observations.
In addition to the data-splits considered in PB20, we introduce the following data splits dedicated to probe the systematics of the polarization angle.
1. 2f/4f angle by obs/bolo: These split the data into observations or detectors by angle of signal synchronous with the rotation of the HWP, to probe the detector-by-detector or observation-byobservation mis-calibration of the HWP angle.2. Time constant by bolo: This splits the data into detectors with larger time constants and those with smaller time constants, to probe the mis-calibration of detector's time constants.

A. Null statistics
We construct three different types of null statistics.First, we take noise-weighted differences as where d t,k and σ t,k are the k-th split data and its error bar.For the data-splits which split the data periodically, we did not use χ DC null for null test because of the risk of unblinding the signal (Table.III).The other null statistics are constructed in the frequency domain, by the difference of LSSA for each data-split as where df,k = dreal f,k + i dimag f,k is the LSSA of k-th split data (Eq.( 23)) and Dmn is the covariance matrix of the difference of LSSAs.χ DC null and χ AC null follows the chi-square distribution with one degree of freedom, χ 2 null approximately follows the chi-square distribution with two degrees of freedom.
As discussed in Sec.III F, the null statistics for different frequency bins are correlated.The pattern of correlation differs for the observation-splits, which introduces imperfect signal cancellation in the null statistics for the observation-splits.This is especially problematic when large sinusoidal signals are present, but is sufficiently small in this study since we are verifying the presence of a sinusoidal signal, which is known to be small.For this reason, we do not apply the transfer function (Sec.III F) for the construction of null statistics (Eq.( 41)).

B. Null test results
Table.III shows the null statistics for each data-splits.The null statistics are compared with the 1440 signflip noise only simulations 5 (Sec.III E).Then we evaluate the probability to exceed value (PTE) by counting the number of simulation whose null statistics exceeds the real data.Table III shows the distribution of PTE for each data split, and Fig. 6 shows the distribution of PTE for each frequency bin.
Next, we compute six representative test statistics for bolometer-splits and observation-splits separately: (1) the average of χ DC null among all data-splits, (2) the average of χ AC null among all data-splits, (3) the most extreme total 5 48 samples from each of the 30 signflip maps.χ DC null with † were blinded during the data validation process, because these partially unblind oscillation signal at certain frequencies.
FIG. 6. PTEs of null statistics for individual frequency bins.The subscript "sp" denotes "data-splits."Although there are correlations between frequency bins, there is no significant deviation from a uniform distribution.χ 2 null by data-splits summed over all frequencies, (4) the most extreme total χ 2 null by all frequencies summed over all data-splits, (5) the most extreme total χ 2 null among all frequencies and all data-splits, (6) the total χ 2 null summed over all frequencies and all data-splits.( 1) and ( 2) probe for the biases, (3)-( 5) probe for the outliers specific to the data-splits or frequencies, and ( 6) probes for the misestimation of our uncertainties.The PTEs for the representative test statistics are also computed by comparison with noise only simulations.
As a pass criteria of the null test, we require the PTE of the lowest PTE value to be larger than 5% to probe the biases and outliers.We also require the PTE of the highest PTE value6 to be larger than 5% to check the mismatch between the real and estimated uncertainties. 7he numerical values for these statistics are given in Table IV.Both bolometer-split and observation-split pass the criteria.

C. Correlation with systematic template
In order to detect possible systematic errors, we additionally test if any significant correlation coefficient is found between the measured Tau A polarization angle and the variation of instrumental conditions or environmental conditions, which could potentially produce systematic time-dependent polarization angle fluctuation.We call the representative fluctuation of instrumental or environmental conditions used in this systematic test as systematic templates.The systematic templates are constructed using the Polarbear's instrumental or environmental monitors.We perform this systematic test using the following systematic templates: the average of PWV, the ambient temperature, the focal plane stage temperature, and the polarization angle of the mirror polarization. 8e define the weighted correlation coefficient between data d t with error bar σ t and the systematic template where is the weighted covariance, and are the weighted average.Figure 7 shows the significance of the correlation coefficients obtained by comparing the correlation coefficients for 1440 noise realizations, and no significant correlation coefficient is found.Earlier in the data analysis, this method was beneficial to identify systematics and determine the optimum calibration and analysis strategies.This method led to the detection of the PWV-dependent mis-estimation of the I2P leakage and the mis-calibration of the polarization angle in the telescope coordinate with correlation coefficients 0.51 and 0.55, respectively, resulting in the improvement of the angle calibration discussed in Sec III D.

V. RESULTS
We report the results of Tau A polarization angle, its frequency spectrum, and upper bounds of its oscillation amplitude.We also report the upper bounds of ALPphoton coupling, assuming that no sources other than ALP are causing Tau A's polarization angle variation.
Our estimated timestream of Tau A polarization angle is shown in the top panel of Fig. 8.The black dots in the bottom panel of Fig. 8 shows its LSSA at the frequency bins below f NYQ (Eq.( 19)).To evaluate the detection significance, test statistics were computed using a total of 24,000 noise realizations generated based on the bootstrap method (Sec.III E) to have sufficient accuracy of the significance level up to 3σ.The blue lines and red lines in the bottom panel of Fig. 8 show the approximated amplitudes with certain local and global significance levels. 9Figure 9 shows the distribution of test statistics.The largest significance level we found is p-value of 0.50%, corresponding to 2.5σ, at 1/61 day −1 .
The black solid line in Fig. 8 shows the LSSA with 10 times finer frequency bins than Eq (19).Another peak at 1/52 day −1 , which is not captured by the sparse frequency bins, is found to have a global significance level of 2.5σ by calculating the test statistics over the 10 times finer frequency bins.We confirmed that the null-test pass criteria were still met with 10 times finer frequency bins.The local significance levels of these two most significant oscillation modes at 1/61 day −1 and 1/52 day −1 are found to be 4.4σ and 4.5σ analytically.We note that the local significance is not relevant as a detection significance because we are examining many frequency bins.
The orange solid line in Fig. 8 shows the LSSA of data minus the best fit sine curve at 1/52 day −1 .The data minus best fit is consistent with the null hypothesis, with the global significance level smaller than 1σ.This means that the two most significant oscillation modes are correlated, and the result can be interpreted as a hint of a single-mode oscillation signal.
If we interpret these two most significant oscillations as axion-like polarization signals, then the corresponding ALP masses are 7.8 × 10 −22 eV and 9.2 × 10 −22 eV, respectively.The possible interpretations of the hint of a signal are discussed in Sec.VII.
We report 95% upper bounds on the oscillation amplitudes and g aγγ , because no point exceeds the significance level of 3σ.The 95% upper bounds on the oscillation amplitudes at the frequency bins below f NYQ are shown in the top panel of Fig. 10.As a measure of our experimental sensitivity, we report the median upper bound below f NYQ as 0.065 • .The upper bound A 95% is obtained with the frequentist approach by the classic Neyman construction [61] as where Â is the maximum likelihood estimator of A f using Eq. ( 25).P ( Â|A) is the probability density for Â, which is obtained by Monte Carlo simulations of Â using Eqs.( 25) and (26).The bottom panel of Fig. 10 shows the 95% upper bounds extended to frequency bins larger than f NYQ (Eq.( 20)).Due to aliasing, the signal seen at f (< f NYQ ) may be aliased from higher frequencies.The red dashed curve represents smoothed approximation of our upper bounds where sinc function approximates the transfer function F f (Eq.(A9)).The 95% upper bounds of g aγγ assuming the stochastic ALP field are also constructed using the Neyman construction.We consider the constraint on the square of the effective amplitude of polarization oscillation, g 2 ≡ g 2 aγγ F 2 f ϕ 2 DM /4, and convert it to a constraint on g aγγ .The detailed derivation of it is described in Appendix C. 10 The median 95% upper bound on the effective oscillation amplitudes below f NYQ is g 2,95% = 0.13 • .Assuming that ALP constitutes all the local dark matter (κ = 1, ρ 0 = 0.3 GeV/cm 3 ), the smoothed approximation of our stochastic bounds is Our bounds are shown in Fig. 11 along with the other selected constraints.Our bounds for individual frequencies and the smoothed approximation are shown using red and black curves, respectively.The black dashed line shows the 95% upper bounds assuming that the amplitude of the ALP field is deterministically given by the density of local dark matter, so called "deterministic 10 The stochastic bounds obtained by this method are 2.2 times conservative than the deterministic bounds using Bayesian inference with a flat prior.The scaling between these two are consistent with our previous result demonstrated in PB23.This section presents the dedicated studies for the evaluation of possible systematic errors.Table V shows the median statistical uncertainty per observation, and the RMS of the systematics, the systematic time-dependent polarization angle fluctuations.The frequency distribution of systematic polarization angle fluctuation is estimated by its LSSA.The overall systematics as well as the individual systematics are shown in Fig. 12.The systematics are shown up to f NYQ , because systematics and statistical error scales the same at frequencies higher than f NYQ .The correlation coefficient between different systematics are found to be smaller than 0.3, as shown in the right panel of Fig. 12.The overall systematics are formed by quadrature sum of individual systematics.We find the total systematics to be subdominant compared to the statistical error at all frequencies.Possible residual systematic errors are discussed in Sec VII.

A. Polarization angle calibration
The uncertainty of polarization angle calibration in the telescope coordinate produces additional statistical uncertainty in the measured Tau A polarization angle.The overall uncertainty of the polarization angle calibration in telescope coordinate is given by the statistical error of the instrumental polarization angle per observation, 0.02 • STD (Sec.III D 3).This uncertainty degenerately includes time constant, timing offset, and HWP angle uncertainties.The upper bound of the uncertainty due to the time constant and timing offset is estimated independently for the cross check from the chopped thermal source calibration.These are found to be < 0.04 • and < 0.03 • , respectively.

B. Intensity to polarization leakage
Mis-estimation of the I2P leakage can produce systematic polarization angle fluctuation.Since the average polarization fraction of Tau A at 150 GHz is 7%, the uncertainty of the I2P leakage coefficient of 0.02% per observation leads to a systematic error of the Tau A polarization angle of ∼0.1 • (Eq.( 14)).The systematics due to the uncertainty of I2P leakage was estimated by the difference of Tau A angle when different models of I2P (Sec.III D 2) are applied.The possible polarization angle fluctuation due to I2P mis-estimation is 0.06 • RMS.

C. Ground contamination
The systematic contamination of data synchronous with the telescope's pointing in ground coordinates is called ground contamination.During Tau A observations, the telescope tracks Tau A by changing azimuth and elevation simultaneously (Sec.III A).Due to the complexity of the scan strategy, the subtraction of the ground contamination is not performed, in contrast to the CMB analysis.Here we evaluate the amount of the ground contamination on the polarization angle of Tau A. First, we generate a ground polarization map for each observation by stacking the timestreams in the ground coordinate, masking timestreams around Tau A FIG. 11. 95% Upper bound of gaγγ.The upper bounds assuming a stochastic ALP field is shown using solid lines.The orange solid curve represents a constraint by searching for an axion-like oscillation from CMB [PB23].The blue solid curve represents a constraint by searching for an axion-like oscillation from various pulsars and Tau A [30].The red solid curve represents our stochastic bounds, the primary results of this study.The black curve represents the smoothed approximation of our stochastic bounds (Eq.( 47)).The upper bounds assuming deterministic ALP field are shown using dashed curves.The green dashed line represents a constraint by searching for an axion-like oscillation of CMB [35].The black dashed curve represents the smoothed approximation of our deterministic bounds (Eq.( 48)).The green shaded bound represents a constraint from the absence of the suppression of CMB polarization due to an axion-like oscillation in the recombination era [32].The upper bound from the CAST experiment [62], an absence of the gamma-ray excess of SN1987A [63], and an absence of the X-ray spectral distortions of the quasar H1821+643 [64] are shown.The lower bound on the ALP mass by the Lyman-α forest [65] and the Milky Way satellite galaxies [66] are also shown.
within 12 ′ diameter.A hit map of Tau A in ground coordinate for each observation is generated simultaneously from the masked timestreams.The maximum ground polarization signal is found to be about 100 µK, which is about 200 times smaller than the polarization intensity of Tau A (≃20 mK).The polarization angle bias due to the ground contamination is evaluated by adding the bias of  The possible systematic polarization angle fluctuation due to the ground contamination is 0.03 • RMS.

D. Residual 1/f
The 1/f noise also produces systematic polarization angle fluctuation, because the statistical uncertainty for each observation is estimated by the detector-by-detector signflip noise estimation, which cancels the 1/f noise, a fraction of which can be common between detectors (Sec.III E).This 1/f noise is modeled by stacking all detectors' timestreams in time domain, masking around Tau A within 12 ′ diameter for each observation.The stacked timestream is projected to the sky using the pointing of each detector to make a residual 1/f map.The polarization angle bias due to the 1/f noise is evaluated by adding the residual 1/f noise map to the Tau A signal as The systematic polarization angle fluctuation due to the residual 1/f is 0.01 • RMS.

E. Mizuguchi-Dragone breaking
The Huan Tran telescope is an off-center Gregorian design which satisfies Mizuguchi-Dragone (MD) condition to cancel the cross-polarization at the primary and secondary mirrors [67,68].The systematic contamination due to the cross-polarization from breaking the MD condition by the HWP being located at the prime focus of the telescope is expected to be negligibly small.Matsuda et al.

F. Pointing
The boresight pointing of the telescope propagates to the polarization angle of Tau A through the transformation of polarization angle from telescope coordinate to sky coordinate.The observation-by-observation fluctuation of the boresight pointing is evaluated by the variation of the celestial position of the Tau A intensity signal, which is 0.1 ′ STD and 0.3 ′ STD in right ascension and declination, respectively.The possible systematic polarization angle fluctuation due to the boresight pointing error is 0.006 • RMS.

G. Filtering & map-making bias
Since we adopt the filtered-binned map-making method, the Tau A angle estimated in map domain is biased (Sec.III E).We evaluated this bias by the full simulation of Tau A signal timestreams.The input Tau A model was constructed by convolving the Tau A map ob-served by IRAM [1] with the beam of Polarbear.We simulated timestreams of Tau A signal scanning the input map, modulating the signal by rotating HWP, and applied the same analysis as the real data to estimate biases.The systematic polarization angle fluctuation due to the filtering and map-making is 0.005 • RMS.

H. Polarization angle of detector
The mis-calibration of the polarization angle of detector averaged over the entire observation period can produce systematic polarization angle fluctuation because the data selection is different for each observation.The amount of mis-calibration of the polarization angle of each detector, calibrated over the entire observation period, is estimated to be 0.1 • RMS, by the variation of the Tau A polarization angle measured by each detector.The possible systematic polarization angle fluctuation is estimated by the average mis-calibration of polarization angle of detectors.The estimated systematic polarization angle fluctuation due to the polarization angle of detector is 0.002 • RMS.

VII. DISCUSSION
Although no point exceeds the significance level of 3σ, the hint of 2.5σ oscillation with an amplitude of 0.11 • and 1/61 day −1 and 1/52 day −1 period deserves further discussion.In this section, we enumerate possible causes and describe our investigations.

A. Statistical Uncertainties
The estimation of statistical uncertainty is verified within the fractional 1σ uncertainty of 6.5% by one of the representative test statistic of the null test, the total χ 2 null summed over all frequencies and all data-splits.The total χ 2 null scales inversely proportional to the square of the assumed statistical uncertainty.We estimated the accuracy of the statistical uncertainty by comparing the measured total χ 2 null and the simulated total χ 2 null .If we assume uniformly larger statistical errors of 6.5% for all 95 observations, the significance level of the hint of signal decreases from 2.5σ to 2.1σ.Therefore it is possible that the statistical fluctuation of the estimated statistical uncertainty is the cause of the hint of a signal.

B. Residual systematic errors
This section enumerates the investigations of the residual systematic errors.We additionally explore the possible systematic polarization fluctuation using a similar method as Sec IV C. We consider additional systematic templates: HWP speed, wind speed, humidity, UV level, As a conservative estimate, we construct the 95% upper bound of the proportionality coefficient as where n l is the l-th noise realization based on the bootstrap method (Sec.III E). Figure 13 shows the spectrum of the possible systematics modeled as ĉ95% T compared with the 1σ upper bound of statistical uncertainty.No significant possible systematics are found around the frequency bins where we find the hint of a signal.Possible systematics are larger than the 1σ upper bound of statistical uncertainty in some frequency bins, however this does not necessarily mean the actual systematics.This is because these are conservative estimates, and their accuracies are limited by the statistical error of the polarization angle measurements of Tau A. Another way to probe the systematics is to check the season-by-season consistency of the hint of a signal as shown in Fig. 14.The hint of a signal at 1/61 day −1 or 1/52 day −1 is consistently seen in both the 4th and 5th season.Additionally, we perform the null tests with 10 times higher resolution frequency bins and also by limiting the frequency bins to the two most significant bins.No residual systematics with a significance level larger than 2σ are found.Finally, we note that we did not do any on-site maintenance activities on 52 to 61 day cycles.
Another possibility of residual systematic error is the aliasing of the diurnal variation of the instrumental systematic error, such as the instrument's rotation around the boresight.Figure 15 shows the measured polarization angle of Tau A as a function of local time at the observation site.The blue points are the unbinned data, the red points are binned by 1 hour grid and the black solid line and dashed line are the fitted 6 day −1 and 7 day −1 diurnal variation, respectively.A possible explanation is that the receiver could be effectively rotated around the boresight by a small amount due to the deformation of the receiver mounting structure caused by the thermal expansion depending on the sunlight exposure, and the measured polarization angle of Tau A varies diurnally.The systematic errors of the boresight pointing due to various environmental conditions including sunlight exposure are evaluated to be smaller than 1 ′ RMS [70], which is negligible for this study.Also, only about half of the observed time is directly influenced by the sunlight because the typical sunrise time at the observation site is 7-8 am.Nonetheless, we do not have a monitor of the instrument's rotation around the boresight to exclude this possible systematic error.
While this may indicate a hint of possible systematic error, it is not statistically significant that the peaks at 1/61 day −1 and 1/52 day −1 can be mapped to 24-hour cycle diurnal signals (Fig. 16).The frequency resolution of our data is set to 1/486 day −1 by the entire observation period, 486 days.On the other hand, the n-th harmonics of 24-hour cycle diurnal variation are aliased to the frequencies of n/366.24day −1 because our observation interval is 1 − 1/366.24day, one local sidereal day.Thus, there is a high chance of such coincidence for any peak found below f NYQ ≃ 0.5 day −1 with the frequency resolution of 1/486 day −1 .In our case, n = 6 and 7 coincides with the 1/61 day −1 and 1/52 day −1 , respectively.
We also point out the possibility of the oscillation of the I2P leakage, which can mimic the signal.There are two known mechanisms of creating the I2P leakage.One is due to the optical elements, in which case the I2P leakage tends to be either constant or to gradually increase with time and not to oscillate.The other is detectors' non-linearity (T17), whose effect would correlate with PWV in our case and thus is investigated in Sec IV C.However, it is worth assessing model-independent constraints we can place on the I2P leakage variability.All the I2P leakage estimation methods discussed in Sec.III D 2 have different drawbacks to evaluate the time variation of I2P leakage.Method 1 may suffer from timedependent systematic errors, method 1 ′ assumes no time variation in I2P leakage, and method 2 is under different observation schedule and conditions than Tau A observations.We evaluate the possible systematic polarization fluctuation of Tau A by adding maximally allowed fluctuations of I2P leakage to the Tau A signal.According to the calibration uncertainty of method 2, the RMS of possible systematic polarization fluctuation is 0.16 • and is not negligible compared to the statistical uncertainty (Table V).Therefore, we cannot fully exclude this possibility.

C. Intrinsic variability of Tau A
This section explores the possibilities of intrinsic polarization angle variability of Tau A. Since we measure the average polarization angle of the entire Tau A, we cannot distinguish the intrinsic variability of Tau A from the ALP-induced polarization angle oscillation.Tau A primarily consists of two components, the Crab Nebula and the Crab Pulsar [46].The Crab Nebula is an expanding synchrotron nebula.The Crab Pulsar is a point source at the center of the Crab Nebula, the energy source of the Tau A system.Explaining the possible single-mode oscillation seen in our data with the intrinsic variability of these components is not very natural, because of the following two reasons.First, the time scale of the single-mode oscillation seen in our data, a period of 52 or 61 days, is not very natural for the Crab Nebula, the expanding synchrotron nebula with the size of several light years.Second, the amplitude of the polarization angle oscillation, 0.11 degrees, is too large to be explained by the Crab Pulsar.If the Crab Pulsar had a polarization variability that is 4 × 10 −3 = arctan(2 × 0.11 × π/180) of the Crab Nebula, that could explain the 0.11 degree polarization angle variability of Tau A. However, at millimeter wavelengths, it is expected to be four orders of magnitude or more darker than the Crab Nebula [71].Also, a single mode oscillation of the Crab Pulsar with a period around 52 or 61 days is not reported at any wavelength.Nevertheless, we explore the possibilities of intrinsic variability of Tau A using Tau A measurements in other experiments, because we cannot completely exclude them.
The time variability of the Tau A has been studied in a wide variety of timescales and wavelengths.To the best of our knowledge, the time variability around 52 or 61 day period has not been observed.Here we summarize the known time variation of Tau A.  4. Large variation in the polarization of Crab Pulsar and Crab Nebula is reported at X-ray wavelengths [9,10].
The polarization angle of Tau A has been measured by a number of experiments.However its time dependence has been explored only recently.The 95% upper bound of oscillation amplitude of ∼1 • was placed by Castillo et al. [30] on periods ranging from 6 × 10 −4 day −1 to 10 day −1 using data from the QUIJOTE experiment between 10 and 20 GHz [72].Tau A has been measured in more experiments for intensity than for polarization.Therefore, we qualitatively investigate whether the single-mode oscillation seen in our data exists in Tau A intensity using the public daily light-curve from the all-sky survey satellites over the observation period of this study.Figure 17 shows the LSSA of the light-curves measured in this study and various all-sky survey satellites.Since these measurements have different cadences than this study, a direct comparison is not possible, although we do not find clear single-mode peaks as seen in our data around the frequency bins where we find the hint of a signal.

D. Consistency with other ALP bounds
We discuss the consistency of our results with other ALP bounds.The consistency analysis discussed in this section is only valid if the hint of a signal we found is purely from ALP.To the best of our knowledge, the second tightest upper bound of an axion-like polarization oscillation of a local ALP field is placed by SPT-3G Collaboration [35], with the median upper bound of 0.071 • from 1/100 day −1 to 1 day −1 , while our median upper bound is 0.065 • .The largest discrepancy is found at the frequencies where we found the hint of oscillation.A more significant discrepancy is found when we translate the polarization oscillation to a constraint on ALPphoton coupling.One of the tightest upper bounds of the axion-photon coupling is placed by the absence of the suppression of CMB E-mode power spectrum due to an axion-like oscillation in the recombination era [32].
The 95.5% upper bound of 5.0 × 10 −13 GeV −1 at m a <10 × 10 −12 eV is also placed from the absence of the distortions of the X-ray spectrum of the quasar H1821+643 due to the ALP-photon coupling [64].The bottom panel of Fig. 18 indicates the significance of the discrepancy.Our bounds are complementary to these bounds because the ALP search by an axion-like oscillation is less dependent on the cosmological model or the model for the galaxy cluster magnetic fields.Several analyses have claimed to rule out the ultralight dark matter in mass region where we found a hint of a signal, m a < 2.0×10 −21 eV by [65] and m a < 2.9×10 −21 eV by [66].Since there are model dependencies in these analyses, our analysis is complementary in searching for an entirely different signal.Our analysis applies even if the ALP constitutes a fraction of the dark matter, the bounds are inversely proportional to the square root of the dark matter fraction, κ.
We also check the consistency with our previous result, the search for the polarization oscillation in the CMB PB23.The hint of a signal we find is below the noise level of PB23, and PB23 is consistent with this study.

VIII. CONCLUSION
We analyze the Tau A observations in the 4th and 5th seasons of the Polarbear experiment, and present a median 95% upper bound of A < 0.065 • on the polarization oscillation amplitude of Tau A over oscillation frequencies from 0.75 year −1 to 0.66 hour −1 .To the best of our knowledge, this work provides the most accurate measurements of the amplitude of the polarization oscillation of Tau A. Our improved analysis of Tau A data introduces new methodologies for investigating possible systematics, which allowed us to improve the polarization angle calibration method, and demonstrate a sensitive search for an axion-like polarization oscillation using the calibration data from the Polarbear groundbased CMB observatory located in the Atacama Desert of Chile.Our measurements constrain potential intrinsic time variability in Tau A, which is an important understanding of Tau A as a polarization calibrator, and may be used to better understand the dynamics of Tau A. Under the assumption that no sources other than an ALP are causing Tau A's polarization angle variation, that the ALP constitutes all the dark matter, and that the ALP field is a stochastic Gaussian field, this bound translates into a median 95% upper bound of an ALP-photon coupling g aγγ < 2.16 × 10 −12 GeV −1 × (m a /10 −21 eV) in the mass range from 9.9 × 10 −23 eV to 7.7 × 10 −19 eV.To the best of our knowledge this is the tightest bound from the measurement of the ALP-induced polarization angle oscillation.This study demonstrates this type of analysis using bright polarized sources are competitive as those using the polarization of CMB in constraining ALPs.Analysis using the polarization of CMB has the advantage of being immune to the complexity of the dynamics of the sources.Using multiple polarized sources and cross-correlating their time variations will lead to a more robust analysis.
We find a hint of a polarization oscillation signal with amplitude 0.11 • at 1/61 day −1 with a global significance level of 2.5σ, and an equally significant correlated polarization oscillation at 1/52 day −1 .We consider whether this hint of a signal is due to residual systematic errors or intrinsic variability of Tau A due to yet unknown dynamics, but so far we find no clear evidence for either of them.
There are possible residual systematic errors we do not exclude.Interpreting this as an ALP signal leads to significant tension with other constraints.We expect that future measurements with more data will not only reduce statistical errors but also allow us to better understand possible systematic errors, leading us to a more conclusive interpretation.A successor experiment, the Simons Array, employs dichroic detectors including 90 GHz and 150 GHz center band [78].Since the polarized intensity of the synchrotron radiation is proportional to ν −0.3 [79], we expect a more sensitive observation from the lower frequency band.Also, since an axion-like polarization oscillation is independent of the optical frequency of light, the control of frequency-dependent systematic errors is possible.
where i(t) is a set of detectors working in the observation performed at t.We calculate the bias from the miscalibration of the polarization angle of detectors by averaging over the phase of oscillation for each frequency as Finally, the calibration of the absolute polarization angle (ψ 0 ) also biases the LSSA as (A8) The overall transfer function F f is evaluated by simulating all the effects at once.The deviation of F f from unity primary comes from F average f , and for the smoothed results (Eqs.( 46), ( 47) and ( 48)) we approximate the overall transfer function F f as Appendix B: Estimator This section compares the methods for estimating g aγγ used in this study and PB23.In this study, we estimate g aγγ by maximizing the probability density for obtaining the LSSA of the data.We call this estimator the "frequency estimator."This is motivated by considering a simplified case where the noise covariance matrix (Eq.( 28)) is diagonal, i.e., when the timestamp is regularly spaced and the error bars are uniformly distributed.In this case, 2| df | 2 / ⟨|n f | 2 ⟩ follows chi-square distribution with two degrees of freedom and Eq. ( 33) is simplified as Eq. ( 32) is a maximum likelihood estimator for this likelihood function.
On the other hand, PB23 estimates the g aγγ by maximizing the probability density for obtaining time-domain data.We call this estimator the "full estimator."The likelihood function is is the probability density for obtaining data when the stochastic oscillation amplitudes of ALP are given.
Figure 19 compares the bias and efficiency of the "frequency estimator" and "full estimator."The "full estimator" is biased when the signal is small.The "frequency estimator" is a good unbiased estimator regardless of whether the noise covariance is diagonal or not, and the efficiency is slightly larger than "full estimator."Another advantage of the "frequency estimator" is its small computational cost.This is because "full estimator" requires convolution for constructing likelihood (Eq.(B2)), while the corresponding convolution is analytically conducted for the "frequency estimator" (Sec.III H).

Appendix C: Upper bound
This section describes the method for estimating the upper bound of g aγγ .First, we consider the constraint on g 2 ≡ g 2 aγγ F 2 f ϕ 2 DM /4, the square of the effective amplitude of polarization oscillation, then we convert it to a constraint on g aγγ .The probability density for ĝ2 ≡ ĝ2 aγγ F 2 f ϕ 2 DM /4 given g 2 true , P (ĝ 2 |g 2 true ), can be generated by Monte Carlo simulations using Eqs.(32) and (33).To quantify the preference of the null-hypothesis over the presence of signal, we consider the ratio between the probability of null-hypothesis P null and the alternative hypothesis P alt .The probability P null is obtained as where the integration range is determined to satisfy following three conditions: The probability P alt is obtained in similar way, P alt = P null (g 2 best (ĝ 2 obs )).The g 2 best (ĝ 2 obs ) is g 2 true which maximizes P null within the physical region.Finally, we obtain the 90% confidence interval of g 2 as CI = g 2 true P null (g 2 true ) P alt ≥ 0.9 (C5) Since we do not find evidence for a signal with significance level larger than 3σ, we only report a 95% upper bound from the upper side of the 90% confidence interval.The upper bound of g aγγ F f ϕ DM /2 is square root of that of g 2 .From the relation between the mean amplitude of ALP (= ϕ DM ) and the dark matter density (Eq.( 4)), the 95% upper bound on g aγγ is where κ is the fraction of the dark matter that ALP constitutes, and ρ 0 is the local dark matter density.

FIG. 1 . 2 DM / 2 .
FIG. 1. Top: Comparison of typical scales of ALP field and the typical scales of Tau A. Bottom: Comparison of the typical time scales of ALP field and the observation time of this study.The gray shaded region represents the ALP mass scale of interest in this study.

FIG. 2 .
FIG. 2. Polarization map of Tau A observed by Polarbear in 4th and 5th season from 2015 to 2016.The color scale represents the polarization intensity of Tau A, Q 2TauA + U 2 TauA .The orientations of white bars in map pixels represent the average polarization angles (ψ0) at each map pixel.The red dashed circle represents the integration region to estimate average Tau A polarization angle or to evaluate systematic errors, which is consistently used in this study.

FIG. 3 .
FIG. 3. Comparison of the PWV dependence of the two intensity to polarization coefficient estimates λ4 for Jupiter observations.Each point is a median of all detectors for one observation.The PWV is scaled by sin(el) to compare the line of sight integrated value.

FIG. 4 .
FIG. 4. Comparison of the each detector's intensity to polarization coefficient λ4 estimated by the Jupiter signal and the atmospheric 1/f signal in Tau A observations.

FIG. 5 .
FIG. 5. Overall transfer function and individual transfer functions over the frequency bins (Eq.(20)).The transfer function is averaged over all the phases of oscillations.Bottom panel shows the transfer function F f , the top panel shows 1 − F f .The calibration of relative polarization angle between detectors using Tau A has the effect of reducing the signal almost uniformly to 99.5% (polcal).The absolute polarization angle (ψ0) calibration has the effect of reducing the signal at certain frequencies that our observing schedule does not favor (abscal).The averaging over one-hour observation time has the effect of reducing the signal at high frequencies (average).

FIG. 7 .
FIG. 7. Correlation coefficient between Tau A angle and systematic templates.The vertical lines represent correlation coefficient between measured Tau A angle and systematic templates.The histograms represent correlation coefficient between the noise simulations and systematic templates.

FIG. 8 .FIG. 9 . 1 . ( 48 )
FIG. 8. Top:The measured polarization angle of Tau A superimposed on the best-fit sine curves.The phase of the oscillation (θ) is defined in MJD.Number of observations is small from January 2016 to July 2016, because angular distance between Sun and Tau A is small in this period, also because of the maintenance of the instruments.There are two equally significant correlated oscillation modes at 1/61 day −1 and 1/52 day −1 .Bottom: Amplitudes of the LSSA of data compared with the approximated amplitudes with certain local and global significance levels.The black dots shows the LSSA of data at the frequency bins of Eq. (19), the black line shows the LSSA of data at 10 times more finely sampled frequency bins.The orange solid line shows the LSSA of data minus the best fit sine curve at 1/52 day −1 .

FIG. 12 .
FIG. 12. 1σ upper bound of the statistical uncertainty compared with the systematics.The overall systematics are formed by quadrature sum of individual systematics.The right panel shows the correlation coefficients between different systematics which are smaller than 0.3.
[69] simulated this effects for the exact same telescope with the Polarbear-2 receiver.The leading effect of the cross-polarization, the dipole component of Q and U mixing, is considered to put the upper bound on systematics.The dipole component of the detectors cancel one another due to their focal plane distribution, but the pattern of cancellation differs from day to day, resulting in time-dependent systematic polarization angle fluctuation.We simulated the systematic polarization angle fluctuation due to the residual dipole component assuming the worst orientation of the dipole component and 2% of Q and U mixing.The upper bound of the systematic polarization angle fluctuation is 0.01 • .

FIG. 13 .
FIG. 13. 2σ upperlimit of the LSSA of the systematic templates.The gray solid line shows the 1σ upper bound of the statistical uncertainty.The vertical red dashed lines show the frequencies where we found a hint of a signal.

FIG. 14 .
FIG. 14. Season by season spectrum of the Tau A polarization angle.The vertical red dashed lines show the frequencies where we found a hint of a signal.The hint of a signal is consistently seen in two seasons.

FIG. 15 .
FIG. 15.The measured polarization angle of Tau A as a function of local time at the observation site.

FIG. 17 .
FIG. 17. Top: Light-curves of Tau A measured in this study and various all-sky survey satellites in arbitrary units [73-77].Bottom: LSSA of the light-curves.The vertical red dashed lines show the frequencies where we found a hint of a signal.No clear peaks are found around these frequency bins.
Fig 18   indicates the significance of the discrepancy.

FIG. 18 .
FIG. 18.Comparison of the two most significant signals found in this study with other bounds.The red and blue line shows the 95% upper bound by other measurements and the black line shows the likelihood of this study.Top panel compares the amplitude of an axion-like polarization oscillation.Bottom panel compares the ALP-photon coupling.
•and Polarbear is at −23 • latitude.Therefore, the highest elevation of Tau A observed by Polarbear is 45

TABLE III .
PTEs of null statistics for individual data-splits NOTE: * indicates observation-split type of the data split.

TABLE IV .
PTEs of representative test statistics used for null test pass criteria.

TABLE V .
Error budget of the polarization angle of Tau A. The median statistical uncertainty per observation and the RMS of the systematic polarization angle fluctuations are shown.