Real‐Time Detection, Location, and Measurement of Geoeffective Stellar Flares From Global Navigation Satellite System Data: New Technique and Case Studies

An alternative approach to detect solar flares and quantify the associated extreme ultraviolet (EUV) solar flux rate was introduced in this journal by the authors: Global Navigation Satellite Systems Solar FLAre Indicator (GSFLAI) was founded on the dependence of the sudden electron content increase of the Earth ionosphere versus the angle regarding the flare source, the Sun, given by a simple but accurate first‐principles‐based model. Such overionization is directly measured from the dual‐frequency phase measurements gathered from hundreds of worldwide permanent receivers of the Global Navigation Satellite Systems, GNSS (like the Global Positioning System GPS), working many of them in real time. In this work we generalize GSFLAI for the very challenging scenario of stellar superflares, with a much weaker expected geoeffectiveness on the Earth ionosphere, making it difficult to distinguish its potential footprint regarding conventional ionospheric variability sources. Indeed, we will show that, unlike GSFLAI for solar flares, the new algorithm presented here (Blind GNSS search of Extraterrestrial EUV Sources [BGEES]) is able to detect EUV flares without the previous knowledge of the position of the source, which is also simultaneously estimated, providing an additional quality check of the detection. It will be first assessed with several case studies of solar flares of different intensities analyzed previously with GSFLAI. Finally, by focusing on the night hemisphere to avoid the Sun's larger effect on the ionosphere, the detection and location with BGEES of two recent stellar superflares, Proxima Centauri (18 March 2016, 08:32UT) and NGTS J121939.5‐355557 (1 February 2016, 04:00UT), are presented, strongly suggesting the extension and applicability of the new technique, also in real time.


Introduction
Stellar flares are sudden electromagnetic emissions on a star's surface that release large amounts of energy. These flares are measured from the Sun by dedicated space missions like the Solar and Heliospheric Observatory (Domingo et al., 1995), and they are directly detected for faraway extrasolar sources by telescopes such as Swift or Fermi placed on Earth orbiting satellites (see, for instance, Troja et al., 2015). However, the sudden increase of the solar flux of photons, mainly in the EUV band but also in the X band, also have an effect on Earth's ionosphere in terms of a sudden electron content increase with an expected spatial pattern in the illuminated part of the Earth ionosphere. This allows, by applying the model and technique introduced in Hernández- , the detection and measurement of the solar flare occurrence and intensity from the dual-frequency trans-ionospheric GNSS receivers measurements gathered from the existing global network of permanent receivers. The technique, GSFLAI, also recently known as Solar EUV flux rate (SOLERA; García Rigo et al., 2018) is accurate enough to detect and measure not only solar flares of high X-ray band high intensity, X-class, but also those of middle-and low-intensity, M-and C-class, respectively (Singh et al., 2015). This is the case specially for those flares not significantly affected by the EUV flux extinction of the solar atmosphere as it was discussed in Hernández-Pajares et al. (2012). GSFLAI has allowed as well to characterize the occurrence statistics of solar flares (Monte-Moreno & Hernández-Pajares, 2014).
In order to extend the applicability of this technique to the detection and measure of stellar flares, the small value of the expected associated electron content increase has made it so far very difficult to distinguish the potential stellar flare footprint from other conventional sources of ionospheric variability (see, for instance, Martínez Cid, 2016). To move forward in this challenging problem, we present a new technique which allows BGEES will be assessed first in several case studies with solar flares with different intensities and geoeffectiveness previously studied in detail. And the feasibility of the detection of such flare events with BGEES for the challenging scenario of faraway stars is studied in two stellar superflares that happened during 2016 in Proxima Centauri and NGTS J121939.5-355557. This technique can open an alternative detection and measurement method of stellar superflares by using free, open, and continuous real-time data streams of GNSS receivers distributed worldwide, like the ones available in the real-time project of the International GNSS Service (see Caissy et al., 2013).

