Investigating the Cosmic and Solar Drivers of Stratospheric ⁷Be Variability

1. Introduction

Since the pioneering work of Libby [1], who hypothesized the production of radioactive isotopes in the development of atmospheric nuclear showers initiated by galactic cosmic rays, cosmogenic radionuclides have attracted growing scientific interest. Among them, ⁷Be has been extensively investigated since its first detection in rainwater samples [2], due to its characteristic half-life, emission properties, and distinctive distribution across atmospheric layers. As a natural radioactive tracer, it is a key radionuclide for investigating physical phenomena ranging from atmospheric dynamics [3,4,5,6] to climate change [7,8], but also for monitoring the fallout of potential anthropogenic releases into the atmosphere, as in the case of the Chernobyl disaster in 1986. As a result, since 1987, environmental laboratories worldwide have started the systematic measurement of ⁷Be concentrations in air (dry deposition) and rain (wet deposition) as part of their routine environmental monitoring programs [9,10]. Consequently, a vast amount of ⁷Be concentration data has been accumulated over the years, providing valuable resources also for other areas of physics, such as planetary science and space weather research.

Cosmogenic ⁷Be is present throughout the Earth’s atmosphere with varying concentrations, with maxima occurring in both production and accumulation layers. Although stratospheric production (12–50 km) represents the predominant source of ⁷Be detected at ground level due to its 53-day half-life, the isotope is also present at higher altitudes. Spectrometric measurements on the outer surfaces of spacecraft returning from low Earth orbit [11], as well as data from balloon-borne experiments [12,13,14], have indicated unexpectedly high ⁷Be concentrations in the mesosphere (above 50 km) and lower thermosphere (above 85 km), suggesting the existence of additional production mechanisms in these regions. In particular, SEPs (Solar Energetic Particles), whose lower energies compared to Galactic Cosmic Rays (GCRs) cause them to interact predominantly in the upper atmosphere, have been proposed as a possible source of additional ⁷Be production. By generating particle showers at high altitudes, they may lead to localized modulation of ⁷Be concentrations, as hypothesized for high-latitude regions [15].

Regarding the lower atmospheric layers, the dominant source of ⁷Be remains associated to the interaction of higher-energy GCRs, which are able to penetrate into the stratosphere and trigger the spallation reactions responsible for its formation. This well-established mechanism explains the statistical dependence and correlation observed between GCR flux and stratospheric ⁷Be concentrations. The GCR flux is modulated by the Sun, decreasing during periods of high solar activity because of the expansion of the heliosphere, which results from enhanced solar wind production and a stronger, more irregular interplanetary magnetic field. Scattering and diffusion processes of incoming GCRs within the expanded heliosphere lead to a reduction of their flux in the Solar System and on Earth[16,17]. As a consequence of this mechanism, there is a well-established anti-correlation between GCR flux and solar activity, which in turn translates into a direct anti-correlation between atmospheric ⁷Be concentrations and solar activity.

In the environmental and atmospheric physics fields the correlation between the stratospheric ⁷Be and solar activity at different time scale has been investigated by some authors [18,19,20] but only a few have also examined its correlation with the galactic cosmic ray flux on Earth at the same time [21].

This study focuses on the correlation and statistical dependence between the monthly stratospheric ⁷Be concentration, measured at the Casaccia Research Centre since 1987, and proxy variables of cosmic ray flux and solar activity, retrieved from open science databases. The use of long-term monthly time series provides an appropriate temporal framework for detecting relationships with cosmic drivers and characterizing multi-annual trends. Dedicated de-trending techniques have been applied to filter out signals induced by atmospheric processes, which are typically characterized by shorter-term modulations, allowing for a clearer assessment of the cosmic component. Phenomena such as the seasonal intrusion of stratospheric air into the troposphere (annual modulation) [22,23], the quasi-biennial oscillation (QBO) in stratospheric wind direction with a period of approximately 2.2 years [24], teleconnections with the El Niño–Southern Oscillation occurring on timescales of about 7 years [25], and semi-decadal oscillations [26] are all recognized contributors to ⁷Be variability, and they can hide the low-frequency cosmic modulation signal. By applying specific analytical techniques to disentangle these Earth-induced effects, this study aims to investigate the underlying correlations and dependencies evaluating whether the stratospheric ⁷Be concentration can be considered as a proxy, or anti-proxy, for cosmic ray flux and solar activity.

