ABSTRACT
The probability of occurrence of extreme solar particle events (SPEs) with proton fluence (>30 MeV) F30 ⩾ 1010 cm−2 is evaluated based on data on the cosmogenic isotopes 14C and 10Be in terrestrial archives covering centennial–millennial timescales. Four potential candidates with F30 = (1–1.5) × 1010 cm−2 and no events with F30 > 2 × 1010 cm−2 are identified since 1400 AD in the annually resolved 10Be data. A strong SPE related to the Carrington flare of 1859 AD is not supported by the data. For the last 11,400 years, 19 SPE candidates with F30 = (1–3) × 1010 cm−2 are found and clearly no event with F30 > 5 × 1010 cm−2 (50 times the SPE of 1956 February 23) has occurred. These values serve as observational upper limits on the strength of SPEs on the timescale of tens of millennia. Two events, ca. 780 and 1460 AD, appear in different data series making them strong candidates for extreme SPEs. We build a distribution of the occurrence probability of extreme SPEs, providing a new strict observational constraint. Practical limits can be set as F30 ≈ 1, 2–3, and 5×1010 cm−2 for occurrence probabilities ≈10−2, 10−3, and 10−4 yr−1, respectively. Because of the uncertainties, our results should be interpreted as a conservative upper limit on the SPE occurrence near Earth. The mean solar energetic particle (SEP) flux is evaluated as ≈40 (cm2 s)−1, in agreement with estimates from lunar rocks. On average, extreme SPEs contribute about 10% to the total SEP fluence.
Export citation and abstract BibTeX RIS
1. INTRODUCTION
Sporadic energy releases on the Sun can accelerate solar energetic particles (SEPs) in the corona and interplanetary medium. Such phenomena often lead to solar particle events (SPEs) observed at Earth; this is an important factor in solar–terrestrial relations, and specifically space weather (Mewaldt 2006; Vainio et al. 2009). A typical quantity for an SPE is the fluence of SEPs with energy above 30 MeV, F30. We have sufficient knowledge of SPEs over the space era since the mid-1950s (Smart et al. 2006), with only several events with F30 = (1–10) × 109 cm−2 and hundreds of weaker SPEs observed. However, it is important to know, both for purely theoretical aspects of solar/stellar physics and for technical applications, the statistics of extreme SPEs with F30 > 1010 cm−2. Such a study is possible only by using indirect proxy data. Some estimates have been obtained from measurements of cosmogenic isotopes in lunar rocks (Nishiizumi et al. 2009), but this gives only the average flux of SEPs over a very long scale without extracting individual SPEs. A list of potential SPEs over the last 500 years was proposed based on nitrate records in polar ice (McCracken et al. 2001), but this result is heavily debated (Wolff et al. 2008, 2012). Data on the cosmogenic nuclides 14C and 10Be measured in terrestrial archives may provide information on SPEs in the past (Lingenfelter & Hudson 1980; Usoskin et al. 2006; Webber et al. 2007), but this possibility has not been fully explored previously. Accordingly, the probability of extreme SPEs has remained exceedingly uncertain (Hudson 2010).
Here we establish a solid observational constraint on the distribution of extreme SPEs using presently available data sets on cosmogenic isotopes (14C and 10Be) measured in terrestrial archives with sufficient time resolution and quality, and modern models of their production in the atmosphere. We note that the result presented is based on terrestrial data and may not well represent the occurrence of solar events, whose geo-efficiency is also affected by the relative Sun–Earth attitude.
2. DATA SETS AND THE METHOD
2.1. Data
The data sets used are as follows:
-
1.IntCal09 Δ14C global series: 11000 BC – 1900 AD, 5 year time resolution (Reimer et al. 2009)
-
2.SB93 Δ14C global annual series: 1511–1900 AD (Stuiver & Braziunas 1993)
-
3.Dye3 10Be Greenland annual series: 1424–1985 AD (McCracken et al. 2004)
-
4.NGRIP 10Be Greenland annual series: 1389–1994 AD (Berggren et al. 2009)
-
5.
-
6.DF 10Be Dome Fuji Antarctic series: 695–1880 AD, quasi-decadal (Horiuchi et al. 2008)
-
7.
-
8.M12 Δ14C Japanese annual series: 750–820 AD, annual/biannual (Miyake et al. 2012)
2.2. Model Computations
In order to evaluate possible SPE signatures in the data, we used model computations of the isotope production by energetic particles, assuming instant injection of SEPs into the atmosphere and calculating the expected isotope response in terrestrial archives.
As the reference event, we considered an extreme SPE with a very hard spectrum (Tylka & Dietrich 2009; Usoskin et al. 2011) from 1956 February 23 (SPE56; Meyer et al. 1956). While SPE56 was the strongest observed ground-level enhancement (GLE), >4000% in the count rate of the Leeds NM, it had relatively modest fluence, F30 = 109 cm−2 (Shea & Smart 1990). However, sometimes SPEs with large fluence but a soft spectrum occur, e.g., a modest (only 10% at the polar Oulu NM) GLE event of 1972 August 4 (SPE72) with high fluence, F30 = 5 × 109 cm−2. Since cosmogenic isotopes are produced by the most energetic part of the SEP energy spectra (>1 GeV), we consider here the SPE56 scenario.
The response of 10Be to SPEs was calculated similar to Usoskin et al. (2006) but using an updated 10Be yield function (Kovaltsov & Usoskin 2010), the corresponding geomagnetic model, and an intermediate atmospheric mixing model (polar tropospheric and hemispheric stratospheric mixing—see Field et al. 2006; Heikkilä et al. 2009). The calculated 10Be production is 7.5 × 104 atoms cm−2 for the SPE56 scenario. Because of the 10Be stratospheric residence time (Heikkilä et al. 2009), a possible SPE peak in 10Be data can be 2–3 years long. Since present models cannot convert the 10Be production into concentrations in ice, they are typically assumed to be directly proportional to each other (e.g., McCracken et al. 2004). For the normalization we used the reference period 1850–1900 AD with moderate solar activity: the galactic cosmic ray (GCR) modulation potential 443 MV (Alanko-Huotari et al. 2007) was close to the mean Holocene value (Usoskin et al. 2007). Accordingly, the reference period 10Be production was 〈Q〉 = 3.45 × 10−2 atoms cm−2 s−1 or 1.1 × 106 atoms cm−2 yr−1. This value was scaled to the mean measured concentration for the reference period in each ice core (1.025, 1.86, 3.47, and 2.76×104 atoms g−1 for the Dye3, NGRIP, SP, and DF series, respectively). The GRIP series was normalized to the entire interval.
The expected radiocarbon response of Δ14C was calculated in two steps. First, the global production Q was computed using a new production model (Kovaltsov et al. 2012), yielding 2.9 × 106 atoms cm−2 for SPE56. Next, the instant atmospheric injection of 14C (Q) was passed through a 5-box carbon cycle model (Damon & Peristykh 2004), and the tropospheric Δ14C was calculated. Figure 1 shows the expected response in Δ14C for SPE56, which is extended over decades and is characterized by a sharp increase and exponential decay; the peak's FWHM is 15–20 years. Accordingly, we look for a signature like that in the Δ14C data. The SPE56 peak in the annual (5 year) Δ14C data is 0.35 (0.19–0.27) permille, viz., below measurement errors of ≈2 permille. Therefore, an SPE needs to be a factor XSPE56 = 10 greater than SPE56 to be observable in the radiocarbon data.
Figure 1. Response of the relative tropospheric concentration Δ14C to SEP56 (t = 0 shown by the arrow), with annual (gray curve) and 5 year (black) time resolutions, depending on the 5 year bin boundaries.
Download figure:
Standard image High-resolution imageAn SPE72-type soft event would require a 40 times larger F30, with respect to SPE56, to produce the same amount of cosmogenic isotopes. Accordingly, the estimates obtained with the reference SPE56 spectrum should be enhanced 40-fold to correspond to a soft-spectrum SPE72 scenario.
3. EVALUATING SEP FLUENCE IN THE PAST
First, we analyzed data series with annual resolution (2–4 in Section 2.1), looking for potential signatures of SPEs.
For the NGRIP series we looked for peaks with magnitude ⩾1.3 × 104 atoms g−1 and duration ⩽3 years. Seven candidates were selected: 1436, 1460, 1650, 1719, 1810, 1816, and 1965 AD. For the 10Be Dye3 series, five candidates (magnitude ⩾0.6 × 104 atoms g−1, duration ⩽3 years) were selected: 1462, 1479, 1505, 1512, and 1603 AD. For the 14C SB93 series, no suitable candidates (peak with sharp rise and gradual decay) were found. Since the candidates found do not coincide in time, we performed a cross-check. For each candidate found in series A we calculated, using the model (Section 2.2), the corresponding isotope production Q and the expected response in series B, and checked if they were consistent with the data within ±2 years and vice versa. For example, a 1436 AD peak in the NGRIP series must be accompanied by a strong peak in the Dye3 series, which was not observed. Also, an NGRIP 1965 AD peak is rejected by direct cosmic ray observations. Using this cross-check, seven candidates can be excluded: 1436, 1650, 1816, and 1965 AD in NGRIP and 1479, 1512, and 1603 AD in the Dye3 series. The SB93 series is insensitive to checking the 10Be-based candidates. Finally, five peaks pass through the check as potential SPE signatures (first block in Table 1).
Table 1. Identified SPE Candidates: Approximate Year, Data Set, and Fluence F30 (cm−2) Evaluated for the SPE56 Scenario
| SPE Year | Series | F30 (XSPE56) |
|---|---|---|
| 1460–1462 ADa | NGRIP(1460) | 1.5 × 1010 (15) |
| Dye3 (1462) | 9.7 × 109 (10) | |
| 1505 AD | Dye3 | 1.3 × 1010 (13) |
| 1719 AD | NGRIP | 1 × 1010 (10) |
| 1810 AD | NGRIP | 1 × 1010 (10) |
| 8910 BC | IntCal09 | 2.0 × 1010 (20) |
| 8155 BC | IntCal09 | 1.3 × 1010 (13) |
| 8085 BC | IntCal09 | 1.5 × 1010 (15) |
| 7930 BC | IntCal09 | 1.3 × 1010 (13) |
| 7570 BC | IntCal09 | 2.0 × 1010 (20) |
| 7455 BC | IntCal09 | 1.5 × 1010 (15) |
| 6940 BC | IntCal09 | 1.1 × 1010 (11) |
| 6585 BC | IntCal09 | 1.7 × 1010 (17) |
| 5835 BC | IntCal09 | 1.5 × 1010 (15) |
| 5165 BC | GRIP | 2.4 × 1010 (24) |
| 4680 BC | IntCal09 | 1.6 × 1010 (16) |
| 3260 BC | IntCal09 | 2.4 × 1010 (24) |
| 2615 BC | IntCal09 | 1.2 × 1010 (12) |
| 2225 BC | IntCal09 | 1.2 × 1010 (12) |
| 1485 BC | IntCal09 | 2.0 × 1010 (20) |
| 95 AD | GRIP | 2.6 × 1010 (26) |
| 265 AD | IntCal09 | 2.0 × 1010 (20) |
| 780 ADa | IntCal09 | 2.4 × 1010 (24) |
| M12 | 4 × 1010 (40)b | |
| DF | 4.5 × 1010 (45)b | |
| 1455 ADa | SP | 7 × 1010 (70)b |
Notes. Scaling to SPE56 is given in parentheses. aDiscussed separately. bOverestimate.
Download table as: ASCIITypeset image
Only one candidate (1460–1462 AD) is present in both series (discussed later).
Next, we considered data with rougher time resolutions. In each of the 10Be data series (5–7 in Section 2.1) with quasi-decadal resolution, we searched for distinctive single peaks. Two candidates were identified in the GRIP series: ca. 5165 BC and 95 AD. One candidate was identified in the SP series: ca. 1455 AD, probably related to the ca. 1460 event in the NGRIP and Dye3 annual series. Two candidates were identified in the DF series: ca. 780 AD and 1805 AD. A search for signatures (sharp ⩽10 year increase, FWHM 15–30 years, magnitude >2 permille) in the IntCal09 14C series yields many candidates (see an example in Figure 3(a)). We also made a cross-check of all the candidates between different series. For instance, the 1805 AD peak in the DF series was rejected based on the annual NGRIP and Dye3 data. About half of the candidates in the IntCal09 series after 7000 BC were rejected based on the 10Be data. Therefore, half of the identified candidates before 7000 BC may also be spurious.
Finally, in Table 1 we list the candidates passing the cross-check. One can see that 10 times SPE56, XSPE56 = 10, is the detection limit in cosmogenic isotope data. Only two candidates appear in more than one series: ca. 780 AD and 1460 AD, which are discussed in detail below.
3.1. Event of ca. 1460 AD
Several series analyzed show a significant peak around 1460 AD (Table 1). This was noticed previously and ascribed to a very strong SPE or to a supernova explosion (e.g., Berggren 2009; McCracken et al. 2004; Delaygue & Bard 2011). We assume here that it is an SPE signature and evaluate its parameters. Annual NGRIP and Dye3 series depict distinct peaks in 1460 and 1462 AD, respectively, which roughly agree with each other, and require F30 ≈ 1010 protons cm−2 for the SPE56 scenario. The more roughly resolved SP series depicts a very strong peak ca. 1455 AD, requiring a fluence seven times greater than that implied from the NGRIP/Dye3 data. None of the other series (GRIP, DF, or IntCal09) show peaks around that date allowing for evaluation of the upper limit. Figure 2(a) shows a cross-check of this event versus the DF series. While the event strength evaluated from the NGRIP and Dye3 series is consistent with no clear signal in the DF data, the huge fluence implied from the SP data contradicts the DF data (cf. Delaygue & Bard 2011). A similar analysis of the GRIP data also suggests that the SP peak is too high. Figure 2(b) shows the IntCal09 data along with signals expected from the NGRIP/Dye3 and SP data, respectively, for the ca. 1460 AD candidate. Because of the steep gradient in Δ14C in 1430–1470, it is impossible to distinguish a small 2.4 permille peak implied by the Dye3/NGRIP data, but the high signal implied by the SP peak apparently contradicts the data. Thus, while all the data sets analyzed are consistent with a hypothesis of a strong SPE ca. 1460 AD, its fluence is likely about 1010 protons cm−2 for the SPE56 scenario. The very high fluence implied by the SP data contradicts the other data sets.
Figure 2. Test of a possible SPE ca. 1460 AD. (a) DF 10Be series, with points corresponding to the SEP production estimated from the NGRIP, Dye3, and SP data series (see the legend). (b) IntCal09 Δ14C series along with the expected signal estimated from the Dye3/NGRIP and SP data series.
Download figure:
Standard image High-resolution image3.2. Event of ca. 780 AD
Another interesting candidate for a strong SPE is ca. 780 AD as manifested through distinct peaks in the DF, IntCal09, and M12 data series (Figure 3). According to these data sets, the event was 25–50 times stronger than SPE56 with fluence F30 = (2–5) × 1010 protons cm−2. The best fit for the IntCal09 data (Figure 3(a)) is obtained for an XSPE56 = 24 event starting ca. 780 AD. The M12 data (Figure 3(b)) are a better fit with an XSPE56 = 40 event starting in 775 AD, while the DF data (Figure 3(c)) are consistent with an XSPE56 = 45 event that occurred between 780 and 790 AD. However, the GRIP series (Figure 3(d)) depicts no peak in that period, which indicates that a fluence greater than F30 = 3 × 1010 protons cm−2 is inconsistent at the 0.03 significance level. Applying the cross-check, we found that the fluence implied by the IntCal09 data (XSPE56 = 24—see Table 1) is consistent with the GRIP data, while both XSPE56 = 40 and 45 based on the DF and M12 data, respectively, yield too high a peak in the GRIP series. Accordingly, we conclude that the event of 780 AD is a possible candidate for a strong SPE with a consensus value of F30 ≈ 3 × 1010 protons cm−2.
Figure 3. Test of a possible SPE ca. 780 AD. Open dots and black stars correspond to the modeled response to an SEP event equal to XSPE56 = 24 and XSPE56 = 45, respectively. (a) IntCal09 Δ14C data series. (b) M12 Δ14C data series. (c) DF 10Be data series. (d) GRIP 10Be data series.
Download figure:
Standard image High-resolution image3.3. Carrington Event of 1859 AD
The Carrington event of 1859 AD is often considered to be an extreme SPE, with F30 = 1.8 × 1010 cm−2 as estimated from the nitrate record in the Greenland Summit core (Dreschhoff & Zeller 1998; McCracken et al. 2001; Shea et al. 2006). Although none of the series analyzed here depict a peak around 1859, we checked if this proposed SPE is consistent with the cosmogenic data. Using the GCR modulation reconstruction (Alanko-Huotari et al. 2007) and the 10Be production model (Kovaltsov & Usoskin 2010), we calculated the expected 10Be concentration along with the additional production for 1859 AD, applying SEP parameters for the Carrington SPE according to Section 8 of McCracken et al. (2001). The model computation was then confronted with the annual data of two Greenland series, NGRIP and Dye3, in Figure 4. The 10Be peak, expected from the nitrate data, is too strong, contradicting the observed data in the two Greenland sites. Radiocarbon data and decadal 10Be series cannot resolve the Carrington peak. Therefore we conclude that the cosmogenic data do not support the hypothesis of a very strong SPE related to the Carrington flare.
Figure 4. Testing a possible SEP event ca. 1859 AD. The gray line depicts a model GCR signal plus the addition of the annual production due to the Carrington 1859 AD event according to McCracken et al. (2001). (a) Time series of the annual NGRIP 10Be data, with 15% error bars. (b) Time series of the annual Dye3 10Be data, with 15% error bars.
Download figure:
Standard image High-resolution image4. DISCUSSION
We propose, by analyzing data from different cosmogenic isotope records and performing cross-checks, a list of candidates for strong SPEs in the past (Table 1). The list is somewhat uncertain, as some events may be missing or some spurious peaks falsely identified. However, the result is quite consistent in a statistical sense.
Four candidates with an XSPE56 factor of 10–15 are identified in the annual 10Be series NGRIP and Dye3 for the last 600 years (first block of Table 1). Events with XSPE56 < 10 cannot be reliably identified. We can also securely say that events with XSPE56 > 20 are not observed during that period (cf. Figure 4). In particular, a strong SPE related to the Carrington flare in 1859 contradicts these data.
We evaluate the probability p of occurrence of such events, applying the Poisson distribution (assuming that SPEs are mutually independent). For example, for four events observed during 600 years, the probability is p = 0.0077+0.0073− 0.0045 yr−1 (90% confidence interval), which is higher than the naively taken 4/600 = 0.0067 yr−1. If no event with F30 > 2 × 1010 cm−2 is observed over 600 years, the corresponding one-sided 90% confidence interval is p = 0–0.0027 yr−1, with median probability 0.0012 yr−1 (once per 850 years).
The last 4–6 centuries cover the full range of solar activity, from Grand minima to the maximum (Usoskin 2008) serving as an archetype of the solar variability. Therefore, we expect that this result is consistent over longer scales. On the multi-millennial scale (the Holocene—11 millennia), 19 candidates are identified (second block in Table 1) with XSPE56 = 10–30. A few candidates with XSPE56 > 50 were rejected by the cross-check. This gives the average occurrence rate of such SPEs as roughly two per millennium. Moreover, we can securely say that during the Holocene there were no events with XSPE56 > 50–100, placing a strong upper limit on the SPE strength. We find no apparent relation between the occurrence of SPE candidates and the solar activity level (Usoskin et al. 2007).
We summarize our findings in a plot (Figure 5) of the integral probability of strong SPE occurrence. The X- and Y-axes give the F30 fluence and the probability of occurrence of SPEs with the >30 MeV proton annual fluence exceeding F30 for different data sets. The measured annual fluences for the space era (1956–2008) are shown by open triangles (Shea & Smart 1990; M. Shea 2012, private communication). Since the space era coincides with the unusually active Sun (Solanki et al. 2004), these probabilities may be higher than those for the typical medium activity Sun. The black triangle reflects the fact that no SPE with F30 > 1010 cm−2 was observed during 53 years. The open star corresponds to the four SPE candidates from the annual 10Be data (first block in Table 1), while the filled star corresponds to no SPE with F30 > 2 × 1010 cm−2 found during 600 years. Open circles represent the candidates found in the rougher time series (second block in Table 1), grouped into four points: F30 > 1, 1.5, 2, and 2.5× 1010 cm−2, respectively. Since the lower fluence value is probably underestimated because of the detection threshold, we ignore the first open circle, using instead the star-symbolled point from the annual 10Be data sets. The filled circle corresponds to no event with F30 > 5 × 1010 cm−2 found over 11 millennia. One can see that there is a break in the distribution at about F30 = 5 × 109 cm−2. This break (cf. McCracken et al. 2001; Hudson 2010) may be related to a streaming limit (Reames 1999) and/or limitations in the energetics of solar flares/coronal mass ejections (Fletcher et al. 2011). Our new results put a stronger constraint on the distribution. For example, an SPE with XSPE56 > 20 is expected only once per millennium, and no XSPE56 > 50 event is expected over tens of millennia. We try to fit the tail of the distribution (F30 > 7 × 109 cm−2) by two simple models: exponential (p∝exp (− 0.33 × 10−9 · F30)) and power law (p∝F−430)—see the gray lines in the figure. The data do not allow us to distinguish between the two shapes.
Figure 5. Probability of occurrence of the annual fluence (>30 MeV) exceeding the given value F30, as evaluated from the data for the space era 1956–2008 (triangles), cosmogenic isotope annual data (stars), and cosmogenic isotope decadal data (circles). Open symbols correspond to the statistics of "observed" data, while filled symbols denote the conservative upper limit (no events greater than X observed during Y years). Error bars correspond to the 90% confidence level evaluated assuming Poisson statistics. Gray dotted and solid lines depict the best-fit exponent and power-law approximations for F30 > 7 × 109 protons cm−2.
Download figure:
Standard image High-resolution image
We also compare our results with the average SEP flux over timescales of thousands to millions of years, estimated from lunar rocks, which range from 21–56 (cm2 s)−1 (e.g., Reedy 1999; Nishiizumi et al. 2009). Time averaging of the results from Figure 5 using the above fits for the tail yields an average flux ≈38 (cm2 s)−1, which is composed of about 35 (cm2 s)−1 for the space era record (e.g., Shea & Smart 2002; Reedy 2012) with an addition of 3.2 (cm2 s)−1 due to the distribution tail based on cosmogenic isotope data (both exponent and power law give similar results), suggesting that extreme SPEs contribute only about 10% to the total SEP fluence. This is totally consistent with assessments based on lunar rocks, giving independent support to the present results. For an SPE with a softer energy spectrum, like SPE72, the tail addition would be 40 times greater (Section 2.2), leading to an average SEP flux ≈150 (cm2 s)−1, in contradiction with the lunar rock data. Therefore, extreme SPEs found in the cosmogenic isotope records must have (on average) hard energy spectra.
5. CONCLUSIONS
We have evaluated the probability of occurrence of extreme SPEs based on data on the cosmogenic isotopes 14C and 10Be in terrestrial archives, spanning the timescale from centuries to 11 millennia (Figure 5). We identified four potential candidates for SPEs with F30 = (1–1.5) × 1010 cm−2 and show that no event with F30 > 2 × 1010 cm−2 existed over the last 600 years using annually resolved 10Be data. In particular, the extreme Carrington SPE of 1859 AD contradicts these data. From more roughly resolved data we identified 19 SPE candidates (Table 1) with F30 = (1–3) × 1010 cm−2, and clearly no event with F30 > 5 × 1010 cm−2, over the last 11,400 years. Two events, ca. 780 AD and 1460 AD, appear in different series making them strong candidates for extreme SPEs. This gives a new strict observational constraint on the occurrence probability of extreme SPEs. Practical limits can be set as F30 ≈ 1, 2–3, and 5 × 1010 cm−2 (10, 20–30, and 50 times greater than SPE56) for the occurrence probabilities of 10−2, 10−3, and 10−4 yr−1, respectively. The mean SEP flux is found as ≈40 (cm2 s)−1, in agreement with estimates from lunar rocks. On average, extreme SPEs contribute about 10% to the total SEP fluence.
We note that the present result tends to represent an upper limit for SPE occurrence, since we explicitly assume that every peak in one data series, consistent with other series, is a signature of an SPE. This may be not exactly correct, as some of the peaks may be spurious. Accordingly, our results should be interpreted as a conservative upper limit on the SPE occurrence near Earth. Application of the results to the Sun is not straightforward, especially for more energetic events with low statistics, because the propagation of SEPs in interplanetary space may greatly affect the geo-efficiency of SPEs and, accordingly, their ability to become detectable in the cosmogenic isotope data series. Given these uncertainties, the present results should be considered with a precision of up to an order of magnitude.
We thank M. A. Shea for the data on observed annual SEP fluences. G.A.K. was partly supported by Program No. 22 of the Presidium RAS and also acknowledges support from the Academy of Finland and Suomalainen Tiedeakatemia (Väisälä foundation).