Before searching for cosmic rays, the raw TODs are pre-processed to remove low-frequency trends and optimize the sensitivity to CR impulses. We apply this pre-processing procedure exclusively during azimuthal scans (in the case of fixed pointing, a simple median is removed from the signal), which involve repeated forward and backward telescope movements in azimuth at a fixed elevation to map specific regions of the sky. The objective is to obtain clean, synchronized sweeps from the raw azimuth signal and prepare these data for subsequent analysis. The pre-processing consists of the following steps:

1. Sweep segmentation:

Raw azimuth data are segmented into forward and backward azimuth scans by analyzing the angular velocity, computed as the time derivative of azimuth. A median filter is applied to reduce noise and irregularities, and sweeps are identified based on threshold crossings of the angular velocity.

2. Sweep labeling and indexing:

Each identified forward or backward azimuth scan is labeled and assigned a unique index, enabling clear identification of the data points corresponding to specific sweeps and directions.

3. Time synchronization and interpolation:

Segmented azimuth data are mapped onto the electronics’ time base through interpolation, ensuring alignment between the azimuth data and the electronic system’s sampling times.

4. Baseline correction:

For each forward and backward azimuth scan, the median value of the signal is subtracted, removing any offset between individual sweeps.

5. Trend removal:

The corrected signals are grouped into macro-bins. Trends and slow drifts within each macro-bin are estimated and subtracted through a second-order polynomial fitting, resulting in a stabilized and centered dataset.

6. Final smoothing:

Residual irregularities at sweep boundaries are corrected by linearly interpolating between the start and end points of each sweep, removing any invalid data points in the process. These invalid points may be caused by scan instabilities occurring at the edges of the back-and-forth motion. An example demonstrating the correction and cleaning method applied to the TOD is provided in Figure 1.

The direction of SQUID coil winding is also taken into account. In cases where coils are wound in the opposite direction (e.g., counterclockwise), the IV characteristics—and presumably the spikes due to CRs, are inverted. Accordingly, the TES signal is sign-corrected to maintain a consistent analysis framework for transient direction.

The signal from a CR candidate can be modeled as a vertical linear growth, followed by a single-pole exponential decreasing trend. Segmentation techniques isolate signals corresponding to potential transient events (CR candidates), ensuring that each analyzed segment accurately represents an individual occurrence. We defined the candidate search region using a high signal-to-noise ratio (SNR) threshold, set to three, four, and five times the signal’s standard deviation plus the average value of the signal, to explore the parameter space of the search algorithm. When a signal exceeds the threshold, it triggers the search for the characteristic CR pattern. The algorithm determines whether the data points immediately following the initial trigger point conform to an exponential decay. This check is applied particularly to the first point above the threshold because the method is designed to detect prominent signal bumps, whose decaying nature is immediately identified after the bump.

The algorithm is sensitive to such trends even if their amplitude falls below the initial detection threshold as long as it remains above the baseline defined by the signal’s average plus its standard deviation.

Concurrently, we explored other parameters, including the minimum number of points required to define the vertical linear growth (1, 2, and 3) and the exponential decreasing trend (2, 3, 4, 5, and 6).

To illustrate the behavior of the detection pipeline and the fitting strategy adopted in this work, Figure 2 presents the examples of candidates that have been identified, removing the restrictive filters explained so far. The signal exhibits a small increase, followed by an exponential decay, and was therefore fitted using the model described in Section 2.3, Equation 5.

Although the fit formally converges, these candidates fail the pipeline because they do not exhibit a pronounced bump in the increasing part of the signal, which is characteristic of a high-energy cosmic ray event, and they may represent signal artifacts. As a result, they do not pass the validation steps required by the algorithm; therefore, the events do not correspond to a cosmic ray-induced transient but are instead compatible with instrumental fluctuations or external interference. Events of this type are not selected by the algorithm, and no CR candidates are found in the dataset analyzed in this work.