2. Materials and Methods

2.1. Stratospheric ⁷Be Production, Deposition and Measurement in Casaccia Research Center

The ⁷Be isotope is produced in the Earth’s atmosphere through interactions of galactic cosmic rays (GCR) predominantly with nitrogen and oxygen nuclei [27]. The production reactions are driven by spallation processes, mainly induced by protons and neutrons generated within the atmospheric shower initiated by a GCR. Examples of such reactions are given below:

$$p\left(n\right)+{}^{14}\mathrm{N}\to {}^7\mathrm{Be}+X$$

(1)

$$p\left(n\right)+{}^{16}\mathrm{O}\to {}^7\mathrm{Be}+Y$$

(2)

where X and Y denote the residual products of the respective spallation reactions. Considering the timescales of atmospheric diffusion processes (typically on the order of several weeks), the production of ⁷Be can be assumed to occur quasi-simultaneously along the entire extent of the hadronic shower.

According to theoretical models (e.g. [28]), the ⁷Be maximum production is reached in the stratosphere ($10-50\mathrm{km}$), where showers exhibit their maximum development [29]. Formed in ionized states as Be⁺ or Be²⁺, beryllium rapidly attaches to sub-micron aerosol particles, whose high specific surface area enhances electrostatic interaction efficiency [30,31]. Once attached, ⁷Be gradually drifts toward the ground due to gravitational settling and atmospheric dynamics, reaching deposition velocities of few $\mathrm{cm}/\mathrm{s}$ [32]. The separation layer between the stratosphere and the troposphere, the tropopause (approximately at $12\mathrm{km}$ of altitude), acts as an accumulation layer for ⁷Be during its drifting. The strong reduction in vertical mixing, caused by the sharp gradient in temperature and density in this layer, effectively inhibits the downward transfer of air masses containing beryllium [33]. Approximately 71% of ⁷Be is removed from the atmosphere by deposition processes, with about 68% of this fraction due to wet deposition (precipitation) and around 3% due to dry deposition (aerosol settling). The remaining 29% is lost through radioactive decay in the atmosphere [34,35]. Due to its characteristic single gamma emission at $477.6\mathrm{keV}$, ⁷Be is readily detectable on ground via gamma-ray spectrometry both in rainwater sample (wet deposition) and in aerosol samples collected on filters (dry deposition).

In this study, only dry deposition data are available because the monitoring network in place at the Casaccia Research Center is designed for nuclear emergency preparedness, aiming to detect potential releases from the on-site nuclear facilities and the following plume passage. Air sampling is performed using three high-volume air samplers dislocate in the center, each processing approximately 800 m³ of air per day. Dry deposition of ⁷Be is measured monthly by performing gamma spectrometry on 90 filters, corresponding to daily samples collected over a month from the 3 sampling points. This method allows the determination of atmospheric ⁷Be concentrations with sufficient sensitivity.

2.2. Datasets and Data Analysis

The upper panel of Figure 1 shows the monthly time series of atmospheric ⁷Be dry deposition, as measured at the Casaccia Research Center.

The variables linked to the space weather conditions around our planet are available through two Open Science databases. In this work, neutron flux data from the Oulu station in Finland, part of the global neutron monitor network [36], are used as a proxy for GCR flux at Earth’s surface. This detector was chosen as the closest to the beryllium-7 sampling site among the stations of the global network, with data available starting from 1987. Neutrons, together with muons, are among the most penetrating particles produced in hadronic showers in the atmosphere and are capable of reaching the ground. Their flux is directly proportional to the incoming GCR flux impacting the Earth.

The sunspot number was used in this study as a proxy for solar activity and the corresponding data were retrieved from the Open Science database provided by [37]. These datasets are shown in the middle and bottom panel of Figure 1.

Figure 1. Time series of monthly ⁷Be concentrations measured at Casaccia (top), neutron monitor counts (middle), and sunspot numbers (bottom). The ⁷Be time-series exhibits a clear higher-frequency modulation driven by stratosphere-troposphere exchange processes and meteorological dynamics. Neutron monitor counts and sunspot numbers are shown as proxies for cosmic ray intensity and solar activity respectively. Vertical dashed lines indicate the boundaries of solar cycles.