Ionosphere
The ionosphere is a layer of the Earth's atmosphere partially ionized that mostly lies within 75-1,000 km above the surface of the planet (Stanford University, 2019). High energy from EUV and X-ray solar radiation cause its atoms to be partially ionized and create a layer of free electrons (NOAA, 2019). These free electrons are capable of affecting radio wave propagation, thus having an effect on GNSS technology in terms of a advance/delay on the carrier phase/pseudorange signals, proportional to the line-of-sight integrated electron number density (i.e., slant total electron content, STEC) and inversely proportional to the squared carrier frequency. This phenomenon allows these satellites and the network of ground receivers to be used as a global scanner for the ionosphere; see the pioneer work of Afraimovich (2000) and, for instance, Hernández-Pajares et al. (2011) for more details.
The main physical quantity used for describing the electron content of the ionosphere is the vertical total electron content (VTEC). The VTEC is the total number of free electrons in radial direction, referred to the geocenter, along a cylinder of base 1 m 2 .
The unique properties of the ionosphere enable us to use the data provided by the GNSS technology to study solar flares and to try to extend it to the challenging scenario of stellar flares, a powerful phenomenon that occurs in many stars across the universe, which can generate a sudden certain overionization with a similar pattern dependent on the cosine of the source-ionospheric pierce point (IPP) angle (Hernández-Pajares et al., 2012;Singh et al., 2015).
The ionospheric pierce points are the locations that we are going to use for our measurements. An IPP is the point where the line between the satellite and a receiver intersect at the ionosphere's effective height, where the VTEC is considered (Fu & Cazenave, 2000).

Stellar Flares
Flares from stars, in particular those that have the Sun as a source (much more noticeable due to its proximity), are sudden flashes of brightness on the surface of stars which release large amounts of energy across the whole electromagnetic spectrum. This happens in particular in X-rays and specially in EUV radiation, the most geoeffective one in order to generate a sudden increase of ionospheric electron content that can be detected and characterized with dual-frequency GNSS measurements.
Flares that have the Sun as a source can increase the electron content of the ionosphere and have an effect on waves passing through it, as it has been mentioned in section 1.1, affecting the satellite signals in a frequency-dependent way. This phenomenon is a key element which enables us to study these events from multifrequency global GNSS measurements (Hernández-Pajares et al., 2012).
The previous phenomenon applies to flares originating from the Sun. Whether a flare that has a star from outside the solar system as a source has an effect on the Earth's atmosphere (detected using GNSS measurements) or not is one of the challenges the new algorithm faces, as the huge distance might dramatically reduce their effect on the ionosphere; at the same time much higher energy than in the Sun can be released by different astronomical sources, such as Supernova1987A and Proxima Centauri, which might have produced signatures in the Earth ionosphere; see Bohringer et al. (1987), Haisch et al. (1977), and Howard et al. (2018).

