QUBIC: an algorithm for detecting cosmic rays

Introduction: QUBIC (the Q $\&$ U Bolometric Interferometer for Cosmology) is an international ground-based experiment designed to observe the polarization of the cosmic microwave background. It has been installed at a high-altitude site in Alto Chorrillos, Argentina (4,869 m above sea level). At this altitude, the cosmic ray flux is high, thus requiring an advanced algorithm for their detection and removal from raw data. Cosmic rays can leave two types of traces in the data: above and below the noise level. This article describes an algorithm for detecting the above-noise traces.

Methods: An algorithm was developed for detecting cosmic ray events in the time-ordered data (TOD) of transition-edge sensor detectors (TES bolometers). Raw signals were pre-processed to obtain de-noised data. Events were searched by applying a threshold to isolate segments showing a rapid increase and subsequent exponential decay. The final goal is to fit each segment to extract the time scale of the candidate and verify the fit quality statistically.

Results: The cosmic ray detection algorithm was applied to datasets acquired in Salta (Argentina) in 2022, during a testing campaign. So far, no candidates have been found after exploring different thresholds for initiating the cosmic ray search, along with various combinations of minimum points required for the sudden increase and exponential decay expected in the signal.

Discussion: We select only high signal-to-noise regions to find the most energetic cosmic ray candidates matching the filters proposed in the method. The null result is not surprising since, for the energy range of cosmic muons of interest here (approximately $[\mathrm{1,100}]$ Gev), the expected energy deposited in our very thin bolometer membranes is small and produces a small signal with respect to the measured noise. However, this methodology could be applied to future longer campaigns to estimate, from the largest (and rare) cosmic ray energy depositions, the TES time constants.

1 Introduction

QUBIC (the Q $\&$ U Bolometric Interferometer for Cosmology) is an international ground-based experiment designed to observe the cosmic microwave background (CMB) polarization, which is one of the major challenges in observational cosmology. Optimized for the measurement of B-mode polarization, QUBIC is designed with a novel approach that combines the advantages of interferometry, in terms of control of the systematic effects, with those of transition-edge sensors (TESs), in terms of wide-band, background-limited sensitivity. QUBIC features an optical system consisting of back-to-back horns that select the relevant baselines and an optical combiner that focuses radiation onto a bolometric focal plane. The optical combiner forms interference fringes, while the bolometers integrate power over timescales that are much longer than the period of the incoming electromagnetic waves. The whole instrument operates at cryogenic temperatures within a large cryostat. The current version, known as the QUBIC technological demonstrator (QUBIC TD), is composed of an array of 64 back-to-back horns and mirrors reduced according to the illumination of the 64 horns. A single array of 256 TES bolometers, readout by two application-specific integrated circuits (ASICs), operates at 150 GHz (Hamilton et al., 2022).

As QUBIC operates in the millimeter wave region of the electromagnetic spectrum, the instrument needs a clean and dry atmosphere. Consequently, it has been installed at a high-altitude site in Alto Chorrillos, Argentina (4,869 m above sea level). At this altitude, there is a significant impact on the measurements caused by cosmic rays (CRs), mainly produced by high energetic cosmic ray muons interacting with the detectors. We expect that most of the detected CRs originate from interactions of primary CRs at the top of the Earth’s atmosphere. They are secondary particles resulting from air showers triggered by hadronic interactions between primary CRs and atmospheric nuclei. There are two types of traces left in the data by CRs: above and below detector noise. CR events are of particular importance for QUBIC for two complementary reasons. From an observational standpoint, CRs represent a source of fast transients that can significantly contaminate the time-ordered data (TOD). Their sudden increase and exponential decay break the statistical stationarity of the noise, introduce non-Gaussian features, and, if not properly identified and removed, can propagate into the map-making stage, ultimately degrading the sensitivity to CMB polarization. At the same time, CRs provide a valuable diagnostic tool for the instrument: each CR impact deposits an almost instantaneous energy impulse in the TES, exciting the thermal response of the detector. The decay that follows this excitation directly probes the effective time constant of the TES. Since the TES response is described by a single-pole exponential in the time domain (equivalently, by a first-order low-pass filter in the frequency domain), the time constant extracted from CR events offers an independent and natural in-flight measurement of the instrument’s transfer function. Knowing this time constant is crucial for recovering signal linearity: once it is measured, the TOD can, in principle, be deconvolved from the detector transfer function, thus correcting for the low-pass filtering introduced by the TES and enabling a more faithful reconstruction of the sky signal. For these reasons, detecting CR events is essential both to preserve data quality and exploit the diagnostic information they carry about the instrument’s temporal response. In this article, an algorithm is described to detect the leftover traces above noise, i.e., in a high signal-to-noise region.