Figure 1. Time series of monthly ⁷Be concentrations measured at Casaccia (top), neutron monitor counts (middle), and sunspot numbers (bottom). The ⁷Be time-series exhibits a clear higher-frequency modulation driven by stratosphere-troposphere exchange processes and meteorological dynamics. Neutron monitor counts and sunspot numbers are shown as proxies for cosmic ray intensity and solar activity respectively. Vertical dashed lines indicate the boundaries of solar cycles.

2.3. De-Trending Methods

The presence of temporal modulations in the ⁷Be time series, induced by atmospheric exchange processes with no counterparts in cosmic drivers, necessitates the application of de-trending procedures to reliably quantify the correlation and the dependency between these variables. Their presence can indeed introduce artificial negative contributions in the covariance structure, potentially obscuring the true physical relationships under investigation.

Short-term variability in atmospheric ⁷Be concentrations primarily reflects dynamic processes occurring in the troposphere and lower stratosphere, including seasonal cycles, quasi-biennial oscillations (QBO), stratosphere-troposphere exchange events, large-scale circulation patterns, and phenomena such as ENSO (El Niño–Southern Oscillation). These terrestrial drivers impart oscillations on timescales ranging from several months to a few years, superimposed upon longer-term variations controlled by solar activity and cosmic ray modulation. To achieve this separation of timescales, three distinct de-trending techniques were implemented and systematically compared in this study, aiming to isolate the long-term variability of atmospheric ⁷Be concentrations from shorter-term meteorological and seasonal perturbations.

The selection of de-trending parameters for each applied method is therefore constrained by the need to suppress these high-frequency, non-cosmogenic fluctuations while preserving the low-frequency signal components associated with solar and heliospheric processes. Accordingly, parameter choices were guided by the typical temporal scales of the known atmospheric oscillations, with the objective of ensuring that the extracted trends accurately reflect variability on timescales consistent with cosmic-ray-related influences.

The first method consists of applying a running average with a fixed time window. This classical approach, widely adopted in geophysical and climate sciences, is used to suppress high-frequency fluctuations and emphasize low-frequency trends [38,39]. In this work, the running average was computed using a centered window, so that each value represents the mean of observations symmetrically distributed around the point in time. This ensures that the smoothed series remains aligned in phase with the original signal, avoiding the temporal shifts that occur when the averaging is performed only over past values. Additionally, near the boundaries of the time series, where a full window cannot be accommodated, the filter adapts by computing the average over the available data points, even if fewer than the nominal window size. This prevents gaps at the start and end of the series and maintains the continuity of the smoothed trend. Overall, this approach preserves the low-frequency component of the series while effectively reducing short-term variability, making it suitable for highlighting temporal patterns associated with atmospheric or solar influences on ⁷Be concentrations.

Secondly, a Seasonal-Trend decomposition based on Loess (STL) was applied [40]. STL is a decomposition technique that separates a time series into trend, seasonal, and remainder components. The method employs locally weighted regression (Loess) to iteratively estimate the trend and seasonal components. Although the seasonal period must be specified in advance, the use of local regression allows the method to accommodate gradual variations in both the amplitude and phase of the seasonal cycle, as well as long-term modulations in the trend component. Finally, a low-frequency Fourier reconstruction was implemented by summing the first n terms of the Fourier series expansion of the ⁷Be time series [41]. The Fourier series expansion preserves the full temporal resolution of the data, avoiding the loss of detail associated with smoothing techniques. In the present work, these de-trending methods were applied exclusively to the ⁷Be time series: the neutron counts and sunspot number series were analyzed in their original form, as they are considered physical drivers of ⁷Be concentration and already represent temporally smoothed indices of cosmic ray flux and solar activity, respectively. This approach avoids the introduction of mathematical artifacts and preserves the physical meaning of the time-series.

2.4. Advanced Measures of Statistical Dependence