For each identified candidate, a two-stage fitting procedure is applied. Initially, the increasing edge is characterized using linear regression to estimate the slope. Subsequently, the following exponential decay model is fitted to the decay profile using non-linear least squares optimization: f t = A + B e − t τ , (5) where A is an offset, B is the amplitude of the event, t is the temporal array, and τ is the time constant. This process yields the time constant associated with the energy dissipation following the event.

The reliability of the fitted models is assessed using chi-square tests, p-value estimations, and residual analyses. These statistical metrics quantify the agreement between the observed data and the model, providing a framework for validating each time-scale estimation.

First, to assess whether a candidate truly follows an exponentially decreasing trend, the natural logarithm of the decreasing points is computed. Then, the difference between consecutive log-transformed values is calculated. If this difference deviates by more than 2 σ from its average value, the candidate is discarded. The 2 σ threshold represents the 95 % confidence interval of the nearly Gaussian noise, rejecting points that are statistically incompatible with a single-exponential decay at a significance level α = 0.05 (where α is the probability of a Type-I error, i.e., falsely rejecting a true exponential model); tighter 1 σ or looser 3 σ bounds would eliminate too many decays or retain too many spurious decays, respectively. Using the logarithm of the decreasing trend provides two main benefits:

It emphasizes small variations, allowing for better control over how individual points of a candidate change.

Regarding Pearson’s test, the Pearson correlation coefficient r and the p-value p are also used as criteria to reject candidates. In particular, the thresholds used in the algorithm are r t = − 0.9 and p t = 0.05 , indicating that it is unlikely that a candidate follows a linear decay if it exhibits r ≥ r t and p ≥ p t . Thus, Pearson’s test answers the question: how likely is it that the correlation between the data and a linear decay (obtained by linearizing the candidate) is due to chance rather than reflecting an actual decay pattern?

An additional goodness-of-fit criterion based on the root-mean-square error (RMSE) of the linear regression applied to the log-transformed decay points is evaluated. Let { ( t k , log s k ) } k = 1 N denote the points belonging to the candidate. These points are fitted with a straight line (Equation 6), log ⁡ s t = a t + b , (6)

obtained through ordinary least squares, and the RMSE is computed as shown in Equation 7: R M S E = 1 N ∑ k = 1 N log s k − a t k + b 2 . (7)

Because the data are expressed in logarithmic units, the RMSE quantifies the typical multiplicative deviation from an ideal exponential decay; for example, an RMSE of 0.02 corresponds to an average amplitude error of approximately 2 % .

A candidate is rejected if either of the following conditions is satisfied:

the fitted slope is non-negative, a ≥ 0 (i.e., the exponential trend is not decreasing);

The RMSE test complements Pearson’s correlation: Pearson’s coefficient assesses the linearity of the log-decay, whereas the RMSE quantifies the absolute misfit, discarding candidates that are formally linear yet deviate from a pure exponential by more than the specified tolerance.

To identify candidates that exhibit a significant deviation from the surrounding signal baseline, the average value of each candidate is calculated. If 20 surrounding points (corresponding to approximately 0.12 s at a sampling frequency of 157 Hz) exceed this average, the candidate is rejected. This criterion ensures that only events increasing above the surrounding fluctuations are classified as CR candidates.

Concerning the uncertainties over candidate points, as an initial rough estimation, they are set to the standard deviation of 20 points before and after the candidate (excluding the increasing edge). This can lead to an overestimation of the uncertainties. To obtain a more accurate estimate, an iterative procedure is adopted. In the first step, the standard deviation of the residual for a given candidate is computed. The fitting procedure is then repeated using this new estimate of the standard deviation. The process continues until the relative change (Equation 8): σ next − σ previous σ previous < ϵ , (8) where ϵ is set to 1%, assuming that candidate points show minimal fluctuations once identified. The final value of the standard deviation is used to compute the reduced chi-square (Equation 9): χ 2 = ∑ res 2 ν σ 2 , (9) where ν is the number of degrees of freedom.

Several visualization tools are used to support data interpretation (i.e., histograms and residual plots). The computational pipeline was implemented entirely in Python, leveraging scientific libraries, including NumPy for numerical computation, SciPy for optimization and statistics, Matplotlib for visualization, and the qubicpack library for initial data handling.

