GLE # 67 Event on 2 November 2003: An Analysis of the Spectral and Anisotropy Characteristics Using Verified Yield Function and Detrended Neutron Monitor Data

During Solar Cycle 23 16 ground-level enhancement events were registered by the global neutron monitor network. In this work we focus on the period with increased solar activity during late October – early November 2003 producing a sequence of three events, specifically on ground-level enhancement GLE 67 on 2 November 2003. On the basis of an analysis of neutron monitor and space-borne data we derived the spectra and pitch-angle distribution of high-energy solar particles with their dynamical evolution throughout the event. According to our analysis, the best fit of the spectral and angular properties of solar particles was obtained by a modified power-law rigidity spectrum and a double Gaussian, respectively. The derived angular distribution is consistent with the observations where an early count rate increase at Oulu neutron monitor with asymptotic viewing direction in the anti-Sun direction was registered. The quality of the fit and model constraints were assessed by a forward modeling. The event integrated particle fluence was derived using two different methods. The derived results are briefly discussed.

ergies allowing their registration by ground-based detectors e.g. neutron monitors (NMs) (for details see Hatton, 1971;Stoker, Dorman, and Clem, 2000;Dorman, 2006, and the references therein). In this case the energy of SEPs is of the order of about 1 GeV/nucleon or even greater (e.g. Aschwanden, 2012;Reames, 2013;Desai and Giacalone, 2016, and the references therein), which is high enough to produce particle shower in the atmosphere leading to an increase of count rate of ground-based NMs. Accordingly, this specific class of SEPs is known as ground-level enhancements (GLEs) (e.g. Stoker, Dorman, and Clem, 2000;Poluianov et al., 2017).
GLE events can be considered as extreme class of SEP events (e.g. Gopalswamy et al., 2016), with the occurrence rate being roughly ten per solar cycle (Shea and Smart, 1990). We note that, at periods related to maximum and decline phase of the solar activity cycle, a slightly greater occurrence probability was reported (Shea and Smart, 1990;Stoker, 1995;Bazilevskaya, 2005).
GLEs have been continuously studied over the years with various instruments, in order to reveal particle transport in the interplanetary medium, as well as their acceleration (e.g. Debrunner et al., 1988;Ryan, Lockwood, and Debruner, 2000;Aschwanden, 2012;Shea and Smart, 2012;Miroshnichenko, Vashenyuk, and Perez-Peraza, 2013;Klein and Dalla, 2017, and the references therein).
The first registration of a GLE was carried out by ionization chambers (Forbush, 1946). Lately in the mid 1950 -1960 the network of NMs, based on improved detectors, was considerably expanded. Nowadays the worldwide NM network is a convenient multi-instrument tool to study GLEs (Simpson, Fonger, and Treiman, 1953;Simpson, 1957;Forbush, 1958;Hatton and Carmichael, 1964;Simpson, 2000). At present, the worldwide NM networks allows one to perform continuous monitoring of the intensity of cosmic-ray particles and it is used also for registration and study of GLEs and the related space weather phenomena (e.g. Mavromichalaki et al., 2011;Miroshnichenko, 2018).
The intensity of high-energy particles of solar origin is not uniform in the vicinity of Earth, specifically during strong or anisotropic SEP events (Bieber and Evenson, 1995;Bütikofer et al., 2009;Mishev, Kocharov, and Usoskin, 2014). Therefore, it is necessary to deploy NMs at different geographic locations in order to derive their characteristics, because NMs are sensitive to SEPs impinging geomagnetosphere and atmosphere from different asymptotic cones. This fact is specifically important during the event's initial phase when GLE particles reveal an essential anisotropic flux (Debrunner et al., 1988;Shea and Smart, 1990;Vashenyuk et al., 2006;Bütikofer et al., 2009). NM records for all available GLE events are collected in the international GLE database (IGLED) (gle.oulu.fi)  and/or NM data base NMDB (nmdb.eu) (e.g. Mavromichalaki et al., 2011). Here, we employed detrended records of NM count rate increases  retrieved from the IGLED. Detrended data accounts for short-term galactic cosmic-ray variability, such as diurnal variations or recovery phase of Forbush decreases and allows more accurate evaluation of the SEP signal in NM records. Employment of the detrended data is especially significant for obtaining reliable scientific results for weak GLE events during disturbed space weather conditions, in particular for GLE # 67, which occurred during the recovery phase of a Forbush decrease after the previous GLE # 66.
The main characteristics of GLEs differ from each other. The spectra, anisotropy, duration, apparent source location as well as geomagnetic and interplanetary transport conditions and their time evolution vary from event to event Moraal and Mc-Cracken, 2012;Raukunen et al., 2018). This usually leads to a study on a case-by-case basis. For the present study, we focus on the third of the sequence of the three so-called Halloween events, namely the relatively poorly studied 2 November 2003 event (GLE # 67).

Experimental Data of NMs During GLE # 67
During the Solar Cycle 23 16 GLE events were observed (Andriopoulou et al., 2011a,b;. Here we focused on the period with increased solar activity during late October through early November 2003 producing a sequence of three GLEs # 65 -67, known as Halloween events. In terms of spectral and angular characteristics of SEPs, the Halloween events have been extensively studied, except for GLE # 67 Pérez-Peraza et al., 2009).
The GLE # 67 event on 2 November 2003 was related to an X8.3/2B solar flare, which was close to the west limb of the Sun, namely at S14 W56 in the solar active region AR10486. The 2 November 2003 event was associated with a fast CME with speed of ≈ 2660 km s −1 (for details see , and the references therein). On ground level, the event onset was observed between 17:30 and 17:35 UT at several NM stations ( Figure 1). The strongest count rate increases were observed at SOPO (36.0%), MCMD (15.2), % TERA (14.2%) and FSMT (12.0%) NMs, compared to the pre-increase levels. In general, the event was very anisotropic in its initial phase, since no significant count rate increase at SNAE NM was observed. However, the angular distribution of the GLE particles could not be very narrow i.e., beam like, because an early NM count rate increase was observed at several stations whose asymptotic viewing direction were in the anti-Sun direction e.g. APTY, OULU, TXBY (see Figure 3).
A non-negligible trend of the NM data was revealed from the experimental records, therefore detrended experimental data of NM count rates during GLE # 67 were retrieved from IGLED (gle.oulu.fi) . This study was based on data of about 30 NM stations, the full list with the corresponding acronyms, cut-off rigidities and altitudes above sea level is given in Table 1, accordingly the map of the used NM stations is shown in Figure 2.

Mathematical Modeling of the Neutron Monitor Response
Using NM records of GLEs allows one to derive the spectral and angular characteristics of quasi-relativistic and relativistic solar protons by modeling the global NM response (e.g. Shea and Smart, 1982). Routinely it is performed employing Equation 1, that is the relationship between the NM count rates and the primary cosmic-ray flux, and appropriate optimization of a set of n model parameters over m experimental data points. The modeling of the NM network is usually carried out in the energy range of ∼ 0.3 -20 GeV/nucleon. In this study, we employed a method initially developed by Shea and Smart (1982), Cramp et al. (1997), Bombardieri et al. (2006), Vashenyuk et al. (2006), whose detailed description and applications are given elsewhere (Mishev, Kocharov, and Usoskin, 2014;Mishev and Usoskin, 2016a;Mishev, Poluianov, and Usoskin, 2017;Mishev et al., , 2021. For the initial guess of the unfolding procedure, we assumed the apparent source position along the interplanetary magnetic field (IMF) line derived from space-borne measure-  Table 1). ments and/or employing procedure described by Cramp, Humble, and Duldig (1995), explicitly considering the time shift of IMF direction similarly to Mishev and Usoskin (2016b), Kocharov et al. (2017), .
The count rate of a given NM can be modeled using the expression where J (P , t) is the rigidity spectrum of the primary cosmic ray (CR) of galactic or solar origin at given moment t , S(P ) is the specific NM yield function (for details see Clem and Dorman, 2000;Dorman, 2006, and the references therein), G(α(P , t)) accounts for the angular distribution of cosmic-ray particles i.e. the pitch-angle distribution (PAD) of SEPs, we noted that for galactic cosmic ray (GCRs) it is assumed to be isotropic. The A(P) is a discrete function, which accounts the magnetospheric transmissivity (A(P ) = 1 for allowed trajectories, accordingly A(P ) = 0 for forbidden trajectories). A is standardly derived during NM asymptotic cone computations. In Equation 1 P cut is the lower cut-off rigidity of the station. For SEPs, we considered P max = 20 GV, whilst for GCRs P max = ∞. The contribution of GCRs was computed employing the force-field model (for details see Gleeson and Axford, 1968;Burger, Potgieter, and Heber, 2000;Usoskin et al., 2005;Vos and Potgieter, 2015), where the modulation is considered according to Usoskin, Bazilevskaya, and Kovaltsov (2011). Thus, the modeled relative increase of a NM station count rate is given by the ratio between the SEPs and GCR contributions, the latter typically averaged over two full hours before the event's onset (e.g. Usoskin et al., 2015). Herein, for the modeling of NM response, we employed newly computed and experimentally verified specific yield function (Mishev, Usoskin, and Kovaltsov, 2013;Mishev et al., 2020), which is in a conformity with latitude surveys and direct space-borne measurements of GCRs by PAMELA (Payload for Antimatter Matter Exploration and Light-nuclei Astrophysics) (Adriani et al., 2016) and AMS 02 (Alpha Magnetic Spectrometer) (Aguilar et al., 2010) (for details see Mishev, Usoskin, and Kovaltsov, 2013;Gil et al., 2015;Mangeard et al., 2016;Nuntiyakul et al., 2018;Koldobskiy et al., 2019a). This new NM yield function is the most suitable for GLE analysis (e.g. Nuntiyakul et al., 2018;Koldobskiy et al., 2019a;Mishev et al., 2021). Moreover, we modeled the response of the corresponding NMs using a specific yield function relevant to the exact station altitude above sea level (for details see Mishev et al., , 2020Mishev et al., , 2021, which allowed us to reduce uncertainties connected to the employment of the double-attenuation-lengths method, that is, rescaling of high-altitude NM count rates to sea level (e.g. McCracken, Rao, and Shea, 1962), eventually leading to an improvement of the optimization procedure.
Here, for the modeling we assumed various functional expressions of SEP spectra, that is a modified power-law rigidity spectrum similarly to Cramp et al. (1997), Vashenyuk et al. (2006): where J (P ) is the particle flux with as a function of rigidity P in [GV], γ is the powerlaw spectral exponent at rigidity P = 1 GV, accordingly δγ is the rate of the spectrum steepening; or an exponential spectrum: where J (P ) is defined in the same way as in Equation 2, accordingly P 0 is a characteristic proton rigidity.
Here, in all cases the PAD was modeled with a complicated distribution, that is, double (Sun-anti Sun) Gaussian: where α is the pitch angle, σ 1 and σ 2 are parameters governing the width of the distribution, B accounts for the contribution of the particles arriving from the anti-Sun direction. One can see that when B = 0, Equation 4 is simplified to a Gaussian distribution along IMF. Therefore, we can model both: simple PAD of GLE particles (B = 0), as well as complicated angular distributions, including particle flux from the anti-Sun direction. In all cases, it is assumed that GLE particles arrive from the Sun along the axis of symmetry whose direction is defined by angles and (latitude and longitude). Here, the optimization was performed over the set of parameters given in Equations 2-4 over the experimental NM records, that is by minimizing the difference between the modeled and measured NM responses. The unfolding was based on Levenberg-Marquardt method (Levenberg, 1944;Marquardt, 1963) using specific and improved regularization (for details see Aleksandrov, 1971;Golub and Van Loan, 1980;Golub, Hansen, and O'Leary, 1999). The main advantage of our method is the possibility to solve ill-posed problem(s), occasionally arising from marginal experimental NM responses and/or complicated spectral and angular distribution(s), e.g. bi-directional flux of SEPs, (see the discussions in Tikhonov et al., 1995;Mavrodiev, Mishev, and Stamenov, 2004;Aster, Borchers and Thurber, 2005;Mishev, Poluianov, and Usoskin, 2017;. The method was recently verified with PAMELA space probe direct observations and good agreement was reported (for details see Koldobskiy et al., 2019b;Mishev et al., 2021, and the discussion therein).

Results of the Magnetospheric Modeling
As the first step we computed cut-off rigidities and asymptotic directions for all NMs considered in the analysis (Table 1). The magnetospheric transmissivity for the GLE particles was modeled with the MAGNETOCOSMICS code . Here, we used the International Geomagnetic Reference Field (IGRF) geomagnetic model (epoch 2000) as the internal field model (Langel, 1987;Macmillan et al., 2003) and the Tsyganenko 89 model  Table 1. The color lines and acronyms and numbers depict the asymptotic directions and rigidities, which are plotted in the rigidity range 1 -5 GV. The lines of equal pitch angles relative to the anisotropy axis are plotted for 20 • , 40 • , 60 • and 80 • for sunward directions (solid lines), 100 • , 120 • , 140 • and 160 • , for anti-Sun direction (dashed lines), respectively. The cross depicts the measured by ACE space probe IMF direction.
as the external field (Tsyganenko, 1989), which allowed us to obtain the needed rigidity cut-offs and asymptotic directions with a reasonable precision (Kudela and Usoskin, 2004;Kudela, Bučik, and Bobik, 2008;Nevalainen, Usoskin, and Mishev, 2013). Selected NM asymptotic directions used for our analysis are shown in Figure 3. We note that here we presented the asymptotic directions in the rigidity range 1 -5 GV (range of maximal NM response), whilst in the analysis we used the 1 -20 GV rigidity range, as described in Section 3.

Unfolding Procedure Criteria
For the goodness of the fit we employed several criteria: i) The main criterion is related to minimization of the normalized square root of the sum of squared relative differences between the observed and modeled increases for each NM station (residual) (e.g. Himmelblau, 1972;Dennis and Schnabel, 1996): A good convergence of the optimization and a robust solution are achieved in the case when D ≤ 5 -10% for strong and moderately strong events, respectively (Vashenyuk et al., 2006;Usoskin, 2016a, 2018). In case when the count rate increase of NMs is smaller (weak events or late phase of an event) D could be as large as about 15 -20% (e.g. Mishev, Poluianov, and Usoskin, 2017;Mishev et al., , 2021, but still leading to reliable results. As additional criteria we employed: ii) the relative difference between the observed and modeled NM increases must be of the order of about 10 -15% or smaller; iii) uniform distribution of the residuals is required; iv) the value of the reduced χ 2 r should be close to unity, where χ 2 r = χ 2 /DoF, DoF corresponds to degrees of freedom. The combination of these criteria allowed us to select the most relevant solution (e.g. see Dennis and Schnabel, 1996;Tikhonov et al., 1995;. The correspond-ing forward modeling allowed us to assess the quality of the derived spectra and angular distribution similarly to Mishev et al. (2021).

Results of the Analysis
Using Advanced Composition Explorer (ACE) measurements as initial guess for the optimization procedure (step 2) and (step 3) we derived the spectral and angular characteristics of GLE particles. For the analysis we considered 5-min detrended NM data retrieved from IGLED . The analysis was performed assuming different spectral and PAD shapes of GLE particles in order to encompass all the possibilities in the model (Equations 2-4). Additional modeling of the interplanetary transport of SEPs (the details are given in Kocharov et al., 2005Kocharov et al., , 2017, allowed us to adjust the derived best fit by selection of appropriate angular distribution of SEPs. Hence, the derived PADs from the interplanetary transport of SEPs were used as initial guess for the optimization as well as a cross check for the derived from NM data analysis PADs. Thus, using the results from interplanetary transport modeling and the global NM response modeling, we derived a set of spectral and angular characteristics of SEPs, which are self consistent. According to our analysis the most appropriate set of parameters corresponds to modified power-law rigidity spectrum (Equation 2) and double Gaussian PAD (Equation 4) i.e. particles arriving from both sunward and anti-Sun directions. This result is consistent with the NM data where an early count rate increase at Oulu station with asymptotic viewing direction in the anti-Sun direction was observed (see Figure 3). The derived rigidity spectra with the anisotropy characteristics of the high-energy SEPs are presented in Figure 4. The left-hand panels of Figure 4 present the derived rigidity spectra during various stages of the event, namely initial phase (Figure 4a), main phase ( Figure 4c) and late phase (Figure 4e). The corresponding pitch-angle distributions are presented in the right-hand panel of Figure 4, accordingly.
The rigidity spectrum was gradually softening throughout the event, being moderately hard at the event onset and soft at the late phase of the event. The steepening of the spectrum δγ varied throughout the event, resulting in a moderate steepening of the spectrum during the early phase. After 18:10 UT, the derived rigidity spectra can be described by a pure power-law. The particle flux rapidly increased during the initial phase of the event, reached the maximum at about 17:55 UT, and thereafter gradually decreased. The anisotropy was relatively high during the event's onset, yet not very narrow, that is, not a beam like GLE particles flux was revealed. Lately, the derived PAD broadened out. The derived spectral and anisotropy characteristics are summarized in Table 2, where the fit quality criteria D and χ 2 r are also presented.
An additional analysis of the event was performed assuming an exponential shape of the rigidity spectrum (Equation 3). We derived similar anisotropy characteristics (the details are given in Table 3). The derived rigidity spectra were slightly harder compared to a power-law model. The quality of the fit of NM response was similar for both used models. Therefore during the initial phase of the event the NM response can be modeled with either a modified power-law or exponential rigidity spectrum with similar precision. However, after 17:55 UT the modeling with a modified power-law rigidity spectrum was more accurate, but the exponential spectrum still provided a reasonable description. After 18:10 UT, the derived rigidity spectra were described by a pure power-law, while the exponential spectrum modeling leaded to significantly greater difference between modeled and observed NM increases.  Table 2. The black solid line denotes the GCR flux, which corresponds to the time period of the event occurrence. Time (UT) refers to the start of the corresponding five minute interval over which the data are integrated.
We would like to point out that the proposed approach of initial (strong anisotropy and modified power-law or exponential spectrum), main (moderate anisotropy, maximal particle flux and modified power-law spectrum) and late (isotropic flux and pure power-law spectrum) phases of the event differs from that of widely discussed prompt and delayed component proposed by Vashenyuk et al. (2006), Vashenyuk, Balabin, and Gvozdevsky (2011), Miroshnichenko, Vashenyuk, and Perez-Peraza (2013. Yet, the initial and early main phases of the event correspond roughly to the prompt component, whilst the late phase corresponds to the delayed component. Both approaches can be interpreted in relation to different episodes of SEP acceleration (for details see the discussion in Pérez-Peraza et al., 2009;Vashenyuk, Balabin, and Gvozdevsky, 2011;Kocharov et al., 2017Kocharov et al., , 2018.

Quality of the Fit
In order to verify the quality of the fit we performed a forward modeling using the derived spectra and PAD. A comparison between modeled and experimental count rate increases during GLE # 67 for selected NMs is presented in Figure 5. The quality of the fit is similar for the other NM stations. For the comparison we used the best fit, i.e., modified power-law rigidity spectra and double Gaussian PAD (Table 2), while in the other cases the difference between modeled and experimental NM responses is considerably greater.
Accordingly, the contour plot of the sum of variances (Equation 5) for the best fit solutions vs. geographic latitude and longitude is presented in Figure 6. The location of the minimal contour of sum of variances slightly differs from the derived apparent source position (details are given in Table 2). This is due to the application of combination of criteria for the best fit instead of using only one criterion as the minimal residual D. The results of the forward modeling demonstrate that the derived solutions have good precision and quality and the model describes reasonably the experimental data ( Figures 5-6.) The existence, uniqueness and stability of the derived solution(s) were also studied. The stability of the solutions is often violated, because the inverse problems are typically illconditioned and/or ill-posed (e.g. Tikhonov et al., 1995, and the references therein). However, the large number of NM stations with different response (more than 2(n − 1)) and the obtained contour plots of D for the best fit solutions, allowed us to assess the uniqueness and stability of the derived solutions. In case when the initial guess is far from the local minimum, the derived solution possessed significantly greater residual compared to that in Table 2 or it was not physical. Note that ill-posed inverse problems typically converge to  the global minimum only if the initial guess is close to the final solution (Dennis and Schnabel, 1996;Tikhonov et al., 1995). Thus, the subsequent forward modeling using also the information from the derived angular distributions from interplanetary transport modeling, allowed us to derive a reliable consistent solution, in which small variations of the input (initial guess) do not altered the derived set of model parameters. An estimation of model accuracy, namely confidence limits at a level of 95% of J 0 and γ is given in Figure 7. The determination of the spectral and angular characteristics during the late phase of the event was with smaller precision compared to the initial and main phase. This is mostlikely due to isotropization of the SEP flux resulting in more uniform distribution of NM responses, eventually leading to a less accurate assessment of apparent source position.

Particle Fluence
The time, energy and angle integrated particle fluence is presented on Figure 8. Note that for the integration an exponential spectrum was used during the initial stage of the event, whilst modified and pure power-law assumptions were used during the main and late phase of the event, respectively. The fluence during GLE # 67 was compared with recent estimation based on simplified "effective rigidity" approach (Koldobskiy et al., 2019a,b), the details are given in Usoskin et al. (2020). The main simplification of this approach is the assumption that SEP fluxes are isotropic. A very good agreement was achieved, specifically in the region of maximal NM response. Yet, a comparison with earlier reconstruction (e.g. Tylka and Dietrich, 2009;Raukunen et al., 2018) agreed within a factor 2 -3. The observed discrepancy is most likely due to different reconstruction methods and model assumptions, since Tylka and Dietrich (2009) and Raukunen et al. (2018) employed simplified analysis of NM data, different spectral shapes and outdated NM yield function (e.g. Clem and Dorman, 2000). The NM yield function used here (Mishev et al., 2020), is more suitable for GLE analysis (for details see Nuntiyakul et al., 2018;Koldobskiy et al., 2019a), because provides a more realistic assessment of high-altitude NMs responses, that is for AATY, BKSN, CALG, LMKS, MXCO, SNAE, SOPO and TSMB, and was recently verified (Koldobskiy et al., 2019b;Mishev et al., 2021). Modeling of NM response with different specific NM yield   . Panel A depicts the particle fluence, whilst panel B depicts the corresponding confidence intervals at a level of 95%.
functions usually leads to a considerable spread of derived SEP characteristics (Bütikofer and Flückiger, 2015). Therefore, we do not expect exact coincidence, because the simplified modeling in Tylka and Dietrich (2009) and Raukunen et al. (2018), where the anisotropy is neglected and a different yield function was used.

Discussion and Conclusions
In the work presented here, we performed a detailed modeling and reconstructed the spectral and angular characteristics of high-energy SEPs in the vicinity of Earth during the GLE # 67 on 2 November 2003. In this study we used ground-based NM and space-borne data and employed the corresponding data analysis. We examined several possible shapes of GLE particles spectra and angular distribution, namely exponential or modified power-law rigidity spectrum and single or double Gaussian PAD. The best fit of the global NM network revealed a complicated PAD and moderately hard SEP spectra. However, an exponential rigidity spectrum, specifically during the event onset and initial phase, depicted a similar quality of the fit. The event was characterized by a bi-directional particle flux and relatively strong anisotropy during the initial phase. The anisotropy gradually decreased in the course of the event, but remained significant throughout the whole event. According to our analysis we found evidence of particles arriving from both sunward and anti-sunward directions. The dynamical evolution of the spectral and angular properties of GLE particles was obtained in the course of the event.
During the event's onset and initial phase, the rigidity spectrum was moderately hard (γ ≈ 5.5), with moderate steepness δγ ≈ 0.1 -0.2. During this stage of the event the derived PAD was relatively narrow for particles arriving from the Sun direction (σ 2 1 ≈ 1.5 -2.0 rad 2 ), while particles arriving from anti-Sun direction revealed considerably wider PAD (σ 2 2 ≈ 8 rad 2 ), which is close to π . During the main stage of the event, when the peak of the SEP flux was observed, the rigidity spectrum was softer (γ ≈ 6.5), the steepness of the spectrum vanished δγ = 0, i.e., a pure power-law spectrum was derived, and the PAD of SEP arriving from both directions broadened. During the late phase of the event, the PAD was almost isotropic and a constant softening of the rigidity spectrum was derived.
Here, we explicitly assessed the uncertainty and confidence limits of the derived fit of the GLE particles characteristics in a truly way, that is by forward modeling of the global NM response using the derived spectra and PADs. The estimation of the errors in the apparent source position ( and ) determination is about 10 degrees. The confidence boundaries of the derived spectra are given in great detail in Figure 7. One can see that the derived error estimates and confidence limits are within the systematic errors of the employed NM yield function, that is within atmospheric cascade development (e.g. Engel, Heck, and Pierog, 2011;Mishev et al., 2020, and the references therein).
The derived characteristics during the GLEs could serve as a good basis to study their origin (e.g. Kocharov et al., 2017) and/or for space weather applications similarly to Usoskin (2018), Mishev et al. (2021). licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.