The investigation of correlations among variables characterized by (quasi-)sinusoidal periodicity, phase shifts, and inherent temporal lags, such as the detrended time series of atmospheric ⁷Be concentrations, sunspot numbers, and ground-level neutron fluxes, requires careful methodological consideration. These properties suggest that analyses confined to linear and monotonic associations, as measured by the Pearson coefficient (r), may be insufficient to capture the complexity of interactions. Therefore, employing additional statistical metrics capable of detecting complex and nonlinear dependencies is essential for a robust and comprehensive assessment of the relationships among these variables.

To address these challenges, a suite of advanced dependence measures has been employed, including Spearman’s ($\rho$) and Kendall’s ($\tau$) rank correlation coefficients for capturing monotonic yet potentially nonlinear associations; Distance Correlation ($\mathcal{R}$), which can detect any type of dependence irrespective of functional form; Mutual Information (I), quantifying shared information content between variables and the Maximal Information Coefficient (MIC), designed to identify a wide range of linear and nonlinear relationships with sensitivity to various dependency structures. Collectively, these methods enable a more comprehensive and robust assessment of statistical associations, offering critical insights into the complex dynamics underlying environmental and cosmic phenomena. For the following analyses, we assume that the two time series, whose correlation and dependency is to be investigated, can be expressed as vectors of observations:

$$X=(x_1,x_2,\dots ,x_i,\dots ,x_n)\mathrm{and}Y=(y_1,y_2,\dots ,y_j,\dots ,y_n).$$

(3)

where n denotes the number of observations, and $x_i$ and $y_j$ represent the values of the variables X and Y at time indices i and j, respectively. These generic representations will be used consistently throughout the discussion of the different statistical dependence measures.

2.4.1. Correlation Metrics

Several correlation metrics can be employed to assess the relationship between two variables. The simplest and most widely used is Pearson’s correlation coefficient, which measures the strength and direction of a linear association. However, it is limited to detecting linear relationships and may fail to capture more complex or nonlinear dependencies. To detect monotonic relationships that may not be strictly linear, rank-based correlation coefficients are often used. Unlike Pearson’s r, these metrics are sensitive to ordinal associations and do not assume a specific functional form for the relationship between variables. Two widely used estimators, Spearman’s $\rho$ [42] and Kendall’s $\tau$ [43,44], are employed to quantify the degree of monotonic dependence. Details regarding these correlation coefficient are given in the following:

• Pearson’s correlation coefficient r - the metric measures the strength and direction of a linear relationship between two variables. It is defined as:

  $$r={\displaystyle \frac{\mathrm{Cov}(X,Y)}{{\sigma }_X{\sigma }_Y}}$$

  (4)

  where $\mathrm{Cov}(X,Y)$ denotes the covariance between X and Y, and ${\sigma }_X$ and ${\sigma }_Y$ are the standard deviations of X and Y, respectively. The covariance is given by:

  $$\mathrm{Cov}(X,Y)={\displaystyle \frac{1}{n}}\sum _{i=1}^n\left(x_i-\overline{x}\right)\left(y_i-\overline{y}\right),$$

  (5)

  where $\overline{x}$ and $\overline{y}$ are the sample means of X and Y respectively. The coefficient r ranges from $-1$ (perfect negative linear relationship) to $+1$ (perfect positive linear relationship), with $r=0$ indicating no linear association.
• Spearman’s rank correlation coefficient $\rho$ - the metric measures the strength and direction of a monotonic relationship by computing the Pearson correlation between the ranked values of the variables:

  $$\rho ={\displaystyle \frac{\mathrm{Cov}\left(\mathrm{rank}\right(X),\mathrm{rank}(Y\left)\right)}{{\sigma }_{\mathrm{rank}\left(X\right)}{\sigma }_{\mathrm{rank}\left(Y\right)}}}$$

  (6)

  where $\mathrm{Cov}$ denotes the covariance and $\sigma$ the standard deviation of the ranked variables. In this context, ranking transforms the original vector of values into a vector of ranks, where each rank indicates the position of a value in the ordered dataset. The coefficient $\rho$ ranges from $-1$ (perfect decreasing monotonic relationship) to $+1$ (perfect increasing monotonic relationship), with $\rho =0$ indicating no monotonic association. Spearman’s $\rho$ is particularly effective in detecting relationships that are nonlinear but still consistently increasing or decreasing.
