Revisited Reference Solar Proton Event of 23 February 1956: Assessment of the Cosmogenic‐Isotope Method Sensitivity to Extreme Solar Events

Our direct knowledge of solar eruptive events is limited to several decades and does not include extreme events, which can only be studied by the indirect proxy method over millennia, or by a large number of Sun‐like stars. There is a gap, spanning 1–2 orders of magnitude, in the strength of events between directly observed and reconstructed ones. Here, we study the proxy method sensitivity to identify extreme solar particle events (SPEs). First, the strongest directly observed SPE (23 February 1956), used as a reference for proxy‐based reconstructions, was revisited using the newly developed method. Next, sensitivity of the cosmogenic‐isotope method to detect a reference SPE was assessed against the precision and number of individual isotopic records, showing that it is too weak by a factor ≈30 to be reliably identified in a single record. Uncertainties of 10Be and 14C data are shown to be dominated by local/regional patterns and measurement errors, respectively. By combining several proxy records, a SPE 4–5 times stronger than the reference one can be potentially detected, increasing the present‐day sensitivity by an order of magnitude. This will allow filling the observational gap in SPE strength distribution, thus enriching statistics of extreme events from 3–4 presently known ones to several tens. This will provide a solid basis for research in the field of extreme events, both for fundamental science, namely solar and stellar physics, and practical applications, such as the risk assessments of severe space‐based hazards for modern technological society.