The GNSS Solar Flare Activity Indicator
GSFLAI can be defined as the VTEC rate due to an EUV flare at the subsolar point (a + b from equation (2)   model discussed in Wan et al. (2002), is correlated by means of a (another way of defining the GNSS Solar Flare Indicator also presented in Hernández-Pajares et al., 2012) with the rate of direct EUV measurements. These can be obtained from space-based photometers during events such as the solar flares (see also Singh et al., 2015).
The model defining GSFLAI has been the starting point from which the new BGEES algorithm has been developed, augmenting from the inital two unknowns associated to the EUV rate, to four unknowns, after including the additional ones associated to the unknown 2-D angular position of the EUV flare source.

Solar-IPP Angle
The solar-IPP angle or in the case of other stars, source-IPP angle (denoted from now onward; see Figure 1) plays a major role when studying this event: it is the angle formed by the star and the Earth's IPP and drives the effect the flare is having on it. It is expected that the cosine of this variable presents a linear correlation with the increase in VTEC with a factor proportional to the source EUV flux rate (Hernández-Pajares et al., 2012).

VTEC
Additionally, the other parameter relevant to the computation is L I = L 1 − L 2 , the ionospheric combination of carrier phases L 1 and L 2 in length units, a direct measurement of the STEC, biased by a constant term for each phase-continuous transmitter-receiver arc and affected by a small predictable wind-up term (see equation 17 in Hernández-Pajares et al., 2011). Its double time difference within each phase-continuous transmitter-receiver is a simple way of providing a detrended VTEC when, for instance, a sampling time of 30 s is considered, versus the direct usage of VTEC rate when the sampling time is 1 Hz, also suitable for flare detection as it will be detailed in the next sections. See Hernández-Pajares et al. (2012) and the following equation (1).
cos Z is the ionospheric mapping function, that is, the inverse of the cosine of the satellite-zenith angle that we have for each IPP at an effective ionospheric height, at central time t and taking in this work the International GNSS Service adopted value of 450 km (Hernández-Pajares et al., 2009). And S andŜ represent, respectively, the STEC at central time t and the smoothed STEC value, computed from the average of the STEC from previous and next time.
As an example, for the major solar flare (X17.2) that took place on 28 October 2003 (Hernández-Pajares et al., 2012), Figure 2 shows the strong correlation between the cosine of the solar-zenith angle and the VTEC rate in the daylight hemisphere. In this plot, each of the points represents the VTEC rate as a function of the source-zenith angle of the corresponding IPP (computed using the known position of the Sun in this case) for a specific time, in this case at the peak of the flare, on approximately 11 hr 05 min GPS time. We can see the VTEC rate linear increase from approximately cos = −0.1 to cos = 1 ( = 90 • ) as expected by the above mentioned first-principles-based model. It can also be observed that the ionosphere is not experiencing the sudden overionization from cos = −1 (180 • ) to cos = −0.1, consistently with the location of the IPPs in the night hemisphere.

BGEES Algorithm
The correlation between the cosine of the source-zenith angle (for the case of the Sun) and the VTEC rate, which can and is already working in real time, has inspired the BGEES algorithm. It will allow, also in real time, to simultaneously solve the flare intensity and the position of the EUV flare source once the corresponding set of linear equations is solved by using the least squares method.
As we have seen, under any solar flare which has had an effect on the ionosphere (the most of strong, middle-intensity, and weak flares of X-class, M-class, and C-class, respectively; Singh et al., 2015), the VTEC rate in the daylight ionosphere follows a linear relationship on the solar-zenith angle cosine, being the slope proportional to the EUV solar flux rate.
In order to try to extend the applicability of this model to potential overionization due to extrasolar EUV flare sources, we do not know the position of the source which should be estimated as well. The flare source position will be a new reliability parameter of the stellar flare ionospheric signature detection, when compared with the source coordinates directly estimated from direct LEO-or ground-based astronomical observations. We take as starting point the above introduced GSFLAI model, for cos ≥ −0.1 (above the border guaranteeing daylight at the ionospheric heights of few hundreds of kilometers). It expresses the VTEC rate value ( . V) directly obtained from the dual-frequency GNSS measurements, as a linear function of the source-IPP cosine (equation (2)): And the cosine of the source-zenith angle can be expressed as the dot product of both associated unit vectors: where X ′ , Y ′ , and Z ′ are the components of the IPP's unit vector and X, Y , and Z are the unknown components of the source's unit vector.
Finding the value of these unknowns simultaneously, the source's unit vector (X, Y , Z) and the coefficients a, b associated to the EUV flare intensity value, is the challenge of this method. Introducing equation (3) into equation (2), we can express the model as In principle the resulting equation (4) is not linear in terms of the unknowns a, b, X, Y , and Z, with X, Y , and Z being tied to be the components of the source position unit vector (i.e., √ X 2 + Y 2 + Z 2 = 1). But we can group them as follows: Being = aX, = aY and = aZ. In this way the resulting equation (5) dependence becomes linear, facilitating very significantly in this way its future applicability also in real-time conditions. After solving initially the previous set of simultaneous linear equations (for instance, in the region of the night ionosphere V values) we can derive the first estimation of a, b, X, Y , and Z for the potential source in case the intensity parameters, specially a, is large enough. This can be done in the following way: knowing that X, Y , and Z are the components of a unit vector, the potential source position. This allows us to, once we know the values of , , and , obtain X, Y , and Z, with an initial sign ambiguity given by the initial sign ambiguity of |a|, by doing and similarly for obtaining Y and Z. The ambiguity can be fixed, for instance, by looking at the predominant sign of . V at the point with the ionospheric point with the stellar source at the zenith (hereinafter subsource ionospheric point, SSP) which should typically provide the corresponding sign of a (predominant in the a + b . V term at SPP) and, hence, the unambiguous values of X, Y , and Z. Moreover, the value of a above a certain positive threshold, greater than the predominant background values, can be used as flare detector, specially at high rate (e.g., 1 Hz) due to the implicit VTEC detrending as it has been commented in section 2.2.
And from these values, which have been obtained deriving the IPP coordinates X ′ , Y ′ , Z ′ from the corresponding right ascension and declination values, ′ and ′ , that is, in a conventional astronomical reference frame, we can derive the first estimated position of the potential source as where s , s are the source declination and right ascension, respectively.
Furthermore, another optimization considered is executing the algorithm multiple iterations. Using the estimated location yielded by each iteration, the algorithm computes the cosine of the angle between the estimated solution and the IPP and discards them below a certain source-IPP cosine threshold: −0.1 in which, as it has been previously shown, the linear model would not apply due to the eclipse of the Earth on the source radiation affecting the ionosphere.
Finally, a refined approach, which significantly reduces the error of the BGEES technique when there was a reduced number of available permanent GNSS receivers or during stellar flares or weak solar flares, can be achieved by combining the measurements of two epochs with a similarly high VTEC rate. Then, by assuming the same source, a common initial model provided in equation (2) can be considered, and the full BGEES approach can be still considered valid but with a better geometry.

Results
Before applying the BGEES algorithm with flares from faraway stars we have analyzed its performance over several case studies of solar flares of different intensities, times, and geoeffectiveness. They have been analyzed in detail in previous publications, to study how it behaves and what model aspects could be tweaked to improve the performence.

BGEES Applied to Solar Flares
A set of nine solar flares, with different intensities, geoeffectiveness, and values relative to the peak, have been selected to perform a first assessment of BGEES, also in challenging conditions (see first three columns in Table 1). They were previously studied in detail by the classical solar-centered technique GSFLAI in Hernández-Pajares et al. (2012) and Singh et al. (2015). The measurements provided by the same or similar distributed set of permanent ground GPS receivers used in such papers, shown for instance in Figure 2 of Hernández-Pajares et al. (2012), have been processed with BGEES for studying the solar flares and stellar superflares.
It can be seen in Figure 3 that, as it might be expected, the final error of the source location, in degrees, is strongly dependent, decreasing on the final ionospheric geoeffectiveness of the solar flare, expressed as the average of the double time difference of VTEC on the daylight hemisphere.
This can be seen in detail in Table 1, where we have approximately ordered the set of representative of strong (X-class), middle-intensity (M-class), and weak (C-class) solar flares from the less to the most challenging cases. The estimated position error of the Sun at the given times with BGEES is included in the fifth column of the table, being less than 10 • for the most part of strong and middle-intensity flares, also when they are analyzed at 10% of the peak intensity. The main exception is the X-class flare with ID 2002.196.72240 in Table 1, with an error of 47 • , coinciding with the very low number of worldwide available permanent GPS receivers at 1 Hz of measurement rate (fourth column). In order to improve this and other results, two epochs of the same solar flare with similar distribution of double time difference of VTEC have been combined in these case studies, as it has been described at the end of the previous section. The source positioning error (sixth column in the above mentioned table) improves in general, specially in such a case: see values in bold and italic fonts for the smaller and second smaller value below 10 • of error (if applicable). This is consistent with the larger number of IPPs when we combine two epochs, specially necessary for such solar flares happened 18 years ago, that is, coinciding with a less extensive ground network of GPS receivers. The situation is completely different nowadays, not only for the existence of a denser global network of permanent ground receivers but also due to the higher number of IPPs per time with four operative GNSS constellations (GPS, GLONASS, Galileo, and Beidou), then multiplying approximately up to more than 3 the number of available IPPs per multiGNSS receiver.

First Results of BGEES With Two Extrasolar Sources
One of the first considerations to be implemented in order to study stellar flares was to avoid the potential interferences of the Sun EUV flux variability in the BGEES functioning. Because the Sun would have a greater effect on the Earth's ionosphere, the GNSS observations within the daylight hemisphere are discarded. Like in optical astronomy this leaves the algorithm with up to approximately one hemisphere, the night one, useful, which could cause potential flares to be missed if they were having an effect on the day hemisphere, that is, if the extrasolar source is not far away from the Sun. This is done by simply taking the IPP's and Sun's right ascension and declination and discarding the IPPs with a cosine smaller than −0.2 (the point at which a linear correlation between the variables can be observed, as can be seen in Figure 2).

Proxima Centauri
The first studied extrasolar flare was presented in Howard et al. (2018). The event corresponds to the first naked-eye-brightness superflare detected from Proxima Centauri, the closer star to the Sun. It was detected on 18 March 2016 with the flare peaking at 8:32 UT. Because this was a more challenging scenario, instead of finding the moment of the peak and only focusing on that time, the algorithm search, over the set of more likely epochs surrounding the moment of the peak, was based on GPS data only.
Taking into consideration that Proxima Centauri has a right ascension of 217.4294 • and a declination of −62.67948 • (WikiSky), at approximately 110 • from the Sun during this event, the algorithm was executed using only the night hemisphere to try to detect this flare from 250 worldwide distributed GPS receivers, using such location as reference to calculate the error of the source position blind estimation with BGEES.
In the previous section the option of using the two best epochs, that is, with higher and similar VTEC double time differences VTEC, was introduced, which yielded improved results for some of the assessments done  Moreover, as described in the last paragraph of section 3, it is convenient to perform an additional iteration to subtract the IPPs not illuminated by the extrasolar source by considering the first estimated location, that is, not following the assumed overionization model (equation (5)).
As we can see (Figure 4), in the first iteration we obtain location errors greater than 20 • and include a time close to the reported peak time of 8:32 but with a location error of around 50 • for this time. The results improve significantly in the second iteration, which provides a source location solution with only 4.2 • of error, when combining the two similarly overionized epochs, 8.5083 and 8.3 hr (the two inverse points seen in Figure 5), very close in time, especially the first, at the peak time of the reported superflare in Howard et al. (2018).
In order to double check the consistency in the detection of the footprint of this superflare of Proxima Centauri on the Earth's ionosphere, we have computed the time evolution of the mean value of VTEC rate at 1 Hz at the SSP, that is, around the coordinates of the source with a threshold of cos > 0.97 (i.e., ≤ 14 •  approximately), once the time and location of the event has been found with BGEES. In spite of the high noise, as expected, it can be seen in Figure 6 that the VTEC rate becomes mostly positive few minutes before 08:32, oscillates right after, and becomes negative; that is, VTEC decreases up to more than 40 min later. This trend is compatible with the light curve shown in Figure 1 of Howard et al. (2018).

NGTS J121939.5-355557
The second stellar superflare studied with BGEES was presented in Jackman et al. (2018) and corresponds to the star NGTS J121939.5-355557, with approximate right ascension and declination of 184.91 • and −35.93 • approximately. The event, which took place on 1 February 2016 with an angular Sun-star distance of 108 • approximately, was detected just less than half an hour after the start of the night. The light curve peak at 04:00 UT (personal communication of the authors), with pulsations around during approximately 30 min (see top plot in Figure 1 of Jackman et al., 2018). In this context we have analyzed around 1 Hr of 1 Hz GPS data surrounding the peak, from 142 worldwide GPS receivers and comprising such time interval of pulsations.
It can be seen that BGEES, during the first and specially the second iteration with the refinement of the ionospheric region (Figures 7 and 8, respectively), provides minima of the NGTS J121939.5-355557 location error below 20 • approximately, during the times of the pulsating stellar flare events, especially around 04:03 with a source location error below 10 • .
Similarly as we have done in the superflare of Proxima Centauri, we have computed the time evolution of the mean value of VTEC rate at 1 Hz at the SSP, once the time and location of the event has been found with BGEES. Figure 9 shows high-frequency VTEC oscillations during approximately half an hour around the peak at 04:00, consistently with the time occurrence of the major pulsations in the stellar flare, as it is shown in the top plot in Figure 1 of Jackman et al. (2018).

Conclusions
The main aim of this work has been developing an algorithm capable of detecting EUV flares without considering the location of the source, first exemplified in the Sun, that is, both in time and space, in order to be applied to the challenging scenario of stellar flares.
The proposed BGEES algorithm has been introduced and tested first with several solar flares, which provided positive results regarding the detection and source location in several solar flares of different intensities and geoeffectiveness, previously studied in detail with the former technique.
From this positive assessment with solar flares, BGEES has been first validated as well with two huge stellar flares detected and characterized with conventional astronomical observations. In the analysis of both flares, the one from Proxima Centauri and the one from NGTS J121939.5-355557, BGEES presents a significant reduction in the source positioning error at the moment of the flare, below 10 • . This can be considered the fulfilment of an important extra condition. This, jointly with the analysis of the time evolution of the Space Weather 10.1029/2020SW002441 high-frequency VTEC rate below the given star consistent with the directly measured light curves, supports the conclusion that the indirect new GNSS-based technique, successfully working with solar flares, has been able to be sensitive as well to these two stellar superflares.
In conclusion, the aim of this work was to develop an algorithm to be able to detect EUV flares, both in time and source location domains. We consider it has provided encouraging results that make us aim to keep working on this direction, studying more recent huge events where we can take profit of the higher number of ionospheric pierce points due to the larger number of GNSS receivers and transmitters, with up to four constellations (GPS, GLONASS, Galileo, and Beidou). In this way we hope to be able to consolidate and extend the number of stellar flare detections with GNSS, which might be helpful to the scientific community in the future, also in terms of a real-time warning service. The last but not the least, the authors consider that this work can contribute to open a new field, the GNSS astronomy, which might be helpful in particular in the study of the habitable zone of exoplanets (see, for instance, Howard et al., 2018).