• Kendall’s rank correlation coefficient $\tau$ - the metric is based on counting the number of concordant and discordant pairs among all possible pairs of observations:

  $$\tau ={\displaystyle \frac{n_c-n_d}{\frac{1}{2}n(n-1)}}$$

  (7)

  where $n_c$ is the number of concordant pairs and $n_d$ is the number of discordant pairs, with n representing the sample size. Two observations $(x_i,y_i)$ and $(x_j,y_j)$ are considered concordant if the ranks of both variables increase or decrease together (i.e., $x_i>x_j$ and $y_i>y_j$, or $x_i<x_j$ and $y_i<y_j$). They are discordant if one variable increases while the other decreases (e.g., $x_i>x_j$ but $y_i<y_j$). The coefficient $\tau$ also ranges from $-1$ (complete inversion of ranks) to $+1$ (perfect agreement in rank order), with $\tau =0$ corresponding to no association.

While all three measures capture associations between variables, they differ in their sensitivity and interpretation. Pearson’s r assesses the strength and direction of a linear relationship and assumes normally distributed data. In contrast, Spearman’s $\rho$ and Kendall’s $\tau$ are rank-based measures that capture monotonic relationships without requiring linearity. Kendall’s $\tau$ tends to be more robust for small sample sizes and provides a probabilistic interpretation, representing the difference between the probability that pairs of observations are in the same order versus different orders. Spearman’s $\rho$ can be slightly more sensitive to variations in the strength of monotonic associations but may be more influenced by extreme values in the ranks. Together, these correlation coefficients offer complementary insights into both linear and non-linear associations in the data. For all three correlation coefficients, significance was assessed through classical statistical tests: the t-distribution for Pearson’s r [39], and asymptotic normal approximations for Spearman’s $\rho$ and Kendall’s $\tau$ [44,45].

2.4.2. Dependence Metrics

To quantify statistical dependence beyond monotonic or linear structures, three complementary approaches were employed: distance correlation [46,47], mutual information [48], and the maximal information coefficient [49] . These metrics are capable of detecting complex nonlinear or non-functional associations between variables and are suitable for identifying dependencies that traditional correlation measures may overlook. Details of these approaches are provided in the following:

• Distance correlation $\mathcal{R}$ - this metric measures the statistical dependence between two variables or vectors, capturing both linear and nonlinear relationships. Unlike Pearson’s correlation, distance correlation is equal to zero if and only if the variables are statistically independent, making it a true test of independence. Mathematically, distance correlation is defined as:

  $$\mathcal{R}(X,Y)={\displaystyle \frac{\mathcal{V}(X,Y)}{\sqrt{\mathcal{V}(X,X)\mathcal{V}(Y,Y)}}}$$

  (8)

  where $\mathcal{V}(X,Y)$ denotes the distance covariance between X and Y. The distance covariance $\mathcal{V}(X,Y)$ quantifies how much the pairwise distances between observations in X are associated with the pairwise distances in Y. To compute it, the first step is to calculate the pairwise Euclidean distance matrices:

  $$a_{ij}=\parallel x_i-x_j\parallel \mathrm{and}{b}_{ij}=\parallel y_i-y_j\parallel$$

  (9)
  These matrices are then transformed through a centering process, where the row means (${\overline{a}}_{i·}$) and column means (${\overline{a}}_{·j}$) are subtracted to the initial value $a_{ij}$ and the grand mean (${\overline{a}}_{··}$) is added back, yielding the doubly centered matrices:

  $$A_{ij}=a_{ij}-{\overline{a}}_{i·}-{\overline{a}}_{·j}+{\overline{a}}_{··}$$

  (10)

  and similarly for $B_{ij}$. The squared distance covariance is computed as:

  $${\mathcal{V}}^2(X,Y)={\displaystyle \frac{1}{n^2}}\sum _{i=1}^n\sum _{j=1}^n{A}_{ij}{B}_{ij}$$

  (11)
  By normalizing the distance covariance with the distance variances of X and Y, the distance correlation coefficient $\mathcal{R}$, defined in eq. 8, is obtained. The coefficient $\mathcal{R}$ ranges from 0 (independence) to 1 (perfect dependence) and can detect complex relationships that are not necessarily linear or monotonic. This makes distance correlation a valuable tool for exploring dependencies in environmental and geophysical data.
