Abstract
We present the results of an investigation of the relation between space-weather parameters and cosmic ray (CR) intensity modulation using algorithm-selected Forbush decreases (FDs) from Moscow (MOSC) and Apatity (APTY) neutron monitor (NM) stations during solar cycle 23. Our FD location program detected 408 and 383 FDs from MOSC and APTY NM stations respectively. A coincident computer code employed in this work detected 229 FDs that were observed at the same Universal Time (UT) at the two stations. Out of the 229 simultaneous FDs, we formed a subset of 139 large FDs(%) ≤ − 4 at the MOSC station. We performed a two-dimensional regression analysis between the FD magnitudes and the space-weather data on the two samples. We find that there were significant space-weather disturbances at the time of the CR flux depressions. The correlation between the space-weather parameters and decreases in galactic cosmic ray (GCR) intensity at the two NM stations is statistically significant. The implications of the present space-weather data on CR intensity depressions are highlighted.
Export citation and abstract BibTeX RIS
1. Introduction
Forbush decreases (FDs) are non-periodic short-term (hours to days) variabilities in the observed intensity of galactic cosmic ray (GCR) flux. This phenomenon that is sometimes regarded as a rapid and transient depression of 2%–30% in GCR intensity (Forbush 1958; Belov 2009; Patra et al. 2011) was first reported by an American Physicist, Scott E. Forbush, over 80 years ago (Forbush 1938). FDs, caused by corotating high-speed solar wind streams (CSWSs) from coronal holes (CHs) which corotate with the Sun, are recurrent FDs, while events initiated by coronal mass ejections (CMEs) and their interplanetary medium signature (interplanetary coronal mass ejections, ICMEs) are regarded as non-recurrent or sporadic FDs (Cane et al. 1993; Richardson 2004; Richardson & Cane 2011; Belov et al. 2014; Richardson 2018). Recurrent FDs typically are symmetric and have small magnitudes while large magnitudes and asymmetric depressions characterize the non-recurrent FD group (Lockwood 1971; Belov et al. 2009; Melkumyan et al. 2019). Generally, solar wind disturbances in the form of CMEs, ICMEs, corotating interaction regions (CIRs) and CSWS structures determine cosmic ray (CR) intensity depressions on Earth (Cane 2000; Dumbović et al. 2012; Okike & Nwuzor 2020).
A large number of research works has intensively discussed the FD phenomenon since its first observation, yet the origin, nature or connection of the phenomenon with space-weather structures is still a subject of interest among CR scientists (Belov 2009; Okike & Nwuzor 2020). Associations between FD amplitude, intensity of the interplanetary magnetic field (IMF), solar wind speed (SWS) and geomagnetic disturbance storm time (Dst) index have been reported (e.g. Burlaga et al. 1984; Cane et al. 1993; Belov 2009; Belov et al. 2001a; Lingri et al. 2016; Okike & Nwuzor 2020; Alhassan et al. 2021) (Alhassan et al. (2021) will be referred to as Paper I hereafter). Belov et al. (2001a) noted that FDs and geomagnetic storms are both traced to disturbances in the interplanetary medium. By extension, disturbances in the solar wind, magnetosphere and FDs are triggered by the same processes in the Sun. The study of these disturbances is important in our understanding of the dynamics of the solar-terrestrial environment since such activities can potentially result in life threatening issues ranging from satellite damage, communication failure to navigational difficulties (Joselyn & McIntosh 1981; Rathore et al. 2011). This underscores the rationale for examining the impact of solar-geomagnetic characteristics on GCR intensity reductions.
While SWS, IMF, Dst and other solar activity data could be sourced from some common websites (e.g., https://omniweb.gsfc.nasa.gov/html/owdata.html, http://wdc.kugi.kyoto-u.ac.jp/, http://cdaw.gsfc.nasa.gov, ftp.ngdc.noaa.gov), FDs have to be calculated or sourced from http://spaceweather.izmiran.ru/eng/fds2005.htm. Only the IZMIRAN group has a comprehensive FD catalog obtained with the global survey method (GSM). The GSM involves computer codes and manual methods. The technique is better described as semi-automated since a lot of manual work is incorporated (Belov et al. 2018). Apart from the IZMIRAN group, other investigators have employed a few FDs either from literature (e.g. Kristjánsson et al. 2008; Lagouvardos et al. 2009; Laken et al. 2011) or from the manual technique (e.g. Todd & Kniveton 2001; Kane 2010; Dragić et al. 2011). The manual method has historically faced a lot of disadvantages such as the time consuming nature of the process, difficulties in identification of low-amplitude FD events, inability to detect a large number of FDs and absence of time domain in the model equation used to calculate FD event amplitude (e.g. Okike 2020b). FD has been found to be the most spectacular variability in GCR flux (Musalem-Ramirez et al. 2013), hence, accurate detection of FDs in a large volume of CR data requires a sophisticated technique. A computer software technique can address the inherent limitations associated with the manual FD selection method.
Emerging efforts to fully automate FD selection are already appearing in the literature (Okike & Umahi 2019; Okike 2020b, a; Okike & Nwuzor 2020; Okike et al. 2021b; Light et al. 2020). With some technical improvements, we extend the code employed by Okike & Umahi (2019) in our attempt to investigate the relation between solar-weather parameters and program selected FDs during solar cycle 23. FDs within this period are unique for study in several ways. It has been thoroughly investigated and is marked by high solar activity (Okike et al. 2021a). Observations during this period are reported to be related to several CHs that trigger high speed wind streams associated with many energetic events, especially during its declining phase (e.g. Zerbo et al. 2013; Lingri et al. 2016).
2. Data
Moscow (MOSC) and Apatity (APTY) CR neutron monitor (NM) stations are among the oldest observatories that have continuous records of such observations. MOSC has online data since 1958 to date while APTY has data since 1961 till date. MOSC NM station is characterized by geographic location of 55.47°N, 37.32°E, geomagnetic cutoff rigidity of 2.43 GV and altitude of 200 m with the 24NM64 detector while APTY is located at 67.50°N, 33.30°E, rigidity of 0.57 GV and altitude of 177 m with the 18NM64 detector. CR pressure-corrected daily averaged data detected at MOSC and APTY NM stations can be downloaded from the IZMIRAN website:http://cr0.izmiran.ru/common. In this study, IMF, SWS and Dst index data are obtained from https://omniweb.gsfc.nasa.gov/html/owdata.html.
3. Algorithm FD Selection Method
3.1. FD Location Program
In this section, we attempt to address some of the weaknesses associated with the manual FD selection technique highlighted in Section 1. The computer code briefly described in this work is an extension of the Okike & Umahi (2019) software. The technique in Okike & Umahi (2019) consists of modules for Fourier transformation where the machine searches for FDs from the transformed CR data. The present algorithm, written in the R programming language for statistical computing, a non-commercial software package (Team 2015), accepts CR raw data as its input signal instead of the high frequency input signal in Okike & Umahi (2019). A similar equation implemented by several authors (e.g. Harrison & Ambaum 2010) that manually calculates FD magnitude is used in the present code for data normalization.
The software is equipped to select FDs from both raw and normalized CR data. It is programmed to search for depressions/turning points as well as the time of the depressions in the raw CR count. The depressions are signatures of FDs. Two subroutines in the software calculate the FD magnitude and the time of occurrences simultaneously. The baseline adopted determines the number of FDs that may be selected in a given period (see Okike & Umahi 2019). Utilizing a very small baseline (CR(%) ≤ − 0.01), as attempted by Okike (2020b), Okike et al. (2021a) and Paper I, the software identified 408 and 383 FDs respectively from MOSC and APTY NM stations.
The IZMIRAN group selected a total of 1346 Forbush Effect (FE) 10 GV CRs with GSM during this period in which they assimilated CR counts from over 50 NM stations spread across the globe. The FDs detected at MOSC and APTY translate to 30% and 29% respectively of the number selected by the IZMIRAN group. The amplitude of the smallest FE is −0.2% whereas the size of the least FD identified by the current code is −0.01%. The fully automated FD technique employed here has been demonstrated to be more efficient than the semi-automated approach. Our software tracks FD minimum depression time while the GSM technique is based on the event onset time. Setting a large threshold defined as (CR(%) ≤ − 4), the number of FE is 75 while FD is 139. This demonstrates that the large number of differences between the two FD catalogs should lie majorly in the small FD population. The observed differences in number of FDs selected with both techniques could be due to a number of factors that may include: altitude, rigidity, atmospheric depth, pressure, temperature, relative humidity, local wind speed, the rotation of the Earth with respect to the acceptance cone of the detectors, latitudinal effects, instrumental variations, a stationʼs sensitivity to CR modulation, equatorial anisotropy, North-South anisotropy, geomagnetic variations, snow, limited cone of acceptance, spurious modulation, magnetospheric effects, differences in data resolution, etc. (see Belov et al. 2018; Okike & Nwuzor 2020).
FD magnitude and time of occurrences selected with the location code are presented in Table 1. In the Table, "S/N" stands for serial number, "Date" for FD dates at MOSC and APTY stations, FDMOSC(%) for FD at MOSC and FDAPTY(%) for FD at APTY.
Table 1. Selected FDs at MOSC and APTY Stations from 1996–2005
| S/N | Date | FDMOSC(%) | Date | FDAPTY(%) |
|---|---|---|---|---|
| 1 | 1998–05–02 | −1.78 | 1998–05-02 | −1.34 |
| 2 | 1998–08–27 | −3.16 | 1998–05–04 | −1.57 |
| 3 | 1998–09–25 | −2.16 | 1998–08–27 | −2.08 |
| 4 | 1999–01–24 | −2.70 | 1999–02–18 | −2.96 |
| 5 | 1999–02–18 | −3.46 | 1999–08–22 | −1.61 |
| 6 | 1999–06–27 | −0.52 | 1999–08–25 | −0.53 |
| 7 | 1999–08–20 | −2.04 | 1999–09–13 | −0.27 |
| 8 | 1999–08–22 | −2.57 | 1999–09–16 | −1.05 |
| 9 | 1999–08–25 | −1.72 | 1999–09–20 | −0.68 |
| 10 | 1999–09–05 | −0.89 | 1999–09–25 | −0.01 |
| 11 | 1999–09–07 | −0.64 | 1999–09–29 | −1.50 |
| 12 | 1999–09–09 | −0.95 | 1999–10–12 | −1.57 |
| 13 | 1999–09–16 | −1.76 | 1999–10–15 | −3.16 |
| 14 | 1999–09–18 | −1.68 | 1999–10–17 | −3.34 |
| 15 | 1999–09–21 | −1.86 | 1999–10–22 | −3.26 |
| .. | ... | ... | ... | ... |
The full table is available at http://www.raa-journal.org/docs/Supp/ms4947Table1.pdf.
3.2. Simultaneous FD list
Simultaneity of FDs has been conceptualized in various ways by different investigators. An FD event is said to be simultaneous when the main phase profile overlaps in Universal Time (UT) (e.g. Oh et al. 2008; Lee et al. 2013). The complete FD event structure (onset, main phase, point of minimum reduction and recovery phase) has been employed by Okike & Collier (2011) in their attempt to define FD event simultaneity. In this submission, we define FD simultaneity with respect to the event time of minimum. The time of minimum of simultaneous FDs would be detected on the same date at MOSC and APTY NMs.
To select FDs that are coincident at MOSC and APTY NM stations, we relied on a simple coincident computer code developed by Okike & Nwuzor (2020). TheFD coincident algorithm is a simple program written to select simultaneous as well as non-simultaneous FDs at two or more NM stations. Event time or magnitude may be input as the key search terms. These coincident FD events were selected with respect to the magnitude of FDs at APTY station. The input data to the coincident algorithm are all the FDs at each of the stations (408 at MOSC and 383 at APTY, see Table 1). The output demonstrates that 229 FDs occurred at the same UT at the two stations. Out of all the 229 simultaneous FDs selected by the coincident program, we formed a subset of 139 large events at MOSC defined as FD(%) ≤ − 4 (see Okike 2020a). These datasets are presented in Tables 2 and 3 alongside their associated solar-geophysical parameters. The columns are organized as follows: "S/N" represents serial number, "Date" stands for date of CR intensity depression, FDMOSC(%) and FDAPTY(%) represent MOSC and APTY FDs respectively, "IMF" signifies interplanetary magnetic field, "SWS" corresponds to solar wind speed and "Dst" is for disturbance storm time index.
Table 2. Simultaneous FD Events from MOSC and APTY and Associated Parameters during 1996–2005
| S/N | Date | FDMOSC(%) | FDAPTY(%) | IMF | SWS | Dst |
|---|---|---|---|---|---|---|
| 1 | 1998-05-02 | −1.78 | −1.34 | 14.50 | 601 | −36 |
| 2 | 1998-08-27 | −3.16 | −2.08 | 14.10 | 630 | −129 |
| 3 | 1999-02-18 | −3.46 | −2.96 | 17.10 | 599 | −84 |
| 4 | 1999-08-22 | −2.57 | −1.61 | 6.00 | 428 | −27 |
| 5 | 1999-08-25 | −1.72 | −0.05 | 7.70 | 538 | −15 |
| 6 | 1999-09-16 | −1.76 | −1.05 | 6.20 | 572 | −46 |
| 7 | 1999-09-25 | −0.75 | −0.01 | 5.60 | 409 | −16 |
| 8 | 1999-09-29 | −1.53 | −1.50 | 6.80 | 539 | −30 |
| 9 | 1999-10-12 | −2.20 | −1.57 | 7.30 | 578 | −48 |
| 10 | 1999-11-01 | −2.93 | −2.86 | 7.20 | 440 | −15 |
| 11 | 1999-11-09 | −2.43 | −2.37 | 6.30 | 615 | −46 |
| 12 | 1999-11-18 | −3.72 | −4.27 | 6.00 | 541 | −31 |
| 13 | 1999-11-20 | −3.63 | −4.83 | 8.10 | 443 | −16 |
| 14 | 1999-12-02 | −3.29 | −3.64 | 10.00 | 344 | 13 |
| 15 | 1999-12-13 | −7.47 | −7.49 | 11.40 | 489 | −46 |
| .. | ... | ... | ... | ... | ... | ... |
The full table is available at http://www.raa-journal.org/docs/Supp/ms4947Table2.pdf.
Table 3. Simultaneous Large FDsMOSC(%) ≤ − 4 Events and Associated Solar-geophysical Parameters during 1996-2005
| S/N | Date | FDMOSC(%) | IMF | SWS | Dst |
|---|---|---|---|---|---|
| 1 | 2000–10–01 | −4.01 | 4.30 | 418 | −37 |
| 2 | 2001–05–25 | −4.02 | 6.60 | 557 | 5 |
| 3 | 2003–12–28 | −4.20 | 9.70 | 508 | −14 |
| 4 | 2001–08–03 | −4.22 | 7.20 | 405 | 0 |
| 5 | 2001–08–06 | −4.24 | 7.00 | 440 | −18 |
| 6 | 2004–08–01 | −4.28 | 6.50 | 471 | −25 |
| 7 | 2000–05–15 | −4.31 | 9.10 | 414 | 7 |
| 8 | 2003–05–22 | −4.33 | 7.00 | 493 | −42 |
| 9 | 2000–09–29 | −4.34 | 5.50 | 378 | −19 |
| 10 | 2000–05–03 | −4.35 | 6.20 | 520 | −12 |
| 11 | 1999–12–27 | −4.35 | 7.90 | 410 | 2 |
| 12 | 2003–09–04 | −4.40 | 8.50 | 612 | −13 |
| 13 | 2001–12–07 | −4.49 | 6.70 | 459 | −19 |
| 14 | 2000–04–24 | −4.53 | 9.40 | 485 | −25 |
| 15 | 2003–09–12 | −4.54 | 4.70 | 593 | −6 |
| .. | ... | ... | ... | ... | ... |
The full table is available at http://www.raa-journal.org/docs/Supp/ms4947Table3.pdf.
4. Results and Discussion
The FD catalogs analyzed in this work are all simultaneous FDs from MOSC and APTY CR stations (229 FDs), and large FD(%) ≤ − 4 (139 FDs) at MOSC with their associated solar wind data, IMF and Dst. The regression and correlation results are presented in Tables 4 and 5. The columns of Tables 4 and 5 are organized as: S/N for serial number, parameters for each of the two continuous variables, R2 indicates coefficient of determination, r is Pearsonʼs product moment correlation coefficient and p-value is chance probability. The analysis of small amplitude FD(%) > − 4 will be taken up in future investigation.
Table 4. Correlation Results of All Simultaneous FDMOSC, FDAPTY and Associated Characteristics
| S/N | Parameter | R2 | r | p-values |
|---|---|---|---|---|
| 1 | FDMOSC-IMF | 0.13 | 0.36 | 3.16 × 10−8 |
| 2 | FDMOSC-SWS | 0.17 | 0.41 | 1.27 × 10−10 |
| 3 | FDMOSC-Dst | 0.16 | 0.40 | 6.51 × 10−10 |
| 4 | FDAPTY-IMF | 0.12 | 0.35 | 8.53 × 10−8 |
| 5 | FDAPTY-SWS | 0.11 | 0.34 | 1.45 × 10−7 |
| 6 | FDAPTY-Dst | 0.14 | 0.37 | 6.20 × 10−9 |
Table 5. Correlation Results of Large FDsMOSC(%) ≤ − 4 and Associated Parameters
| S/N | Parameter | R2 | r | p-values |
|---|---|---|---|---|
| 1 | FDMOSC-IMF | 0.14 | 0.37 | 5.37 × 10−6 |
| 2 | FDMOSC-SWS | 0.23 | 0.48 | 2.13 × 10−9 |
| 3 | FDMOSC-Dst | 0.13 | 0.36 | 1.34 × 10−5 |
For the two FD catalogs (N = 229, N = 139), we found a connection between FD versus IMF. The R2, r and associated p-values for all the simultaneous events at MOSC (Fig. 1(a)), APTY (Fig. 1(b)) and strong FDs at MOSC (Fig. 2(a)) are 0.13, 0.12 and 0.14, −0.36, −0.35 and −0.37, and 3.16 × 10−8, 8.53×10−8 and 5.37 × 10−6 respectively. They all show a high level of significance at the 95% confidence level. The results show further that 13%, 12% and 14% of CR intensity variations observed in the two catalogs could be linked to the effect of IMF intensity. While it may be inferred from Table 2 that the magnitudes of the simultaneous FDs vary appreciably between the two stations, the close variation rate with respect to IMF observed (see Table 4) at MOSC and APTY CR stations seems to suggest that the simultaneity of FDs may not be dependent on the point of observation. Nevertheless, several other parameters such as CME speed or transit speed, magnetospheric effects, sunspot number, current sheet tilt angle, solar magnetic turbulence level, CR anisotropy, heliospheric magnetic sector, rigidity, differences in the local time of the two stations, instrumental variations and so on, which affect CR time-intensity variations, should also be tested before reaching any definitive conclusions (Cliver & Cane 1996; Singh et al. 1997; Smith 1990; Wibberenz et al. 2001; Owens et al. 2014; Okike 2019, 2020a; Okike & Nwuzor 2020). In some past articles that analyzed the relationship between FDs and these physical parameters, the magnitude and timing of Forbush events were manually estimated. The present large event catalog is an indication that only FD subsamples (see Laken et al. 2012, for a detailed discussion on the bias implications of a small sample of FD events) were employed in some of the previous submissions that manually calculated FD event magnitude and timing (e.g. Barouch & Burlaga 1975; Richardson 2004; Kane 2010). These automated event catalogs provide an opportunity for re-assessments as well as statistical analyses of the previously reported relations. Some of these investigations may be the focus of a future publication.
Fig. 1. Plots of MOSC, APTY FD and related solar wind parameters and geomagnetic activity index.
Download figure:
Standard image
Fig. 2. Plots of large FDsMOSC(%) ≤ − 4 and related solar wind parameters and geomagnetic storm index.
Download figure:
Standard imageOne of the major pitfalls in the current algorithm is its inability to account for the influence of CR anisotropy on the timing and magnitude of FDs. Relying on a combination of fast Fourier transform and an FD location code, Okike (2020a) provided, for the first time, empirical evidence of significant differences that might exist between FDs identified from unprocessed data and those selected from Fourier transformed CR data. Before making a firm statement on the similar correlation between FDs and IMF at APTY and MOSC stations as well as other relations reported here, the present analysis should be repeated using FD data calculated from Fourier transformed CR flux variation.
Based on strong FDs selected using a large CR(%) ≤ − 3 baseline, Okike et al. (2021b) found FD versus IMF correlation coefficient r ∼ − 0.34 statistically significant at the 90% significance level. The present regression analysis result compares favorably with their finding. This is consistent with the fact that enhanced IMF hinders CR intensity propagation to the Earth and thus results in high magnitude FD detection (Lingri et al. 2016).
The plots of FDMOSC-SWS, FDAPTY-SWS for all simultaneous FD catalog and FDMOSC-SWS for large FDs are displayed in Figure 1(c), Figure 1(d) and Figure 2(b) respectively. Note that their respective R2, r and p-values are 0.17, 0.11 and 0.23, 0.41, 0.34 and 0.48, and 1.27 × 10−10, 1.45 × 10−7 and 2.13 × 10−9. The p-values indicate that the results are statistically significant at the 95% confidence level. The values of R2 suggest that variations (17%), (11%) and (23%) in SWS play some roles that cannot be ignored in comparison with other solar agents that might influence large amplitude FDs.
In an earlier work, Iucci et al. (1979), considering averaged hourly NM data from Alert, Deep River, Goose Bay and Inuvik from 1964–1974, found that FD amplitude is not directly related to SWS. Our results on the two datasets are contrary to this submission. We found evidence of the FD versus SWS relation significant at the 95% confidence level. This disparity could be a pointer to the merits of fully automated FD data over the manual approach employed by Iucci et al. (1979). In separate works, Belov et al. (2014), Badruddin & Kumar (2015), Lingri et al. (2016) and Okike et al. (2021b) reported FD-SWS correlation coefficients of 0.65, 0.58, 0.40 and 0.81 respectively. Our results for the two datasets on average tend to reflect the submissions of Badruddin & Kumar (2015) and Lingri et al. (2016) but are lower than those of Belov et al. (2014) and Okike et al. (2021b). The present result, which is significant at the 95% confidence level, indicates a strong impact of solar wind structure on GCR flux modulations.
FDMOSC-Dst and FDAPTY-Dst for all simultaneous FD data and FDMOSC-Dst for strong FD data are plotted in Figure 1(e), Figure 1(f) and Figure 2(c) respectively. The analysis yields their corresponding R2, r and p-values as 0.16, 0.14 and 0.13, 0.40, 0.37 and 0.36, and 6.51 × 10−10, 6.20 × 10−9 and 1.34 × 10−5. These results indicate that 16%, 14% and 13% of CR flux modulation at MOSC and APTY NM stations can be attributed to Dst index. The chance probability value further affirms that the results are highly significant at the 95% confidence level.
Belov et al. (2001a) performed a regression analysis for 1428 events and obtained a correlation coefficient of < 0.42 between the FD amplitude and the Dst index. The results we find for FD versus Dst parameters for the two datasets are in close agreement with their finding. Recently, Paper I reported a correlation coefficient of 0.46 for the FD-Dst relation based on 129 large CR(%) ≤ − 3 data from the Oulu NM station. While we note that the result obtained from the current data is in close agreement with the result of Paper I, we remark that with respect to Dst, the properties of CR flux modulation appear to be dependent on the point of observation.
5. Summary and Conclusions
Investigation of solar-weather parameters responsible for GCR variations on Earth remains a subject of interest. Several works on this subject employ correlation and regression statistical approaches but often fail to carry out a significance test that accounts for autocorrelation in solar-geophysical data (Belov et al. 2001b; Lingri et al. 2016). The selection of FD as a key parameter is of primary importance in the study of solar-terrestrial connection. Many investigators employ either manual or semi-automated methods in identifying FDs. In this work, we have deployed a fully FD location algorithm to both calculate the amplitude of FDs as well as detecting the number of FDs that occurred between 1996-2005. The FD location code employed here is capable of selecting FD magnitude of −0.01%. The correlation and regression results demonstrate that there were significant space-weather perturbations at the time of the CR flux depressions. All the relations are statistically highly significant. The two-dimensional analysis carried out in this work reveals that that there are other factors than IMF, SWS and Dst that control the amount of CR flux arriving at Earth.
The observed high statistical significance of the correlation between FDs, solar wind data and geophysical parameters could imply that SWS, IMF intensity and Dst have the same causative agent. Their impact on GCR intensity variations is significant. Program selected FDs have shown that solar activity parameters and geomagnetic storm characteristics investigated in this work are important factors that drive variations of FDs (Richardson 2004; Richardson & Cane 2011).
We remain indebted to the teams that host repositories at the websites http://cr0.izmiran.rssi.ru/ and https://omniweb.gsfc.nasa.gov/html/owdata.html from where we freely downloaded the data for this work. In a very special way, we want to acknowledge the non-commercial R software developers. The prompt and enriching input of the anonymous referee is hereby acknowledged. His comments have had a significant effect on the manuscript. Our friends in the R4DS learning community, especially Shamsuddeen who first introduced JAA to the platform, are heartily acknowledged for their invaluable assistance.