The algorithm was applied to datasets acquired under different scanning strategies, including azimuthal scans. The datasets were further categorized by measurement site, but for this analysis, we focus on a single dataset acquired during the testing campaign conducted in Salta (Argentina; 1,150 m a.s.l., 2022).

However, regardless of the site or scanning strategy, the goal for each dataset is to extract the time scale of the candidates for individual TES detectors, each accompanied by an estimate of the statistical uncertainty. These values allow for direct comparisons between detectors. It is important to note that, following the analysis of datasets from a given campaign, the algorithm determines whether the time scales of the candidates associated with a specific TES are of the same order of magnitude. If not, these time scales are excluded from the final overall estimation. This criterion is based on the assumption that if CR candidates are related to the intrinsic time constants of the TESs, then the time scales detected for a given TES should be consistent in order of magnitude. Residual analysis was performed on each candidate event to validate the exponential fit. Finally, an overall analysis is performed across all datasets, including residual histograms, amplitude distributions, elevation and candidate frequency plots, and best-fit models.

If the algorithm detects candidates across all analyzed datasets, an overall analysis of the TESs on the focal plane would be performed:

Estimation of candidate time constant.

The reliability of the testing framework was evaluated by varying all major algorithmic parameters, including the high SNR threshold (3, 4, 5 σ ) and the minimum number of points required for the increasing and decaying phases. Across these configurations, the algorithm exhibited stable behavior: no detections occurred in the high SNR regions.

The input CR flux was derived from the integrated muon flux over the solid angle of the upper Earth’s hemisphere measured at Bogotà [2,657 m a.s.l., Borja et al. (2022)], Φ Bogotá = 255.1 ± 5.8 , m − 2 s − 1 , extrapolated to the Salta site at 1,150 m a.s.l. Following published measurements in Gadey et al. (2020), the flux increases by approximately 10 % per 300 m of altitude. For the altitude difference Δ h = 1507 m between the two sites, the multiplicative factor is presented in Equation 10: C = 1 . 1 Δ h 300 ≈ 1.61 , (10)

leading to an estimated flux at the QUBIC site of Equation 11: Φ QUBIC ≈ 411 m − 2 s − 1 . (11)

The geometric area of a single bolometer membrane is taken to be A TES = 3 m m × 3 m m (Piat et al., 2022), and accounting for the 256 TESs, the effective focal-plane area is A FP = 2.3 ⋅ 1 0 − 3 m 2 . The expected CR rate is therefore presented in Equation 12: R QUBIC = Φ QUBIC × A FP ≈ 0.94 s − 1 , (12)

corresponding to approximately 3409 events per hour on the 256 TESs. If we perform the same calculation at the QUBIC site (Alto Chorrillos 2,657 m a.s.l), we obtain a factor C ≈ 2.019 , a flux Φ QUBIC ≈ 515 m − 2 s − 1 , and thus a rate of R QUBIC = 1.4 s − 1 , which translates into approximately 5,191 events per hour on the 256 TESs.

We expect CR muons with energies in the range 1 − 100 G e V (Borja et al., 2022). In the 500 n m -thick silicon nitride membrane (Piat et al., 2022), muons in this energy range deposit on the order of 320 e V (Workman, 2024).

For each sampled energy, the corresponding cosmic ray amplitude in ADU was computed by converting the deposited energy from eV to Kelvin through Equation 13: Δ T = E deposited C , (13) where E _deposited is first converted into Joules and C = 1 0 − 11   J   K − 1 is the detector heat capacity. Using the electric responsivity R = 2 × 1 0 8   A D U   K − 1 , the expected amplitude of the cosmic ray spikes in ADU is then obtained as Equation 14: S CR A D U = Δ T R = 1026 . (14)

Events of this amplitude, even occurring at a rate of approximately 3409 events per hour, nonetheless, remain buried in the instrumental noise, quantified as the standard deviation of the TOD (10509.6 ADU), as evident in the TOD shown in Figure 1.