• Mutual information I - the metric quantifies the amount of information shared between two variables and is rooted in information theory. It is defined as:

  $$I(X;Y)=\sum _{x\in X}\sum _{y\in Y}p(x,y)log\left({\displaystyle \frac{p(x,y)}{p\left(x\right)p\left(y\right)}}\right)$$

  (12)

  where $p(x,y)$ is the joint probability distribution and $p\left(x\right)$, $p\left(y\right)$ are the marginal distributions of X and Y, respectively. The index is always non-negative and equals zero if and only if X and Y are independent. Unlike correlation coefficients, mutual information has no fixed upper bound, and its magnitude depends on the entropy of the individual variables.
• Maximal Information Coefficient MIC - the metric allows to capture complex and potentially nonlinear relationships between time series. This method belongs to the family of maximal information-based nonparametric exploration (MINE) statistics and is designed to detect a wide range of association types, regardless of their functional form or strength. Conceptually, MIC measures how well one variable can be used to predict another by exploring how data points distribute over different grid partitions of the scatter plot. The algorithm divides the plane defined by variables X and Y into multiple grids of varying resolutions (up to a maximum resolution determined by the sample size). For each grid, it computes the mutual information—a measure quantifying how much knowing one variable reduces uncertainty about the other. The MIC is then defined as the highest normalized mutual information observed across all tested grids, formally expressed as:

  $$\mathrm{MIC}(X,Y)=\underset{xy\le B\left(n\right)}{max}\left\{{\displaystyle \frac{I^{*}(X,Y;x,y)}{logmin(x,y)}}\right\}$$

  (13)

  where $I^{*}(X,Y;x,y)$ denotes the maximum mutual information achievable with a grid of $x\times y$ cells, and $B\left(n\right)$ limits the maximum grid resolution based on the sample size n. This normalization ensures that MIC scores remain comparable across different grid sizes. A key advantage of MIC is its capacity to identify not only linear and monotonic relationships but also more complex patterns—such as periodic, exponential, or piecewise associations that traditional measures like Pearson or Spearman correlations might fail to detect. MIC values range between 0 and 1, where values close to 1 indicate strong dependence (irrespective of shape), and values near 0 suggest little or no association.

These measures of statistical dependence expand the analytical framework beyond linear or monotonic relationships, allowing the detection of complex, non-linear, and non-functional associations between variables. By identifying dependencies that traditional correlation coefficients might overlook, they offer a more comprehensive characterization of the data structure and are crucial for investigating subtle interactions in environmental and space-weather studies.

3. Results

3.1. Correlation on the Raw Data

The correlation and statistical dependence between ⁷Be, neutron flux, and sunspot number were initially evaluated using raw, non-detrended data. This preliminary analysis was conducted to establish a baseline reference for comparison with the results obtained after applying various de-trending techniques. The left and middle panels of Figure 2 show that the association between ⁷Be concentration and both neutron flux and sunspot number is modest across all metrics. In contrast, the right panel reveals a strong negative correlation between neutron flux and sunspot number, indicating that these series are closely anti-correlated and reflect the same solar modulation signal. This strong interdependence supports the use of both neutron counts and sunspot number as external drivers of ⁷Be variability, even without de-trending.

Metrics index and related p-values are summarized in Table 1.

The analysis of raw data indicates that the relationship between atmospheric ⁷Be concentration and both solar and cosmic proxies is relatively weak and method-dependent, in contrast to the strong and consistent anticorrelation observed between neutron counts and sunspot number. This suggests that ⁷Be does not respond to solar and cosmic drivers as directly as neutron monitors, but rather integrates additional atmospheric processes that may obscure or modulate the primary signal. These findings highlight the necessity of applying appropriate de-trending techniques in order to isolate the intrinsic, large-scale modulation of ⁷Be by solar and cosmic factors, which may otherwise be masked by dominant seasonal patterns and atmospheric variability.

3.2. De-Trending Methods Results