Introduction
Sun is an active star, which sporadically produces eruptive events, such as flares, coronal mass ejections (CMEs), interplanetary shocks, etc., in a wide range of energy release. Most powerful events, particularly related to solar particle events (SPEs) with a huge enhancement of near-Earth radiation environment, can be hazardous for the modern technological society, especially for high-and polar low-orbit spacecraft and space missions (Miyake, Usoskin, & Poluianov, 2019). Such SPEs have been monitored since the 1940s by ground-based detectors and later by space-borne ones. So far, 72 SPEs have been recorded by ground-based detectors (called ground-level enhancement, GLE -https://gle.oulu.fi), ranging from barely detectable to very strong ones. The strongest SPE recorded so far took place on 23-Feb-1956 (GLE #5) and had a huge enhancement of the radiation environment (50-fold over the galactic cosmic-ray (GCR) background, as recorded by Leeds neutron monitor) and a very hard energy spectrum of solar energetic particles (Asvestari et al., 2017). Quite often the strength of a SPE is quantified in terms of the event-integrated (over the entire duration ∆t of the event) omnidirectional fluence of particles with energy above a given threshold E * : where J(E, t) is the differential energy spectrum of solar energetic particles at time t. For production of cosmogenic isotopes by solar energetic particles in the Earth's atmosphere, F 200 is used (G. A. Kovaltsov, Usoskin, Cliver, Dietrich, & Tylka, 2014), which is the event-integrated fluence with energy above 200 MeV (E * = 200 MeV).
Statistics of the recorded event occurrence during the last decades (Raukunen et al., 2018) are summarized in Figure 1 as open triangles. The plot shows the integral occurrence probability density function (IOPDF) of SPEs, which is the probability of an event with the F 200 fluence greater than the given value to occur within a year on average (irrespective of the solar cycle phase). For example, the event of 23- Feb-1956 or stronger (the rightmost open triangle) has the probability to occur ranging between once per 13 and 170 years (see details in G. A. Kovaltsov et al., 2014;. The IOPDF appears quite flat, and its is hardly possible to say whether even stronger (extreme) events may occur and how often. Straight- Integral occurrence probability density function of the annual SPE F200 fluence, according to the space era data (triangles -data from Raukunen et al. (2018) and this work) and cosmogenic proxy data (circles, data from Miyake et al. (2019)). The plot is modified after G. A. Kovaltsov et al. (2014). Open and filled symbols correspond to the known SPEs and conservative upper limits, respectively. Error bars bound the 68% confidence interval. The greyhatched area displays an estimate based on the lunar rock data (S. Poluianov et al., 2018). The range of sensitivity of 14 C-and 10 Be-based reconstructions (see Section 3) is shown by blue and yellow-red gradient filled boxes.
forward extrapolations of this distribution lead to uncertainties as large as several orders of magnitude, because of the insufficient statistics (Lingenfelter & Hudson, 1980). Thus, until recently, estimates of the extreme SPEs parameters and occurrence probability remained highly unreliable (Gopalswamy, 2018).
The situation has changed in 2012, when an abnormal peak in carbon-14 14 C was discovered in data samples corresponding to the year 775 AD (Miyake, Nagaya, Masuda, & Nakamura, 2012;. The peak was so strong that its solar nature was initially rejected and exotic astronomical phenomena were proposed (Miyake et al., 2012). However, it was soon confirmed in other cosmogenic isotope records (beryllium-10 10 Be and chlorine-36 36 Cl) and proven to be caused by a giant SPE (Mekhaldi et al., 2015;Sukhodolov et al., 2017;Usoskin et al., 2013), which was a factor 50 -100 stronger than the strongest directly observed SPE of 23-Feb-1956. Although the 14 C signal peaked in 775 AD, the event occurred around summer 774 AD (Büntgen et al., 2018;Güttler et al., 2015;Sukhodolov et al., 2017;Uusitalo et al., 2018). It appeared to be a black swan -something which is not expected to exist but it does. It was argued that the 774 AD event was the strongest for the Sun over the last twelve millennia Holocene (Usoskin, 2017), as depicted by the rightmost open circle in Figure 1. This sets up the class of extreme solar events. Subsequently two other, slightly weaker (60 -80 % of the one on 774 AD) but still extreme SPEs were found in cosmogenic isotope data, viz., 993 AD (Miyake, Masuda, & Nakamura, 2013) and 660 BC (O'Hare et al., 2019;Park, Southon, Fahrni, Creasman, & Mewaldt, 2017). One more event candidate was found on 3372 BC (Wang, Yu, Zou, Dai, & Cheng, 2017) that still awaits to be confirmed. Thus, at present we know 3 -4 extreme SPEs (each a factor of 30 -100 stronger than that of 23-Feb-1956) to occur during the last several millennia, as reflected by the circles in Figure 1.
The statistics based on cosmogenic data imply a roll-off of the occurrence probability of extreme SPEs. However, there is a large gap of a factor of ≈30 between the space-era and the proxy-based datasets. The gap is caused by observational techniques: on one hand, no extreme SPEs have occurred during the space era, and on the other hand, the sensitivity of the cosmogenic isotope method is presently insufficient to detect 'regular' events, but only extreme SPEs. This leaves a question open of whether the extreme SPEs form a special class of events with different distribution or it is a tail of the same distribution as all SPEs. This question is crucial to study eruptive events and their terrestrial/planetary impacts, and can be answered only by detection of intermediate events, which are a factor 3 -10 greater than that of 23-Feb-1956. However, at present, the sensitivity of the cosmogenic isotope method to SPEs has not been studied in detail, the observational threshold is not defined, and the recommendations on improvements are not set. Some earlier estimates (McCracken & Beer, 2015;Usoskin, Solanki, Kovaltsov, Beer, & Kromer, 2006) were made in an oversimplified way, based on outdated cosmogenic isotope yield functions and neglected details of the atmospheric transport and deposition, which are crucial for reliable detection of weak signals.
Here we present the first systematic assessment of the sensitivity of the terrestrial system to record extreme SPEs based on up-to-date realistic isotope production and transport models and newly available high-resolution data. First, we fix the parameters of the reference SPE of 23-Feb-1956, revisited using the most up-to-date models described in Section 2. Then, we discuss the sensitivity of the cosmogenic isotope method to SPE detection in Section 3, including detailed modelling of the isotope production and atmospheric transport/deposition, for both 10 Be in ice cores and 14 C in tree rings. Finally, we set up the sensitivity of the cosmogenic-isotope method and provide recommendations on filling the observational gap, as well as discuss further prospectives.
2 Reference event of 23-Feb-1956 (GLE #5) revisited Extreme SPEs were recorded in terrestrial cosmogenic isotope archives, which provide energy-and time-integrated responses rather than direct measurements of the particles spectra. However, for modelling of the events and their impacts, energy spectra need to be determined or assumed. It has been shown, using an analysis of strong space-era SPEs, that an event's strength correlates with the spectrum hardness, viz. the stronger SPE is, the harder its spectrum tends to be (Asvestari et al., 2017). This is related either to the efficiency of particle acceleration in extreme conditions, which require simultaneous appearance of several favorable factors, or to interplanetary transport effects, such as the so-called "streaming limit" (Reames & Ng, 2010). In particular, the spectrum of the strongest directly observed SPE of 23-Feb-1956 was one of the hardest for the directly recorded events (e.g., Tuohino, Ibragimov, Usoskin, & Mishev, 2018;Vashenyuk, Balabin, & Miroshnichenko, 2008). Estimates of the extreme SPEs, made using a multi-isotope approach, imply that their spectra were also hard and consistent with that of 23-Feb-1956 in the energy range below a few hundred MeV, and scaled up by different factors (Mekhaldi et al., 2015;O'Hare et al., 2019). Accordingly, the SPE of 23-Feb-1956 is typically used as a reference SPE so that it is simply scaled up to match indirectly found extreme SPEs (Mekhaldi et al., 2015;Sukhodolov et al., 2017;Usoskin et al., 2013). Here we also consider this SPE as the reference one, denoting it further as SPE 1956 . Until recently, the SPE 1956 spectrum was estimated to be in the lower energy range from ionospheric data (Webber, Higbie, & McCracken, 2007) and in higher energy range from neutron monitors, NMs, (Raukunen et al., 2018). However, these earlier estimates of the spectrum, especially in the high-energy range, were based on an obsolete NM yield functions (Clem & Dorman, 2000) and a simplified approach. Here we revise the spectral estimate of the SPE 1956 using a recent verified NM yield function Kovaltsov (2013) updated as Mishev, Koldobskiy, Kovaltsov, Gil, andUsoskin (2020)) and an improved effective-rigidity method (Koldobskiy, Kovaltsov, Mishev, & Usoskin, 2019).
First, we performed a full re-analysis of all the available NM data for the reference SPE, viz. GLE #5. The GLE was observed by 14 NMs (Table 1) located at different geomagnetic cutoff rigidities, from 1 -13.4 GV, as shown in Figure 2. The highest peak count rate (5117% above the background due to GCR) was recorded by Leeds NM, while the greatest event-integrated count-rate increase (see definition in Asvestari et al., 2017) of 5276 %*hr was observed by Ottawa NM. Event-integrated responses of other NM ranged from 14.5 %*hr for the equatorial Huancayo NM to >5000 %*hr for high-latitude stations (see column I in Table 1). Next, we applied the methodology developed by Koldobskiy et al. (2019), so that the effective rigidity R eff and the event-integrated fluence of SEP with rigidity above it F (> R eff ) were evaluated for each NM (see the last two columns of Table 1). These points are shown as black dots with error bars in Figure 3.
The reconstructed spectrum was fitted with a prescribed spectral shape of the event-integrated rigidity fluence of SEP. Here we used a modified Band-function, which includes also exponential roll-off, similar to the Ellison-Ramaty spectral shape   Table 1), green stars represent low-energy fluence from Webber et al. (2007). The blue dashed line is the earlier spectral estimate by Raukunen et al. (2018). The red curve with the shaded area depicts the best-fit functional estimate (J1 = 1.40 · 10 8 cm −2 , γ1 = 1.822, R1 = 0.823 GV, J2 = 1.023 · 10 8 cm −2 , γ2 = 4.207, R2 = 4.692 GV and R b = 2.380 GV -see Equation 2) of the spectrum with the 95% confidence interval, reconstructed here. Table 1. Parameters of the neutron monitors and their responses for the GLE #5 on 23-Feb-1956. Columns are: 1 -name and 2 -standard acronym; 3 -vertical geomagnetic rigidity cutoff for the location and date; 4, 5 and 6 -geographical longitude, latitude, and altitude, respectively of the NM location; 7 -reference barometric pressure; 8 -GLE integral response I; 9 -effective rigidity R eff ; and 10 -the event-integrated fluence of SEP with rigidity above R eff .

Name
Acr  (Ellison & Ramaty, 1985): where F (> R) is the omnidirectional fluence of particles in units of cm −2 , R is the rigidity expressed in gigavolts, parameters γ 1 , γ 2 , R 1 , R 2 and J 2 are defined by fitting, and other parameters are calculated as This function is constructed so that it and its first derivative are continuous. The data points were fitted with the prescribed spectrum using the following Monte-Carlo procedure. First, the exact pair of values of R eff and F (> R eff ) was randomly taken from the values with uncertainties shown in Table 1 independently for each NM. The fit (Eq. 2) was obtained by minimizing the logarithmic discrepancy, using NM data for higher-rigidity range R ≥ R b (viz. parameters J 2 , γ 2 and R 2 ), while the lower-rigidity range was fitted using the data-points from Webber et al. (2007). For the low-energy data points, originally provided by Webber et al. (2007) without error bars, a constant 10% error was artificially applied. The formal value of χ 2 was computed for the fit and saved. Then, a new set of R eff and F (> R eff ) values was randomly obtained, and the fit repeated. Overall, 1000 fits were performed. The one corresponding to the minimum χ 2 was considered as the best fit, but all fits lie very close to the best one. The best-fit spectrum, in the form of Equation 2, is shown in red in Figure 3 along with its 95% confidence interval. The parameters of the best-fit spectrum are given in caption of Figure 3. The 95% confidence-interval uncertainties of the fitting are within 10%.
We used this best-fit spectrum as the reference SPE spectrum in the subsequent computations.

Imprints in cosmogenic records
Cosmogenic isotopes (most useful being 14 C, 10 Be and 36 Cl) are produced as a subproduct of the nucleonic-muon-electromagnetic cascade induced by energetic particles in the terrestrial atmosphere (e.g., Beer, McCracken, & von Steiger, 2012). The bulk of terrestrial cosmogenic isotope is produced by GCR which are always present near Earth, but extreme SPEs with hard spectra also can produce a detectable amount of them (McCracken & Beer, 2015;Usoskin et al., 2006;Webber et al., 2007). This process is well studied and can be modelled via the yield-function of the isotope production (see, e.g., G. Kovaltsov, Mishev, & Usoskin, 2012).
Here we modelled production of the 10 Be and 14 C cosmogenic isotopes due to GCR and the reference SPE 1956 using the set of the yield functions, which provides a full 3D isotope production pattern, as computed recently by S. V. Poluianov, Kovaltsov, Mishev, and Usoskin (2016).

10 Be in polar ice cores
The beryllium 10 Be isotope is unstable with the half-life of about (1.387 ± 0.012) · 10 6 years (e.g., Korschinek et al., 2010). In the atmosphere, it is produced mainly by spallation of nitrogen and oxygen by nucleons of the cosmic-ray induced atmospheric cascade. After production it becomes attached to atmospheric aerosols (e.g., Raisbeck et al., 1981) and relatively quickly (within a few years) precipitates, mostly at mid latitudes but reaches also polar regions. The exact atmospheric path of beryllium is affected by scavenging, stratosphere-troposphere exchange and intertropospheric mixing (e.g., McHargue & Damon, 1991). Deposition of beryllium to ice/snow, where it is subsequently measured, can be dry-or wet-dominated, and its concentration in ice may be additionally affected by post-depositional processes, depending on the snow accumulation rate. Accordingly, the concentration of 10 Be measured in a given ice core does not precisely reflect the global production rate, but may be influenced by the local/regional climate, especially on a short time scale (Usoskin, Horiuchi, Solanki, Kovaltsov, & Bard, 2009). Moreover, strong volcanic eruption may affect beryllium deposition in polar areas (Baroni, Bard, Petit, Magand, & Bourlès, 2011;Baroni, Bard, Petit, & Viseur, 2019).
Accordingly, we performed an analysis of 10 Be series from different ice cores from both hemispheres and model them individually.

Data
We used published data of 10 Be from eight ice cores or snow pits covering the decade around the SPE 1956 -four from Antarctica and four from Greenland, see Table 2. Since the datasets have different resolution, we reduced them all to annual 10 Be concentration, as shown in Figure 4. While tree rings are absolutely dated within the annual accuracy, ice-core dating is less accurate, with uncertainties being relatively small for the last decades, viz. 1 year for high snow-accumulation sites and 1 -2 years (depending to the proximity to tie points such as major volcanic eruptions or nuclear tests, as, e.g., the well-known peak in 1955) for low-accumulation sites (see, e.g., Baroni et al., 2011), but may reach 7 years a millennium ago (Sigl et al., 2015) or up to 70 years for the early Holocene (Adolphi & Muscheler, 2016).

Modelling
We modelled production and deposition of 10 Be for the period 1950 -1960, covering roughly one solar cycle including the maximum phase of the strongest cycle #19 (1954 -1964) and the studied event. Variability of GCR was modelled using the heliospheric modulation potential φ (Usoskin, Alanko-Huotari, Kovaltsov, & Mursula, 2005;Usoskin, Gil, Kovaltsov, Mishev, & Mikhailov, 2017) adopted for the solar forcing incorporated in the CMIP6 (Coupled Model Intercomparison Project, stage 6) project (Matthes et al., 2017). The flux of solar energetic particles during the SPE 1956 was modelled using the fluence spectrum reconstruction described above. All SPErelated 10 Be was assumed to be produced evenly within the day of 23-Feb-1956. 3D spatial distribution (location and altitude) was modelled for each day using the input spectra of GCR and SPE 1956 by applying the 10 Be yield function (S. V. Poluianov et al., 2016).
Transport and deposition of the produced 10 Be were modelled using chemistryclimate model (CCM) SOCOL v.3 (Sukhodolov et al., 2017). The model configuration was the same as applied for the study of 774 AD event (Sukhodolov et al., 2017), but all boundary conditions (sea surface temperature, sea ice concentration, greenhouse gas concentration, land use, geomagnetic field strength etc.) were prescribed for the considered period of time (1940 -1960). The model was spun-up for 10 years starting from 1940 to reach 10 Be distribution similar to 1950. GCR-and SPE-produced 10 Be was traced as different species allowing us to separate their production, transport and deposition.
Depositional flux of 10 Be (in units of atoms cm −2 sec −1 ) was recorded with daily resolution for the sites listed in Table 2, as shown in Figure 5 separately for the GCRand SPE 1956 -related species (blue and red curves, respectively). One can see that the simulated GCR-related depositional signal is very noisy because of the regional/local weather patterns (curves in Figure 5 are smoothed with a 31-day running mean filter), but several features are clearly observed: the seasonal cycle, mostly related to the period of intrusion of stratospheric air (where a major fraction of 10 Be is produced) into the troposphere; and the (inverted) solar cycle, so that almost all the blue curves have a maximum around 1954 -1955 corresponding to the minimum of solar cycle #19 and a minimum around 1958 (maximum of the cycle #19). These cycles are more visible in the globally averaged 10 Be deposition (Figure 5i). It is important to mention that the seasonal variability (roughly a factor of two) seen in the modelled 10 Be deposition flux ( Figure 5) is greater than the 11-cycle one in most of the series, while short sub-annual fluctuations (in particular spikes) may be also of that order of magnitude. The modelled 10 Be deposition related to the SPE 1956 is shown by red curves in Figures 4 and 5, and it would be totally lost among the random short-term variability of the GCR-related signal and the seasonal wave pattern. This fact makes high-resolution 10 Be data not promising in recognition of SPE signals, which can be mimicked by these spikes. A natural smoothing period for the 10 Be data is annual, where the seasonal variability is averaged out. On the other hand, the 10 Be response to the SPE of 774 CE was spread over several years (cf. Sukhodolov et al., 2017) so that about half of the produced 10 Be is deposited within the first year (viz. 1956 here), another quarter -during the second year, and further on, so that each next year half of the previous is deposited, implying that the mean residence time of 10 Be in the atmosphere is roughly one year (e.g., Raisbeck et al., 1981).

Model-vs-data
In order to compare the model results with the measured concentration of 10 Be, the annually-accumulated modelled flux was divided by the simulated annual water (snow) precipitation in the same grid cell. Post-depositional effects were not taken into account. In order to match the exact levels, the simulated values were scaled to the mean measured 10 Be concentration over the period of 1950 -1960. The scaling factors appear to range from 0.8 -1.6 for individual sites, which implies good agreement considering the roughness of the model spatial grid (2.79 • ×2.79 • ) not accounting for the very local meteorology/orography, and the neglect of post-depositional effects. The resultant modelled concentrations are shown in Figure 4 as colored curves (red and blue ones depict results for GCR+SPE and GCR-only related depositional flux) plotted over the measured 10 Be concentrations in different ice cores (see Section 3.1.1). One can see that the SPE 1956 -related signal (the difference between the red and blue curves for the year 1956) is barely visible on the GCR-related background.
Although the agreement between the modelled and measured data is far from perfect, there is a clear tendency of agreement regarding the shape of the 11-year cycle: most of the measured data sets exhibit a decline in the last years of the series. Some series exhibit a 1-year time mismatch between the data and model results, that may be related to the meteorological "noise" or to slightly imprecise dating. The best agreement is observed for the DSS series. It is clear that the SPE 1956 cannot be detected even in the simulated data, as the corresponding 10 Be increase is only ≈5% of the GCR background. The situation is even worse for the actual noisy data with other types of possible errors not existing in the model, such as measurement and dating uncertainties (cf. McCracken & Beer, 2015). Even a 5-fold SPE 1956 (SPEx5) would be not detectable over the noise in a single ice-core series (except for the DSS one) as one can see with the green curve in Figure 4. However, a 10-fold stronger event (SPEx10) would probably be detected even in a single series. On the other hand, the standardized mean of the eight individual 10 Be series (Figure 4i) is much smoother than any individual one, and the SPEx5 signal stands out there at a detectable level. Here we quantify this.
First, for each series we have analyzed the discrepancy between the (normalized) modelled X mod and observed X obs data, as shown in Figure 6a: where X mod is the GCR+SPE model. For most of the series the standard deviation of δ 10Be is 0.2 -0.3, while the DSS series depicts a significantly smaller discrepancy of ≈ 0.1. Figure 6b shows the distribution of the δ 10Be values for all data series together. The distribution can be well fitted with a Gaussian with zero mean and σ = 0.22. Thus, for the subsequent estimates we consider that the uncertainty of an annual 10 Be data point is ≈ 20%. This is significantly greater than the 1σ measurement uncertainty of 10 Be concentration that is estimated as 2 -7% (e.g., Pedro Table 2) and the globally averaged (panels a-h and i, respectively), simulated here using the SOCOL CCM for the period 1950 -1960. Blue curves denote beryllium deposition due to GCR, while red ones -contribution from the SPE1956. All curves are 31-day running mean smoothed. 1997) and is largely dominated by the 'meteorological noise', dating uncertainties or site peculiarity rather than the precision of the measurements.
The modelled response of polar-ice 10 Be to the 2SPE 1956 is ≈4.8% of the mean level (see Figure 5), which is much below the uncertainty of a single annual data point of 20%, making the signal of the reference SPE 1956 in a single 10 Be series at the 0.25σ level. Thus, a 8x reference event would have been detected at a 2σ level, viz. significantly, even in a single 10 Be series, or highly significant, 3.8σ if there were two independent ice cores. The real datasets shown in Figure 4 are poorly correlated with each other so that only 2 out of 28 possible cross-correlations between individual series pairs appear highly significant (confidence level > 0.95), implying that the noise is quite local. The overall sensitivity of the 10 Be data to the reference event is shown in Figure 7a for the scaling factor of the SPE 1956 (X-axis) and the number of independent ice cores (Y-axis). In a realistic case of 4 -5 ice cores, a 4x reference SPE 1956 can be detected at a significant level of 2σ (0.95 confidence). We note that increasing resolution of an isolated ice-core measurement would not enhance the sensitivity because of the seasonal cycle and meteorological noise in the data, but a larger number of far separated ice cores with coarser resolution could help reducing the noise.
For the existing 10 Be opportunities, an optimal detection limit can be set as 4 -5x reference SPE 1956 . This is illustrated in Figure 4i, which shows the standardized mean of all the eight individual series (black histogram with errors bars) along with the model curve for three scenarios: only GCR (blue), GCR+SPE (red), and GCR+SPEx5 (green). Despite a very strong diversity of the individual ice-core data series (panels a -h), the combined time profile is smoother and agrees well with the modelled scenario GCR+SPE: the Pearson's correlation coefficient is r=0.86 (p≈6·10 −4 ), χ 2 (10) ≈ 17. We note that Antarctic sites depict a more consistent 11-year cycle pattern (r = 0.88 between measured and modelled spatially averaged datasets) than Greenland-based series (r = 0.66). The reference SPE cannot be identified in a combined dataset of eight series, since its departure from the GCR-only scenario is only 0.8σ. On the other hand, a 5-fold reference SPE (GCR+SPEx5 scenario) would have led to a 3.7σ enhancement in the year 1956, implying a significant identification of the peak. This is in agreement with the modelled sensitivity shown in Figure 7a, where the expected detection level of a SPEx5 event with 8 different ice-core series is 3.5σ. A SPEx10 scenario (not shown) would have been clearly identified with high statistical significance (7.7σ).

Radiocarbon 14 C in tree rings
Radiocarbon 14 C is produced as a result of (n, p)-reactions (neutron capture) by nitrogen, which makes the main sink of thermal neutrons in the atmosphere. Radiocarbon gets oxidized to 14 CO 2 and takes part in the global carbon cycles. If absorbed by a living organism, e.g., a tree, it remains there decaying in time (half-life is about 5730 years). The relative abundance of 14 C in independently dated samples (e.g., tree rings) provides a measure of cosmic-ray flux at the time of the sample growing (see, e.g., Beer et al., 2012;Usoskin, 2017, and references therein).
Unfortunately, direct 14 C data cannot be used as an index of cosmic-ray intensity for the middle of the 20th century, because of the Suess effect (burning of fossil 14 C-free fuel, which dilutes radiocarbon in the terrestrial system) and nuclear bomb tests that produced anthropogenic 14 C well above the natural level (Beer et al., 2012). On the other hand, the accuracy of modern models is sufficient to reproduce the 14 C signal from known source (e.g., Usoskin, 2017). Figure 8 presents the modelled 14 C time profile for the period under investigation. Computations were performed using the 14 C production model by S. V. Poluianov et al. (2016) and the 22-box carbon cycle model by Büntgen et al. (2018). The input from GCR and the SPE 1956 were treated similar to 10 Be (see Section 3.1.2).
The modelled globally averaged 14 C production rate Q 14C is shown in Figure 8a as reduced to the annual resolution. The (inverted) 11-year solar cycle with the Q−values ranging from 1.3 to 1.9 atom/cm 2 /sec (viz. varying by a factor of ≈1.5) is clear. Annually averaged production of 14 C by the SPE 1956 is ≈0.1 atom/cm 2 /sec (3.04 · 10 6 at/cm 2 in total), viz. 5 -7 % of the GCR-related level, similar to 10 Be.
However, the measurable quantity, which is the relative normalized abundance of 14 C with respect to 12 C, considering isotope fractionizing, ∆ 14 C (e.g., Beer et al., 2012), is not a direct projection of the production signal Q, because of the complicated carbon cycle which acts as a complex signal filter (Bard, Raisbeck, Yiou, & Jouzel, 1997). The 11-year solar cycle is shifted by 3 years (≈ 1 / 4 of the period) and attenuated by a factor of ≈100 to the magnitude of 3 -5 . Since the carbon cycle effectively damps high-frequency signals, the measurable response to the SPE 1956 is attenuated even stronger, by a factor of ≈ 300, down to 0.2 (note that the red curve in Figure 8b is ten-fold enhanced for visibility) and spread over many years. Since the modern AMS (Accelerator Mass Spectrometry) technique allows reaching the measurement uncertainties of 1.5 in a single year measurement (e.g., Güttler et al., 2015;Park et al., 2017), the reference SPE 1956 yields a signal at 0.13σ level, which is undetectable. A 15x enhancement of the reference event is required to reach the level of significant detection in a single 14 C series, which is double compared to 10 Be. The sensitivity of radiocarbon data to the reference event is shown in Figure 7b. One can see that the sensitivity of 14 C to the reference SEP event is nearly half of that for 10 Be. This is mainly due to the fact that, while beryllium, produced by a SPE, precipitates within a few years leading to a recognizable peak (e.g., Sukhodolov et al., 2017), SPE-produced 14 C resides in the atmosphere for many years, spreading the signal over a long period.
On the other hand, since carbon is gaseous and well mixed in the terrestrial system, uncertainties of ∆ 14 C are largely defined by the measurement errors rather than by the regional meteorological noise. Accordingly, almost 'independent' series can be obtained by re-measuring the same samples without requiring a trees from different regions.

Summary and conclusions
Here we present a study of the sensitivity of the cosmogenic isotope method to solar energetic particles, based on the SPE of 23-Feb-1956, which was the strongest directly observed SPE with a very hard spectrum, and is standardly used as the reference SPE for analyses of the proxy-based extreme events. First, we have revised this reference event and reconstructed the rigidity spectrum of energetic particles for it (Figure 3), using the new state-of-the-art method, and have shown that earlier studies underestimated the strength of the event by a factor of up to three in the mid-energy range of a few GeV. By applying this newly reconstructed spectrum, we have evaluated the expected response of the cosmogenic isotopes 10 Be and 14 C to the reference event (Figures 4, 5 and 8). We found out that the reference SPE cannot be reliably identified in a single proxy record (cf., McCracken & Beer, 2015;Usoskin et al., 2006). On the other hand, we have estimated the sensitivity of the proxy method to extreme SPEs, using a multi-proxy approach, as shown in Figure 1. By combining several independent proxy records, the uncertainties can be reduced, so that a SPE by a factor 4 -5 stronger than the reference one, viz. an order of magnitude weaker than the strongest event of 774 AD, can be potentially detected in a multi-proxy record.
Since the uncertainties of an ice-core based 10 Be record are dominated by the local/regional pattern rather than by the measurement errors, combining results from different ice-cores may help significantly reducing the uncertainties. The uncertainties of the 14 C annual data are largely determined by the measurement errors, and represent the GCR and reference SPE1956 scenarios, respectively. The SPE curve is offset by 2 (at cm −2 sec −1 ) for better visibility. Panel b) depicts the modelled ∆ 14 C for the same period. The solid blue, dashed red and dotted green curves correspond to GCR, GCR+SPEx10 and GCR+SPEx20 scenarios, respectively. Error bars in panel b correspond to the best modern measurement uncertainty, 1.5 of ∆ 14 C (Güttler et al., 2015;Park et al., 2017).
thus, more precise and repetitive measurements of the same sample can improve the sensitivity to SPEs.
The ice-core based data of 10 Be is estimated to have approximately double sensitivity to SPEs with respect to 14 C mostly because of the strong attenuation of the fast SPE signal (a factor of three compared to the 11-yr solar cycle) and the spread of the SPE-produced 14 C over many years. However, this estimate does not account for possible dating errors in ice cores (Sigl et al., 2015), assuming that all proxy records are well dated and can be easily superposed. The real situation is not as good, and the signal can be smeared.
Overall, it is possible to increase the sensitivity of the proxy method to SPEs by an order of magnitude (Figure 1), filling the observational gap between events directly observed during the space era and the extreme events discovered recently, and to increase statistics of extreme events during the last millennia from 3 -4 known today to several tens. This will provide a solid basis for research in the field of extreme events, both for fundamental science, viz. solar and stellar physics, and practical applications, such as the risk assessments of severe space-based hazards for modern technological society.