When an amount of energy $\mathrm{\Delta }E$ is deposited in a bolometer (for example, due to the impact of a CR) in a time much shorter than the bolometer time constant, the temperature increase in the bolometer is $\mathrm{\Delta }T=\mathrm{\Delta }E/C$, where C is the heat capacity of the detector. This peak level is achieved in less than one time constant; the subsequent decay, instead, is regulated by the time constant itself. If the noise of a bolometer corresponds to physical temperature fluctuations smaller than $\mathrm{\Delta }E/C$, most of the CR-induced spikes will be evident in the bolometer TOD; otherwise, and in any case, the presence of small signal spikes due to CRs and hidden in the noise increases the noise variance (Masi et al., 2010). In voltage-biased TES, the effective time constant is not set solely by the intrinsic ratio (Equation 1)

${\tau }_0=\frac{C}{G},(1)$

where $G$ is the thermal conductance to the bath. Biasing in the superconducting–normal transition establishes a strong negative electro-thermal feedback (ETF) with loop gain (Equation 2):

$\mathcal{L}=\alpha \frac{P_{\text{bias}}}{GT},\alpha =\frac{T}{R}\frac{\partial R}{\partial T},(2)$

where $P_{\text{bias}}$ is the Joule power and $\alpha$ is the logarithmic temperature sensitivity of the TES resistance. Under typical operating conditions, the detector behaves as a first-order system, with a shortened effective thermal time constant (Equation 3),

$\tau =\frac{{\tau }_0}{1+\mathcal{L}},(3)$

often an order of magnitude faster than ${\tau }_0$. Consequently, the decay phase of CR-induced spikes is governed by ETF dynamics rather than by the bare thermal capacity and must be modeled according to Irwin and Hilton (2005). The decay constant $\tau$ measured from the exponential relaxation of CR-induced transients corresponds directly to the effective thermal time constant introduced in Equation 3. A CR deposits energy impulsively, driving the TES out of equilibrium on a timescale much shorter than its response time, so the subsequent recovery is entirely determined by the ETF dynamics and bath temperature that define $\tau$. Thus, the decay observed in the TOD is a direct expression of the TES single-pole response modeled using Equation 4:

$f\left(t\right)=e^{-\frac{t}{\tau }}.(4)$

Knowing $\tau$ also provides access to the instrument’s transfer function, which behaves as a low-pass filter in the frequency domain. The readout system adds another low-pass filter, which is expected to be much faster than the detector time constant. This creates a fast rise-time, which is not analyzed here.

Another essential component of the QUBIC readout system is the superconducting quantum interference device (SQUID). The direction of the SQUID coil winding, which couples the SQUID to the TES, affects the current–voltage (IV) response and TES polarity and must be considered during data analysis to correctly interpret the direction of CR-induced transients.

2 Methods

To investigate the impact of CR candidates on TODs, a detection algorithm was developed. The description of the methodology used to perform this analysis would focus on a single TES of the focal plane; however, through multiprocessing, the algorithm developed in Python can perform the analysis of the entire focal plane.

The dataset analyzed in this work was acquired during the QUBIC technical run of 2022 in Salta (1,150 m a.s.l.). It consists of TODs from the 256 TES detectors of the focal plane, sampled at 157 Hz, and includes azimuthal scanning observations performed under stable cryogenic and readout conditions. Prior to analysis, all datasets underwent the standard QUBIC quality-control procedures, including verification of TES stability through IV-curve checks. Only TES passing these controls were included in the final CR analysis.

2.1 Data pre-processing

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.

Figure 1

Figure 1. TOD of TES 33 before (blue) and after (orange) the pre-processing procedure.

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.

2.2 Candidate identification

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.

Figure 2

Figure 2. Detection of cosmic ray candidates without the restrictive filters.

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.

2.3 Signal fitting procedures

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\left(t\right)=A+B{e}^{-\frac{t}{\tau }},(5)$

where $A$ is an offset, $B$ is the amplitude of the event, $t$ is the temporal array, and $\tau$ is the time constant. This process yields the time constant associated with the energy dissipation following the event.

2.4 Statistical validation

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\sigma$ from its average value, the candidate is discarded. The $2\sigma$ 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 $\alpha =0.05$ (where $\alpha$ is the probability of a Type-I error, i.e., falsely rejecting a true exponential model); tighter $1\sigma$ or looser $3\sigma$ 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.

• It enables the use of Pearson’s test, which is not suitable for testing non-linear relationships between points.

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\ge r_t$ and $p\ge 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\left(t\right)=at+b,(6)$

obtained through ordinary least squares, and the RMSE is computed as shown in Equation 7:

$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{\frac{1}{N}{\displaystyle \sum _{k=1}^N}{\left[\log s_k-\left(a{t}_k+b\right)\right]}^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\ge 0$ (i.e., the exponential trend is not decreasing);

• $\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}>{\epsilon }_{\text{RMSE}}$, where the default threshold is ${\epsilon }_{\text{RMSE}}=0.02$;

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):

$\frac{{\sigma }_{\text{next}}-{\sigma }_{\text{previous}}}{{\sigma }_{\text{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):

${\chi }^2=\frac{\sum {\text{res}}^2}{\nu {\sigma }^2},(9)$

where $\nu$ is the number of degrees of freedom.

2.5 Visualization and mapping

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.

3 Results

3.1 Cosmic ray detection

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.

• Histogram of the amplitudes fitted with a Landau distribution, which is an asymmetric probability distribution—most commonly used in high-energy physics to describe the energy loss of charged particles as they pass through a thin layer of material.

• Distribution of residuals between the fitted and observed decays, fitted with a combined normal–Cauchy distribution.

• Two plots to show the relationship between the time scales obtained for the candidates of TESs, which, at the end of the overall analysis, are of the same order of magnitude for that TES, along with $V_{\text{bias}}$ and the temperature of the dataset from which they were extracted;

The reliability of the testing framework was evaluated by varying all major algorithmic parameters, including the high SNR threshold (3, 4, $5\sigma$) 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.

3.2 Expected cosmic ray rates

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)], ${\mathrm{\Phi }}_{\text{Bogotá}}=255.1\pm 5.8,{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{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 $\mathrm{\Delta }h=1507$ m between the two sites, the multiplicative factor is presented in Equation 10:

$C=1.1^{\frac{\mathrm{\Delta }h}{300}}\approx 1.61,(10)$

leading to an estimated flux at the QUBIC site of Equation 11:

${\mathrm{\Phi }}_{\text{QUBIC}}\approx 411{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}.(11)$

The geometric area of a single bolometer membrane is taken to be $A_{\text{TES}}=3\mathrm{m}\mathrm{m}\times 3\mathrm{m}\mathrm{m}$ (Piat et al., 2022), and accounting for the 256 TESs, the effective focal-plane area is $A_{\text{FP}}=2.3\cdot 10^{-3}{\mathrm{m}}^{\mathrm{2}}$. The expected CR rate is therefore presented in Equation 12:

$R_{\text{QUBIC}}={\mathrm{\Phi }}_{\text{QUBIC}}\times A_{\text{FP}}\approx 0.94{\mathrm{s}}^{-\mathrm{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\approx 2.019$, a flux ${\mathrm{\Phi }}_{\text{QUBIC}}\approx 515{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$, and thus a rate of $R_{\mathit{QUBIC}}=1.4s^{-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\mathrm{G}\mathrm{e}\mathrm{V}$ (Borja et al., 2022). In the $500\mathrm{n}\mathrm{m}$-thick silicon nitride membrane (Piat et al., 2022), muons in this energy range deposit on the order of $320\mathrm{e}\mathrm{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:

$\mathrm{\Delta }T=\frac{E_{\text{deposited}}}{C},(13)$

where $E$_deposited is first converted into Joules and $C=10^{-11}\mathrm{J}{\mathrm{K}}^{-\mathrm{1}}$ is the detector heat capacity. Using the electric responsivity $\mathcal{R}=2\times 10^8\mathrm{A}\mathrm{D}\mathrm{U}{\mathrm{K}}^{-\mathrm{1}}$, the expected amplitude of the cosmic ray spikes in ADU is then obtained as Equation 14:

$S_{\text{CR}}\left[\mathrm{A}\mathrm{D}\mathrm{U}\right]=\mathrm{\Delta }T\mathcal{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.

4 Discussion

In this work, we have proposed a method for identifying above-noise CR transients in QUBIC’s TES TODs by combining threshold-based event segmentation, two-stage fitting, and statistical validation.

The algorithm has been applied to data from the testing phase in Salta (Ar), and at the current stage of the analysis, no CR candidates have been found. The methodology developed here is directly applicable to future QUBIC longer campaigns to estimate, from the largest (and rare) cosmic ray energy depositions, the TES time constants, and to other bolometric experiments operating in high–cosmic-ray environments. At Salta’s altitude, the CR spectrum is dominated by secondary muons, with an energy distribution extending from 1 to 100 GeV. In a $500\mathrm{n}\mathrm{m}$ layer of silicon nitride (Piat et al., 2022), muons in this energy range deposit on the order of $320\mathrm{e}\mathrm{V}$, which produces an amplitude in ADU of approximately 1026. Events of this amplitude, even occurring at a rate of approximately 3409 per hour, nonetheless, remain buried in the instrumental noise.

The algorithm can, nonetheless, be integrated into offline or near-real-time pipelines with negligible computational overhead, providing a robust framework for automated transient rejection. Moreover, the ability to extract TES time constants from CR hits offers a valuable in-flight calibration opportunity. The modularity of the approach, combined with its statistically driven validation, makes it suitable for deployment in future instruments aiming to characterize transient energy depositions without compromising science data or map-making performance.

Data availability statement

The data analyzed in this study is subject to the following licenses/restrictions: access to the datasets is restricted by the QUBIC collaboration’s internal policy. Requests for access can be evaluated on a case-by-case basis. Requests to access these datasets should be directed toaGFtaWx0b25AQVBDLklOMlAzLkZS.

Author contributions

SF: Data curation, Methodology, Software, Formal analysis, Visualization, Writing – original draft, Writing – review and editing. EBt: Writing – review and editing. SM: Writing – review and editing. EBr: Writing – review and editing. AC: Writing – review and editing. BC: Writing – review and editing. PD: Writing – review and editing. GD: Writing – review and editing. MG: Writing – review and editing. J-CH: Writing – review and editing. NM: Writing – review and editing. GI: Writing – review and editing. CO’S: Writing – review and editing. AP: Writing – review and editing. MP: Writing – review and editing. FS: Writing – review and editing. CS: Writing – review and editing. ST: Writing – review and editing. MZ: Writing – review and editing.

Funding

The author(s) declared that financial support was received for this work and/or its publication. QUBIC is funded by the following agencies. France: ANR (Agence Nationale de la Recherche) 2012 and 2014, DIM-ACAV (Domaine d’Intérêt Majeur-Astronomie et Conditions d’Apparition de la Vie), CNRS/IN2P3 (Centre national de la recherche scientifique/Institut national de physique nucléaire et de physique des particules), CNRS/INSU (Centre national de la recherche scientifique/Institut national et al. de sciences de l’univers). Italy: CNR/PNRA (Consiglio Nazionale delle Ricerche/Programma Nazionale Ricerche in Antartide) until 2016, INFN (Istituto Nazionale di Fisica Nucleare) since 2017. Argentina: MINCyT (Ministerio de Ciencia, Tecnología e Innovación), CNEA (Comisión Nacional de Energía Atómica), CONICET (Consejo Nacional de Investigaciones Cieníficas y Técnicas).

Acknowledgements

The computational work reported in this article was performed on the supercomputers at CINECA, Italy, under the 2024-2025 INFN-CINECA agreement and on the CC-IN2P3 computing cluster operated by CNRS/IN2P3, France. The authors gratefully acknowledge the computing time and technical support provided by both facilities.

Conflict of interest

The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Borja, C., Ávila, C., Roque, G., and Sánchez, M. (2022). Atmospheric muon flux measurement near earth’s equatorial line. Instruments 6, 78. doi:10.3390/instruments6040078

CrossRef Full Text | Google Scholar

Gadey, H., Chatzidakis, S., and Farsoni, A. T. (2020). Monte carlo characterization of the cosmic ray muon flux in shallow subsurface geological repositories intended for disposal of radioactive materials. Appl. Radiat. Isotopes 163, 109209. doi:10.1016/j.apradiso.2020.109209

PubMed Abstract | CrossRef Full Text | Google Scholar

Hamilton, J.-C., Mousset, L., Battistelli, E. S., de Bernardis, P., Bigot-Sazy, M.-A., Chanial, P., et al. (2022). Qubic I: overview and science program. J. Cosmol. Astropart. Phys. 2022, 1–34. doi:10.1088/1475-7516/2022/04/034

CrossRef Full Text | Google Scholar

Irwin, K. D., and Hilton, G. C. (2005). Transition-edge sensors. Top. Appl. Phys. Cryog. Part. Detect. 99, 63–150. doi:10.1007/10933596_3

CrossRef Full Text | Google Scholar

Masi, S., Battistelli, E. S., de Bernardis, P., Lamagna, L., Nati, F., Nati, L., et al. (2010). On the effect of cosmic rays in bolometric cosmic microwave background measurements from the stratosphere. Astronomy & Astrophysics 519, 9. doi:10.1051/0004-6361/201014065

CrossRef Full Text | Google Scholar

Piat, M., Stankowiak, G., Battistelli, E. S., de Bernardis, P., D’Alessandro, G., Petris, M. D., et al. (2022). Qubic iv: performance of tes bolometers and readout electronics. J. Cosmol. Astropart. Phys. 2022, 037. doi:10.1088/1475-7516/2022/04/037

CrossRef Full Text | Google Scholar

Workman, R., Isildak, B., Dogru, A., Aydogan, R., Bayrak, B., and Ertekin, S. (2024). Passage of particles through matter. Prog. Theor. Exp. Phys. 2024, 083C01. doi:10.1093/ptep/ptae068

CrossRef Full Text | Google Scholar