The main parameters of each de-trending method were selected based on their effectiveness in isolating low-frequency components while minimizing shorter-term atmospheric variability. The specific configurations adopted are summarized below:

• Running Average: A window length of 70 months was chosen to effectively suppress high-frequency atmospheric oscillations, including the annual cycle ( 12 months), the quasi-biennial oscillation ( 2.2 years), and other meteorological variations such as ENSO ( 7 years), while preserving lower-frequency components linked to solar activity, particularly the 11-year solar cycle. This window was identified as an optimal compromise, ensuring that the extracted trend predominantly reflects variability of cosmic origin with minimal contamination from terrestrial atmospheric dynamics.
• Fourier Series: The trend component of the ⁷Be time series was reconstructed using the first five Fourier harmonics, aiming to isolate variability on multi-annual to decadal timescales. This spectral resolution was selected to capture the dominant long-term modulation associated with the 11-year solar cycle, while effectively filtering out higher-frequency fluctuations driven by meteorological and seasonal processes. The resulting low-frequency component defines the primary large-scale variability inherent in the ⁷Be record.
• STL Decomposition: The STL method was applied with a seasonal period of 36 months to suppress intra-annual and inter-annual fluctuations driven by meteorological processes and stratosphere-troposphere dynamics. This configuration isolates the low-frequency trend component of the ⁷Be time series, preserving variability potentially linked to solar-cycle modulation, while reducing the influence of shorter-term atmospheric oscillations.

The results obtained are shown in Figure 3.

The consistency of the results obtained from the three de-trending methods was assessed by calculating the periods of the resulting carrier waves. To estimate systematic uncertainties, the key parameters of each method were varied within a plausible range: ±12 months for the running average window, ±2 harmonics for the Fourier reconstruction, and ±12 months for the STL seasonal period. Results are shown in Table 2.

Considering the associated uncertainties, the periods obtained with the three de-trending methods are consistent with each other and compatible with the typical duration of the solar cycle, which drives the main modulation on long timescales.

3.3. Correlation and Dependence Metrics Results

The correlation and dependence among the investigated time series were assessed using the metrics described in sub-section 2.4. The results of all dependence measures and the associated p-values computed for the different variable pairs and de-trending methods are summarized in Table 3. The pairwise analysis includes also the “Neutron vs Sunspot” relationship as a benchmark, given its well-known anti-correlation, allowing all metrics to be compared against an established physical connection.

Figure 4 reports the numerical results obtained after applying the three de-trending approaches to the ⁷Be concentration data, namely running average, Fourier series reconstruction, and STL decomposition. The plots allow comparison of the statistical relationships between ⁷Be and both neutron counts and sunspot number under different assumptions of trend extraction.

4. Discussion

The analysis of the raw time-series data (see Table 1), performed to establish reference values in the present study, reveals only weak statistical dependencies and low correlation coefficients between the atmospheric ⁷Be concentration and both neutron counts, a proxy variable for the cosmic ray flux reaching Earth, and the sunspot number, a proxy for solar activity. Almost all the dependency and correlation indices show modest values, although with a high level of statistical significance. Conversely, the same analysis performed on the neutron counts and sunspot number time series consistently shows the well-known anti-correlation between these two variables, with high absolute values across all correlation and dependency metrics. These findings support the conservative hypothesis adopted in this work to avoid de-trending the neutron count and sunspot number time series, as they already capture almost all the information regarding the two cosmic phenomena of interest. In the ⁷Be monthly concentration time-series, the presence of high-frequency modulations (see Figure 1) due to atmospheric dynamics, substantially hides the low-frequency modulations which carry the information related to the cosmic drivers phenomena. For this reason, three de-trending techniques have been applied to the ⁷Be signal to extract only the low frequency part. An assessment on the obtained low-frequency carrier signal has been performed evaluating the period of the sinusoidal patterns, finding values consistent with the 11 years period of the Solar cycle (see Table 2). The impact of de-trending on the considered metrics is shown in Figure 5.

The analysis of the de-trended ⁷Be series highlights a significant strengthening of the statistical dependence and correlation with both neutron counts and sunspot number, compared to the raw data results. Although all three de-trending techniques adopted (Fourier carrier, running average, and STL decomposition) lead to similar outcomes, we focus the discussion on the results obtained through the running average method, as this approach systematically yields the highest correlation and dependence coefficients. Focusing first on the relationship between ⁷Be and neutron counts (GCR proxy), dependence metrics after de-trending reach 0.69 for Distance Correlation (+48%), 0.93 for Mutual Information (+202%), and 0.60 for MIC (+140%), highlighting a strong shared low-frequency variability linked to cosmic ray modulation. Correlation indices likewise increase, attaining 0.65 for Pearson (+75%), 0.75 for Spearman (+51%), and 0.55 for Kendall (+61%), indicating a pronounced linear and monotonic association once high-frequency atmospheric fluctuations are removed. Overall, the results indicate that the low-frequency component of the ⁷Be time series exhibits substantial statistical association with neutron counts, consistent with shared variability driven by cosmic ray modulation. The increase observed in dependence metrics suggests the presence of shared information potentially extending beyond linear associations, while higher Spearman’s and Kendall’s coefficients point to strengthened monotonic relationships. Furthermore, the rise in Pearson’s correlation implies that a significant part of this dependency may be described through linear models, although statistical evidence leaves open the possibility of additional non-linear contributions, underscoring the mixed nature, linear and possibly non-linear, of the connection between ⁷Be and cosmic ray flux.

For the relationship between ⁷Be and sunspot number (solar activity proxy), dependence metrics after de-trending reach 0.51 for Distance Correlation (+56%), 0.68 for Mutual Information (+224%), and 0.39 for MIC (+137%), indicating moderate shared variability mediated by cosmic ray modulation. Correlation indices show also an absolute value increase, attaining -0.48 for Pearson (+54%), -0.58 for Spearman (+49%), and -0.40 for Kendall (+51%), consistent with a strengthened inverse monotonic relationship, though less pronounced than for neutron counts. These results suggest that, while the connection between ⁷Be and solar activity remains weaker than for neutron counts, the low-frequency component of ⁷Be still reflects a measurable imprint of solar modulation. This is consistent with the physical mechanism whereby solar activity does not directly regulate the stratospheric ⁷Be production, but rather modulates the intensity of galactic cosmic rays, which constitute the primary source of ⁷Be through spallation processes in stratosphere. The increased dependence metrics indicate potential non-linear associations, while enhanced correlation indices point to a clearer monotonic and partially linear relationship.

5. Conclusions

In this study, we investigated the correlation between ⁷Be concentrations, measured at the ENEA Casaccia site from 1987 to 2024, and proxies of cosmic ray flux and solar activity. The application of de-trending techniques to the atmospheric ⁷Be time series significantly improves the detection of its statistical associations with cosmic drivers such as GCR flux and solar activity. By removing high-frequency atmospheric variability, the analysis reveals a low-frequency signal in ⁷Be concentrations that shows stronger dependence and correlation metrics, particularly with neutron counts (GCR proxy), reflecting the shared influence of galactic cosmic ray flux. Although the associations with solar activity are comparatively weaker, they remain significant and align with the indirect role of the Sun in modulating cosmic ray intensities, which are the direct drivers for the ⁷Be production through spallation processes in the stratosphere.

The study employed a comprehensive statistical framework, combining classical correlation coefficients (Pearson, Spearman, Kendall) with dependence metrics such as Distance Correlation, Mutual Information, and the Maximal Information Coefficient (MIC), enabling the detection of both linear and non-linear associations as well as shared information content.

These results highlight the mixed linear and potentially non-linear nature of the connections between ⁷Be and cosmic phenomena, underscoring the importance of employing advanced signal processing techniques on the ⁷Be time series. From a practical standpoint, this work supports the potential use of atmospheric ⁷Be as a valuable proxy for monitoring long-term variations in galactic cosmic ray flux and, to a lesser extent, for inferring aspects of solar activity. However, careful consideration must be given to the choice of de-trending method, because the ⁷Be concentration is subjected to complex atmospheric transport and deposition processes.

Finally, given the extensive global network of atmospheric ⁷Be monitoring stations and the large datasets accumulated over decades, ⁷Be records represent a promising resource for reconstructing past trends in galactic cosmic ray flux on Earth and may potentially contain valuable information related to solar activity variations.