January 31, 2026
Direct observation of \(\rm H\alpha\)emission from extended sources in the local Universe often proves difficult due to contamination from geocoronal \(\rm H\alpha\)emission in the sky. Usually, sky spectrum subtraction on small angular-sized objects, such as extragalactic objects or stars, is done by subtracting the pure sky spectrum observed at nearby pointings, which would contain a similar level of geocoronal \(\rm H\alpha\). Unfortunately, when the target is spatially-extended and is larger in angular size than the scale on which geocoronal \(\rm H\alpha\)varies, the sky spectrum observed at a separated sky location would not be sufficient for subtracting geocoronal \(\rm H\alpha\)cleanly.
\(\rm H\alpha\)emission from ionised interstellar medium (ISM) in the Milky Way is one prime example of a spatially extended source impacted by the sky subtraction challenge mentioned above. \(\rm H\alpha\)emission line is widely used as a tracer for the active star-forming (SF) region in galaxies, as young massive stars produce large amounts of ionising photons, ionising surrounding clouds. \(\rm H\alpha\)emission from diffuse ionised gas (DIG) is also important since it contributes to 30 % to 60 % of the total \(\rm H\alpha\)flux we observed from galaxies [1]–[4], which will greatly bias the metallicity measurement of the regions in SF galaxies as DIG introduces significant bias to the observed line flux [5].
In order to study the interaction of regions and surrounding DIG in greater detail, a high spatial resolution of pc scale will be needed to resolve the small-scale structure of the ISM. However, common spatially-resolved spectroscopic surveys on galaxies such as SDSS MaNGA (in \(\sim\) 1-2 kpc scale, [6]) or PHANGS-MUSE (in \(\sim\) 50 pc scale, [7]) cannot provide us with the high spatial resolution required. Surveys like SDSS-V/LVM and AMASE will be able achieve this. Before their data are available, one may be able to use background fibres from surveys like SDSS-IV/MaNGA or MaStar to study DIG. However, we have to find a way to subtract geocoronal emission cleanly.
If the spectral resolution is sufficiently high, the velocity separation (VLSR) between the local standard of rest (LSR) frame and the local rest frame of the observatory can be used to separate the geocoronal \(\rm H\alpha\)(\({\rm \sim 7\;km\;s^{-1}}\) FWHM) from the Galactic \(\rm H\alpha\)(\({\rm \sim 20\;km\;s^{-1}}\) FWHM) [8], [9] by simple double Gaussian fit. Fabry-Perot interferometers have long been used to observe \(\rm H\alpha\)emission from the Galactic ISM, given their high spectral resolution. Some examples include the Wisconsin \(\rm H\alpha\)Mapper Northern Sky Survey (WHAM-NSS) installed at Kitt Peak [8] (R\(\sim\) 25000, \({\rm 12\;km\;s^{-1}}\)), the Wisconsin \(\rm H\alpha\)Mapper Southern Sky Survey (WHAM-SSS) installed at Cerro Tololo [10] (R\(\sim\) 25000, \({\rm 12\;km\;s^{-1}}\) ), and the Dual Etalon Fabry-Perot Optical Spectrometer (DEFPOS) installed at TÜBÍTAK National Observatory (TUG) [11] (R\(\sim\) 11000, \({\rm 27\;km\;s^{-1}}\)). In WHAM, this method worked well because the VLSR is usually offset by \({\rm \sim 18 - 42\;km\;s^{-1}}\) from solar peculiar motion and Earth’s orbital motion, which provides sufficient separation of two emission lines for their spectral resolution. However, the spatial resolution of a typical Fabry-Perot Spectrometer is poor, for example, 1 \({\rm deg^2}\) for WHAM, rendering them unsuitable for the high spatial resolution required to study any internal structure of ISM.
Others have tackled this issue by using the WHAM data as a constraint on their flux calibration [12], [13]. Newer technique independent of WHAM data uses the flux ratio between geocoronal \(\rm H\alpha\)emission and OH \({\rm \lambda 6554}\) sky line emission, and its relation to the sun altitude has been used [14], however, a need for large VLSR \({\rm >\;25\;km\;s^{-1}}\) is used in their fitting given their limited data set consists of both Galactic \(\rm H\alpha\)and geocoronal \(\rm H\alpha\).
We developed a method that can reliably estimate the geocoronal \(\rm H\alpha\)flux from the spectrum based on the shadow altitude with the much larger data set from MaStar, a spectroscopic survey for stars. It does not rely on the large VLSR between geocoronal \(\rm H\alpha\)and Galactic \(\rm H\alpha\) and can work on a relatively low spectral resolution. This method can be used to subtract geocoronal \(\rm H\alpha\)contamination to reveal the underlying Galactic \(\rm H\alpha\).
This paper is structured as follows, we describe the dataset and data preprocessing in Sect. 2, the sky background model in Sect. 3, fitting of the spectrum in Sect. 4, the result in Sect. 5, and discussion in Sect. 6.
MaNGA Stellar Library (MaStar) [15] is part of the fourth generation of the Sloan Digital Sky Survey (SDSS-IV) [6], which uses the Integral Field Unit (IFU) [16] from Mapping Nearby Galaxies at APO (MaNGA) [17] to obtain the spectra of stars with a spectral resolution of R\({\rm \sim 1800}\) during bright time.
The star located at the centre of each IFU can only contaminate a finite angular area, hence the fibres at the periphery of the IFU pointing at the background can be used to study the Galactic ISM, and their spectra will be the data used in this work.
However, directly using these MaStar data to obtain the Galactic spectra comes with a few challenges. Since the observations were done during bright times, significant sky background contamination was present in our data. Although the MaStar data reduction pipeline did handle the sky background subtraction, it was not designed for measuring the Galactic gas emission, as these were not the original targets planned for MaStar. The sky fibres targeting empty sky around the IFUs would contain different levels of the Galactic gas emission and geocoronal \(\rm H\alpha\)to the IFUs due to spatial variation. This is especially true when observing extended targets such as the regions and DIG within the Milky Way. Therefore, we shall construct a model of the sky background first in order to separate the Galactic \(\rm H\alpha\)from the geocoronal \(\rm H\alpha\), the sky background model is explained in detail in Sect. 3.
First, all MaStar visits of all stars are chosen, where a visit is a single night of observation for a star, and a star may have one or more visits. Quality flags tabulated in Table 1 are used to remove any visits that will affect our peripheral fibres, we consider the filtered visits to be good visits in this study. Moreover, small 7-fibre mini-IFUs pointing at standard stars for calibration purposes are not included, given that the starlight could easily contaminate all of the 7 fibres.
| Bit name | Bit number | Description |
|---|---|---|
| HIGHSCAT | 2 | High scattered light in one or more frames. |
| BADFLUX | 3 | Bad flux calibration. |
| POORCAL | 5 | Poor throughput. |
| SEVEREBT | 12 | Severe blowtorch artifact. |
| NONSTELLAR | 16 | Non-stellar spectrum based on visual inspection. |
| CRITICAL | 30 | Critical failure in one or more frames. |
Then, we use the full sky Galactic \(\rm H\alpha\)map from [13]. This map combined observational data from WHAM [8], the Virginia Tech Spectral-Line Survey (VTSS) [18], [19], and the Southern H-Alpha Sky Survey Atlas (SHASSA) [12] together to cover the entire night sky. The WHAM provided a high spectral resolution (\({\rm R \sim 25000}\), \({\rm 12\;km\;s^{-1}}\)) in 1 \({\rm deg^2}\) spatial resolution, while both VTSS and SHASSA are narrow band observations with arc-minute scale spatial resolution (1.6\(\arcmin\) pixel width for VTSS, 0.8\(\arcmin\) pixel width for SHASSA), with their flux calibrated by WHAM observation. We can treat the full sky Galactic \(\rm H\alpha\)map as the ground truth, but with a lower spatial resolution than the MaStar IFU. The full sky Galactic \(\rm H\alpha\)map with the location of each visit filtered for the quality flag in Table 1 over-plotted is shown in Fig. 1.
We utilise this full sky \(\rm H\alpha\)map to separate the MaStar pointings at low or high Galactic \(\rm H\alpha\)locations in the following way. For every MaStar pointing with good visits, we compare its pointing to the Finkbeiner full sky \(\rm H\alpha\)map for the expected Galactic \(\rm H\alpha\)surface brightness (SB), then we separate the IFUs into their corresponding group based on the \(\rm H\alpha\)SB threshold, which is chosen at \({\rm 10^{-17.5}\;erg\;cm^{-2}\;s^{-1}\;arcsec^{-2}}\). This value is near the peak of the distribution of \(\rm H\alpha\)SB (see Fig. 2) while providing enough IFUs for the low Galactic \(\rm H\alpha\)group for statistical studies. Notice this threshold value is arbitrary and can be chosen as any other reasonable value, this will be further discussed in Sect. 6.1.
The distribution of the Galactic \(\rm H\alpha\)SB extracted from the full sky \(\rm H\alpha\)map of each unique pointing is shown in Fig. 2. The low Galactic \(\rm H\alpha\)pointings are assumed to have no or little Galactic \(\rm H\alpha\)SB and are used to calibrate our sky background model, including the geocoronal \(\rm H\alpha\)flux. The high Galactic \(\rm H\alpha\)IFUs can be used to study the Galactic \(\rm H\alpha\)emission for later science purposes.
We use the version of the mgCFrame file with linear wavelength grid (the LIN version) provided by the MaStar data reduction pipeline3. This is a
data product of the pipeline containing the sky-subtracted and flux-calibrated spectrum for each fibre. However, as mentioned, the sky subtraction procedure of the standard pipeline would subtract away both the Galactic line emission and the geocoronal
line emission, we add back the subtracted sky model to each spectrum for the subsequent analysis. There is one mgCFrame file per exposure. Any exposure with an exposure time of less than \(\rm{890\;s}\) and a
transparency of less than 0.8 is excluded in our analysis due to their potential low signal-to-noise ratio (S/N). We also exclude astronomical twilight exposures to prevent significant sunlight scattering and refraction, which is difficult to subtract
correctly.
The wavelength window is limited to be from 6530 Åto 6600 Åfor computational speed, which includes the \(\rm H\alpha\)line, \(\mathrm{[NII] \lambda \lambda 6550,6585}\) doublet lines and eight skylines.
We have in total 10378 exposures with 1146 unique pointings in the low Galactic \(\rm H\alpha\)group and 65373 exposures with 7514 unique pointings in the high Galactic \(\rm H\alpha\)group after all the data preprocessing steps and spectrum fittings described in this section, and Sects. 2.3, 2.4, 3, and 4.
The fibres containing stellar light within each IFU are first excluded since the star is not our target and is a source of contamination. We simply exclude the fibres if the median science flux, \(f_{m}\), of that fibre within a wavelength window without any Galactic emission line (6450 Å- 6500 Å) is larger than a flux threshold.
The flux threshold is determined by sky fibres, which are single fibres targeting empty sky regions around the IFUs. We first calculate the median flux of each sky fibre within the wavelength window in the same way as we do for science IFUs. Then, for each exposure of each plate, we compute the median of these median fluxes among all sky fibres on that plate, \(m\), and their standard deviation, \(s\). We define \(m+2s\) to be the flux threshold to exclude fibres with potential stellar light contamination, which obviously will have larger median flux than that of the sky fibres. To take into account the effect of starlight leaking into nearby fibres, we also exclude any neighbouring fibres of fibres with large median flux. An example of this fibre exclusion step is illustrated in Fig. 3.
IFUs pointing at standard stars for calibration purposes are not included since there are only 7 fibres within these IFUs, and the star would contaminate the entire IFU.
To ensure a high S/N, the spectra within each IFU are stacked together before any fitting. We take the mean spectrum among all the remaining fibres of each IFU for each exposure. Thus, we effectively have a single spectrum with a larger field of view (FOV). Indeed, this stacking will lower our effective spatial resolution. For MaStar IFU, the FOV is from 12\(\arcsec\) to 32\(\arcsec\) (\({\rm \sim 0.12\;pc\;-\;0.31\; pc\;at\;2\;kpc}\)) depending on the size of the IFU, but it still yields higher spatial resolution compared to the Finkbeiner full sky \(\rm H\alpha\)map.
The mgCFrame files contain sky background and must be subtracted away before fitting for Galactic lines at the high Galactic \(\rm H\alpha\)plates. The sky background is modelled in three major components as
shown below,
\[\label{eq:skybg95model} \begin{align} &q(\lambda)=\\ &af(\lambda - k_f) + b + \sum_{i}A_{i, sky}*exp(-\frac{(\lambda - (\lambda_{i, sky} + k_{sky}))^2}{2\sigma_i^2}), \\ \end{align}\tag{1}\] where each term is labelled as follows.
The first term represents those components in the sky background that are due to scattered sunlight by the moon and the atmosphere, where \(a\) is the amplitude, \(f(\lambda - k_f)\) is the shifted solar template, and \(k_f\) is the shifting in wavelength of the solar template.
The second term, \(b\), is the remaining sky background continuum after excluding any components with the solar spectrum.
The third term is the skylines, where \(A_{i, sky}\) are the amplitude of the Gaussian peaks, \(\lambda_{i, sky}\) is the skyline centres, \(k_{sky}\) is the shifting in wavelength of the skylines, \(\sigma_i\) is the post-pixelised instrumental line spread function (LSF) used as the width of the skylines, which varies slightly with wavelength, therefore \(\sigma_i\) located at \(\lambda_{i, sky} + k_{sky}\) is used.
\(q\) is the observed flux.
The details of each component are explained in the following sections.
The majority of sky background contamination comes from the reflection of sunlight from the moon and scattering by the atmosphere, as the MaStar data were observed during bright time. Therefore, we use the observed solar spectrum from Jungfraujoch observatory [20], provided by the BASS 2000 website [21] with a spectral resolution of 0.002 Åfrom 3000 Åto 10000 Åto be the template of our solar spectrum.
Since the template has a much higher spectral resolution and finer sampling than MaStar, we first apply a Gaussian convolution to approximate the instrumental effect, in which the mean of the standard deviation of the pre-pixellised MaStar LSF in the
fitting wavelength window (6530 Å- 6600 Å) is adopted as the standard deviation of the kernel. Then a box kernel of 1 Åwide is used to lower the sampling to match the wavelength grid used in the mgCFrame file. The wavelength window used in
these convolution steps is extended by 10 Åto avoid any boundary effect.
Moreover, the template provided is in air wavelength so it is converted into vacuum wavelength using the equation from VALD3 derived by N. Piskunov [22]. Notice the wavelength grid will not be uniform after this conversion, hence linear interpolation is used to enforce a uniform wavelength grid.
To remove any wavelength calibration issue, the wavelength of the solar template is allowed to shift by \(k_f\), which is limited from -1 Åto 1 Å. At last, the processed solar spectrum template, \(f(\lambda)\), with the same wavelength grid as the mgCFrame file is shown in Fig. 4.
The multiplicative free parameter, solar amplitude, \(a\), is used to represent the strength of the solar template, which could be affected by the illumination fraction of the moon or the angular separation between the moon and the IFU.
The distribution of solar amplitude, \(a\), at different lunar angular separations and lunar illumination fractions is shown in Fig. 5. A clear trend of increasing \(a\) towards the lower right corner can be seen, which corresponds to a small angular separation and large lunar illumination fraction, verifying the idea that the strength of the sky background is correlated to the reflected sunlight from the moon.
The additive free parameter, \(b\), is used to represent any sky background continuum other than the solar continuum. Given the small wavelength window we used, this sky background continuum contribution can be approximated as a flat contribution.
Skylines are extracted from the sky emission spectrum obtained by UVES 4, a high-resolution optical spectrograph of the VLT from the ESO [23]. Not all skylines measured in UVES are used, due to the lower spectral sampling of mgCFrame data. It is noted that the MaStar
90-th percentile FWHM is about 4 Åat \(\rm H\alpha\), therefore, if two neighbouring skylines have their line centre difference < 1 Å, they are first combined as a single line, with their average wavelength to be the new
line centre of the combined line. This is to reduce the number of skylines to be fitted and prevent degeneracy in the spectrum fitting. The fluxes of these two combined lines provided by UVES are simply summed together.
We further removed any skylines, combined or not, with the peak flux smaller than 2 \(\mathrm{\times 10^{-17}~ erg~ cm^{-2}~ \AA^{-1}~ s^{-1}~ arcsec^{-2}}\) in the UVES observation, as there is a group of relatively weak skylines before combining the skylines (see appendix 9). The skylines provided by UVES are also in air wavelength therefore we converted it into vacuum wavelength following the same scheme as in Sect. 3.1.
The skyline components selected are located at 6534.85 Å, 6542.33 Å, 6545.84 Å, 6555.03 Å, 6570.60 Å, and 6579.11 Å, 6585.92 Å, 6598.46 Å, all of which except skyline 6585.92 Åare OH skylines, and the geocoronal \(\rm
H\alpha\)line at 6564.60 Å(wavelength presented in vacuum wavelength). Each line is fitted with a single Gaussian. The centre of the Gaussian component of all skylines shares the same shifting, \(k_{sky}\), with
respect to the original wavelength centre, by 1 Åat most, to address any potential wavelength calibration problem in the mgCFrame file. The line width observed is dominated by the instrumental effect at our spectral resolution, hence we use
the post-pixelised instrumental \(\sigma_i\), as the width of the Gaussian, which varies slightly with wavelength.
The low Galactic \(\rm H\alpha\)IFUs are assumed to have no Galactic components, hence we can use Eq. 1 as the spectrum model, simply fitted by optimize.curvefit package
from scipy. However, the fitting can still be improved by minimising the error introduced in the solar template, and we implemented an iterative refitting process to reduce the error, explained in the following section.
Both the solar template and skylines are allowed to shift freely in wavelength by 1 Åat most in the first fitting as both the solar template and MaStar data may have imperfect wavelength solutions. The median difference, \(\Delta \lambda_{\rm median} = {\rm median}(k_{sky} - k_f)\), between their wavelength shift of all low Galactic \(\rm H\alpha\)IFU pointings is determined to be the systematic difference of the wavelength solution between the solar template and skylines. This median difference, \(\Delta \lambda_{\rm median}\), is added as a constant shifting at the remaining refitting process.
The BASS 2000 solar template is from direct observation of the sun, which could be modified by the moon’s reflection and the scattering of the atmosphere. Part of the fitting residuals due to imperfect template should scale linearly with the solar amplitude, \(a\). Thus, we fit the residual as a linear function of \(a\) with a Huber linear fit at each wavelength grid as shown below,
\[\label{eq:huber95refit} \begin{align} r(\lambda) = m(\lambda)a + c(\lambda), \end{align}\tag{2}\] where \(r(\lambda)\), \(m(\lambda)\), \(c(\lambda)\) are the fitting residual, slope, and y-intercept of the Huber fit at wavelength grid \(\lambda\), respectively.
Huber fit is preferred over the usual linear fit to negate the effect of outliers having a disproportionally large impact on the resulting \(m(\lambda)\) and \(c(\lambda)\) [24]. The threshold used in Huber fit is 1.35.
\(m(\lambda)\) and \(c(\lambda)\) obtained from Eq. 2 are then used to modify the solar template and remaining sky background as shown below,
\[\label{eq:solar95template95corr} \begin{align} f_{\rm corr}(\lambda) &= f(\lambda - k_{\rm sky} + \Delta \lambda_{\rm median}) \\ &+ m(\lambda - k_{\rm sky} + \Delta \lambda_{median}), \\ \end{align}\tag{3}\] where \(f_{\rm corr}(\lambda)\) is the corrected solar template.
\[\label{eq:skybg95corr} \begin{align} b_{\rm corr}(\lambda) &= b + c(\lambda - k_{\rm sky} + \Delta \lambda_{\rm median}), \\ \end{align}\tag{4}\] where \(b_{\rm corr}(\lambda)\) is the corrected sky background continuum.
The solar template, \(f(\lambda - k_f)\), and remaining sky background continuum, \(b\), in Eq. 1 are replaced by \(f_{\rm corr}(\lambda)\) and \(b_{\rm corr}(\lambda)\) respectively. The updated sky background model is then used in refitting to remove any residual, \(r(\lambda)\). This refitting process can significantly reduce the median and standard deviation of the residual at each wavelength after a few iterations. Five iterations are used in this study.
The refitting result is shown in Fig. 6. The median residuals are much closer to zero after refitting, and the root mean square of all medians of residuals at each wavelength decreases from 0.04509 to 0.00195 after refitting. The distribution of the residuals is more symmetric around the median as well.


Figure 6: The median of the fitting residual of all exposures at each wavelength (red cross), the inner 50th-percentile (dark red shade) and the inner 90th-percentile (light red shade) of the residuals of all exposures at each wavelength are plotted. The root mean square of all medians of the residuals of all exposures at each wavelength is shown at the lower right of the figures..
After this refitting process, we then have a final corrected sky background. An example of fitting result using the corrected sky background with Eq. 1 is shown in Fig. 7.
Geocoronal \(\rm H\alpha\)emission originates from the excited hydrogen at the upper thermosphere and exosphere, which is mainly contributed by Lyman \({\rm \beta}\) excitation by solar radiation, and minor contribution from higher energy level (n > 3) cascade process that populates n = 3 (see Sect.8 at [9] for details).
The geocoronal \(\rm H\alpha\)emission rate depends on many factors such as the strength of Lyman \({\rm \beta}\) flux, the profile of hydrogen density in the atmosphere, and viewing geometry. Especially, the geocoronal \(\rm H\alpha\)emission decreases with altitude as the density of hydrogen decreases [25]. The observed geocoronal \(\rm H\alpha\)intensity mostly captures the entire column of geocoronal \(\rm H\alpha\)emission along the line of sight.
The shadow altitude is a geometric viewing parameter [25], defined as the length from the surface of the Earth, in the radial direction, to the intersection between the line of sight (LOS) and the shadow cone of the Earth pointing, shown as the blue arrow in Fig. 8 and blue solid line in Fig. 9. In general, the shadow altitude depends on the location of the observatory, observation time, and celestial coordinate (or equatorial coordinate, which can be easily converted by the observation time), as shown more clearly in Fig. 9. We used a 3D vector calculation detailed in [26] and implemented in [27] to calculate the shadow altitude. We additionally add 102 km to the base radius of the shadow cone since Lyman \(\beta\) radiation would be fully absorbed by atmospheric \({\rm O_2}\) under 102 km [28].
The shadow altitude can be interpreted as the height of the base of the hydrogen emission column, and therefore related to the density of hydrogen at the top of the emission column, with a larger shadow altitude corresponding to a lower density of hydrogen. With this reasoning, we can characterize the observed geocoronal \(\rm H\alpha\)intensity with the shadow altitude of a particular line of sight.
The geocoronal \(\rm H\alpha\)flux observed is defined to be the flux fitted from the sky background model in Eqs. 1 and 3 , \({\rm F(H\alpha)_{obs} \equiv \sqrt{2\pi}\sigma_{H\alpha}A_{H\alpha}}\). The uncertainty, \({\rm \delta F(H\alpha)_{obs}}\), is estimated as the diagonal term of the covariance matrix returned by
scipy.optimize.curve_fit.
We only use \({\rm F(H\alpha)_{obs}}\) of low Galactic \(\rm H\alpha\)pointings, assuming the detected \(\rm H\alpha\)is only coming from geocoronal sources. We further choose pointings with S/N of \({\rm F(H\alpha)_{obs}}\) > 3, which left us with 9202 pointings (88.7% of the original pointings). Data points after the S/N cut are plotted against the shadow altitude, shown below in Fig. 10.
A clear trend of decreasing geocoronal \(\rm H\alpha\)flux with increasing shadow altitude can be seen, this is expected as larger shadow altitude corresponds to a line of sight (LOS) nearly opposite from the Sun, with this LOS passing through the geocoronal layer at very high altitude and low gas density, and therefore less geocoronal \(\rm H\alpha\)emission. This result is consistent with the results found by WHAM [9], [25], [28], [29].
We fit \({\rm F(H\alpha)_{obs}}\) as a function of shadow altitude using a two-segments broken linear model (pwlf, fitting package from [30]) in the log-log scale, inverse variance weighting of \({\rm F(H\alpha)_{obs}}\) is used and the break point is treated as a free parameter. The fitting result is shown in Figs. 11 and 12 and summarised in Table 2. Physically, \({\rm
F(H\alpha)_{obs}}\) never goes to zero because of multiple scattering of Lyman \({\rm \beta}\) in the atmosphere into Earth’s shadow cone [25], hence explaining the gentler slope at higher shadow altitude.
| Segment | 1 | 2 | |
|---|---|---|---|
| Slope | \(-0.652\) | \(-0.334\) | – |
| Intercept | \(4.537\) | \(2.476\) | – |
| Break/End Point | \(5.560\) | \(6.482\) | \(8.693\) |
The same problem of removing geocoronal \(\rm H\alpha\)at a medium spectral resolution was encountered in study based on the observations from LAMOST [14]. They established a linear relationship between the flux ratio of geocoronal \(\rm H\alpha\)to OH\(\lambda\)6554 (6555.03 Åin vacuum) and the solar altitude (be mindful that the solar altitude is a measurement in terms of angles, whereas the shadow altitude is measured as a length).
We repeat their exercise with our data and use the same two segments broken linear fitting as before on high S/N of \({\rm F(H\alpha)_{obs}}\) IFUs detailed in Sect. 5.1. An additional requirement of large VLSR was used in [14], which is not implemented here to provide a fair comparison between both methods.
The result of their method is tabulated in Table 3, and plotted in Figs. 13 and 14. With our much larger data set, we find an overall consistent result with [14], except for small deviation at smaller solar altitudes.
| Segment | 1 | 2 | |
|---|---|---|---|
| Slope | \(-0.0211\) | \(-0.00939\) | – |
| Intercept | \(1.282\) | \(0.877\) | – |
| Break/End Point | \(18.300\) | \(34.616\) | \(80.228\) |
We compare the predictive power of two methods, \(\delta\), which is defined using the root mean square of the fractional residual of all data points, shown below in Eq. 5 ,
\[\label{eq:frac95uncer} \begin{align} \delta &= {\rm rms}\left(\frac{\Delta {\rm F(H\alpha)}}{\rm F(H\alpha)_{fitted}}\right), \end{align}\tag{5}\] where \(\Delta {\rm F(H\alpha) \equiv F(H\alpha)_{obs} - F(H\alpha)_{fitted}}\), and \({\rm F(H\alpha)_{fitted}}\) is the geocoronal \(\rm H\alpha\)flux predicted from the broken linear fitting of each method.
The fractional residual of \({\rm F(H\alpha)_{obs}}\) of each method at different shadow altitudes or solar altitudes is also shown in Figs. 15 and 16. The fractional residuals of each method are fairly stable across different shadow altitudes or solar altitudes. Our method shows a more symmetric and smaller fractional residual distribution compared to the method of [14].
To determine if there is any intrinsic scatter for the models, we have also evaluated the root mean square of the ratio between the fitting residual, \(\Delta {\rm F(H\alpha)}\), and the observed uncertainty of geocoronal \(\rm H\alpha\)flux, \(\delta {\rm F(H\alpha)_{obs}}\), shown below in Eq. 6 ,
\[\label{eq:intrinsic32scatter} \begin{align} {\rm rms} \left(\frac{\Delta {\rm F(H\alpha)}}{\rm \delta F(H\alpha)_{obs}}\right). \end{align}\tag{6}\] This number should be unity if the fitting residual, \(\Delta {\rm F(H\alpha)}\), is limited by the uncertainty of the data points, \(\delta {\rm F(H\alpha)_{obs}}\).
The result is tabulated in Table 4. The error analysis is detailed in the appendix 8.
| Parameter | Shadow Altitude | [14] |
|---|---|---|
| \(\delta\) | \(23.52\%\) | \(49.22\%\) |
| \(\mathrm{rms}\left( \frac{\Delta F(\mathrm{H}\alpha)}{\delta F(\mathrm{H}\alpha)_\mathrm{obs}} \right)\) | \(3.756\) | \(5.231\) |
The root mean square of fractional residual, \(\delta\), from our method is about half of the one based on the method of [14]. Moreover, both method shows a \(\mathrm{rms}(\Delta \mathrm{F(H\alpha)} / \delta \mathrm{F(H\alpha)_{obs}}) > 1\), showing that intrinsic scattering is still present and contributing to \(\delta\).
Despite the success of predicting the geocoronal \(\rm H\alpha\)flux from shadow altitude, the intrinsic scatter is still a major contributing factor to our standard deviation of fractional residual, \(\delta\). We discuss some of the potential contributors to the intrinsic scatter.
First, one may naturally question that, since different LOS correspond to different emission columns sharing the same shadow altitude, shadow altitude is not the only geometric parameter affecting the observed intensity of geocoronal \(\rm H\alpha\). This is also discussed in [28], which suggested a variation of \(\sim\) 1 Rayleigh, and up to 2 Rayleigh at low shadow altitude. We investigate this by plotting the residual of the broken linear fitting against the angle between the shadow altitude line (blue solid line in Fig. 9) and LOS at the shadow cone surface, \({\rm \psi}\), in Fig. 17. However, no significant relation between the residual and \({\rm \psi}\) can be found. This is somewhat expected given the emission column will only change significantly at extremely large \({\rm \psi}\), which is rare in our data.
Second, solar activity could also contribute to the intrinsic scatter. The amount of UV flux from the sun will obviously affect the strength of geocoronal \(\rm H\alpha\), \({\rm F(H\alpha)_{obs}}\), as it originates from Lyman \(\beta\) excitation from the sunlight. We investigated the effect of the 11-year long-term solar cycle and the total solar irradiance (TSI) on \({\rm F(H\alpha)_{obs}}\).
The MaStar data we used span from 11/12/2014 to 30/07/2020, which roughly corresponds to the second half of the solar cycle 24. We split the data set into before and after 1/06/2017, corresponding to the solar maximum period and solar minimum period respectively. We find that there is indeed a stronger geocoronal \(\rm H\alpha\)flux during the solar maximum period compared to that of the solar minimum as expected. A simple K-S test shows that the probability that the two subsamples come from the same parent distribution is almost zero. The result is shown in Fig. 18.


Figure 18: Normalised distribution of fitting residual in log space (panel a). Cumulative distribution of fitting residual in log space (panel b). Blue solid line and red solid line represent data obtained after June 2017 (solar minimum) and before June 2017 (solar maximum), respectively, in both panels. Dashed lines represent the median value of the fitting residual in log space in panel a. The black dashed line represents the location of maximum separation between the two datasets in panel b..
For the TSI, we use the hourly data obtained by the VIRGO Experiment on the cooperative ESA/NASA Mission SoHO (version 8.0, using PMO6V-A+B(Fused) with the new scale [31]. The effect of TSI on the fitted residual of \({\rm F(H\alpha)_{obs}}\) with shadow altitude is shown in Fig. 19. We find an overall increasing trend of the residual of shadow altitude fitting with the TSI. However, the trend is weak when compared to the dispersion of the residual (slope obtained is 0.0519, and \(R^2\) is 0.0394). A two-segment multivariate broken linear fitting of \({\rm log_{10}(F(H\alpha)_{obs})}\) with both \({\rm log_{10}(shadow\;altitude)}\) and TSI together is also tested, \(\delta\) decreases slightly from 23.52% to 22.82%. Hence, we conclude that the TSI alone is not enough to explain all the intrinsic scatter of \({\rm F(H\alpha)_{obs}}\) observed.
Third, the distance between Sun and Earth varies by \(\sim\) 3.3% within a year, with perihelion at early January and aphelion at early July, which causes the solar radiation received by Earth and therefore the observed geocoronal \(\rm H\alpha\)flux, \({\rm F(H\alpha)_{obs}}\), varies as well.
We selected observations during December, January, and February to be the perihelion ones, and observations during June, July, and August to be the aphelion ones. There is a detectable difference between the two datasets, with the perihelion observations having a stronger observed geocoronal \(\rm H\alpha\)flux as expected. A simple K-S test shows that the probability that the two observation sets come from the same parent distribution is almost zero. The result is shown in Fig. 20.


Figure 20: Normalised distribution of fitting residual in log space (panel a). Cumulative distribution of fitting residual in log space (panel b). Blue solid line and red solid line represent data obtained perihelion (i.e., December, January, and February) and aphelion (i.e., June, July, and August), respectively, in both panels. Other symbols are the same as those in Fig. 18..
Another potential contributor to the intrinsic scatter is the imperfect sky background template reconstruction detailed in Sects. 3 and 4.
The majority of the exposures have a lunar reflection of sunlight as the main factor in the modification of the solar spectrum continuum. However, atmospheric scattering and refraction would also modify the solar spectrum continuum. Since these secondary effects would have a minimum impact on the refitting procedure for correcting the solar spectrum continuum (see Eq. 2 ), the reconstruction of the modified solar spectrum continuum would be flawed when these secondary effects become dominant.
Tropospheric scattering could also affect the sky background and strength of geocoronal \(\rm H\alpha\)by about 13% to 18%, as suggested by [28].
These issues have been mitigated by the use of high transparency and pre-dusk selection. However, it is noted that a more sophisticated sky background template can be used in the future.
DIG contamination can also contribute to the intrinsic scatter. We made a simple assumption that the pointings in the low Galactic \(\rm H\alpha\)group (refer to Sect. 2.2) have no Galactic \(\rm H\alpha\)contribution. Given we stacked the spectrum for higher S/N, we have extremely good detection limit on the order of \({\rm 10^{-20}~ erg~ cm^{-2}~ s^{-1}~ arcsec^{-2}}\), depends on the number of fibres in the IFU, therefore, it is inevitable to have DIG emission detected in our stacked spectrum. This DIG contamination is difficult to separate even with a significant VLSR, given the low spectral resolution and freedom of the \(\rm H\alpha\)absorption feature in the solar spectrum template. Fortunately, our fitted results for \(\rm H\alpha\)flux show a typical flux about one order larger than the flux threshold we chose for the low Galactic \(\rm H\alpha\)group, and we do not find any correlation between the Galactic \(\rm H\alpha\)flux provided by the Finkbeiner full sky \(\rm H\alpha\)map and our fitted geocoronal \(\rm H\alpha\)flux. The exact fraction of DIG contamination depends on the shadow altitude, the Galactic altitude, and other parameters discussed above.
In the future, more sophisticated techniques such as PCA can help investigate the sources of intrinsic scatter of the fitting residual between \({\rm F(H\alpha)_{obs}}\) and shadow altitude. Data with minimum sky background contamination and higher spectral resolution can also be used to alleviate the problem of imperfect sky background reconstruction and DIG contamination.
Any intermediate spectral resolution observations around the \(\rm H\alpha\)window that run into the problem of removing geocoronal \(\rm H\alpha\)contamination will benefit greatly from this work, especially when the sky background is significant, or the observed targets are extended on the sky.
A similar approach can be done on higher-order hydrogen recombination lines, albeit the lower S/N from these recombination lines will make this more challenging.
MaStar exposures from the high Galactic \(\rm H\alpha\)group can also be corrected for geocoronal \(\rm H\alpha\)contamination using the calculated shadow altitude and the fitted parameters from Table 2. This will provide us with high spatial resolution and wide wavelength coverage observation of the Milky Way’s interstellar medium, enabling us to study the spatial distribution, line ratios, and even velocity profiles of regions and DIG.
Future high spatial and spectral resolution integral field spectroscopy surveys such as Affordable Multiple Aperture Spectroscopy Explorer (AMASE) [32], [33] and Local Volume Mapper from SDSS-V [34], [35] also targeting the ionised gas in the Milky Way can be used. Geocoronal-\(\rm H\alpha\)-subtracted MaStar observations from the high Galactic \(\rm H\alpha\)group can give us a glimpse of the expected results from AMASE and LVM.
We have demonstrated the ability to remove geocoronal \(\rm H\alpha\)emission from intermediate-resolution spectra from the MaStar background fibres using a simplistic sky background model and the shadow altitude. By selecting pointings with low expected Galactic \(\rm H\alpha\)emission, we are able to characterize the relationship between shadow altitude and geocoronal \(\rm H\alpha\)flux using a two-segment broken linear fit in the log-log space. The fitted result is tabulated in Table 2.
The root mean square of the fractional uncertainty, \(\delta\), a proxy for the predictive power, is 23.52% for our shadow altitude method, which shows a significant improvement when compared to the previous method from [14].
We investigated a few potential sources of intrinsic scatter, including geometric factors from the LOS, solar activity, orbital distance, and an imperfect sky background model. We found tentative evidences that the geocoronal emission depends on the level of solar activity and the distance between the Earth and the Sun. DIG contamination may also be an important contribution to the intrinsic scatter.
Future work can improve the modelling of the sky background and further investigate sources of intrinsic scatter. Overall, this technique provides reliable subtraction of geocoronal \(\rm H\alpha\)and enables Galactic \(\rm H\alpha\)studies at intermediate spectral resolution surveys.
We would like to thank Dr. Sabyasachi Chattopadhyay at SAAO for helpful suggestions on this work. We acknowledge the financial support by grants from the National Natural Science Foundation of China (No. 12425302 and 12373008) and Research Grant Council of Hong Kong (Project No. 14302522 and 14303123). RY acknowledges support by the Hong Kong Global STEM Scholar Scheme (GSP028). The research work described in this paper was conducted in the JC STEM Lab of Astronomical Instrumentation funded by The Hong Kong Jockey Club Charities Trust.
Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS website is www.sdss4.org.
SDSS is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, Center for Astrophysics | Harvard & Smithsonian (CfA), the Chilean Participation Group, the French Participation Group, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU) / University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatório Nacional / MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.
The residual of \({\rm F(H\alpha)_{obs}}\) from \({\rm F(H\alpha)_{fitted}}\) should increase as \({\rm F(H\alpha)_{obs}}\) increases, hence a fractional residual from \({\rm F(H\alpha)_{fitted}}\) is better to determine the predictive power of the fit. Because each data point of \({\rm F(H\alpha)_{obs}}\) yields a fractional residual and the root mean square of all fractional residual of each data point will yield \(\delta\). This value is irrespective of the independent variable, allowing us to compare different methods.
The shadow altitude method evaluates the \({\rm log_{10}(F(H\alpha)_{obs})}\) as a function of \({\rm log_{10}(shadow\;altitude)}\). Since the shadow altitude is calculated from time and celestial coordinates, the uncertainty of the shadow altitude is negligible.
The residual from the fitting is in log space. To convert it into linear space, we can simply assume a Gaussian distribution of residual of \({\rm log_{10}(F(H\alpha)_{obs})}\) centred at \({\rm log_{10}(F(H\alpha)_{fitted})}\) given as follow,
\[\begin{align} \frac{\Delta {\rm F(H\alpha)}}{\rm F(H\alpha)_{fitted}} \approx {\rm log_{10}(\frac{F(H\alpha)_{obs}}{F(H\alpha)_{fitted}})*ln10}, \end{align}\] where \(\Delta {\rm F(H\alpha) = F(H\alpha)_{obs} - F(H\alpha)_{fitted}}\), which is assumed to be much smaller than \({\rm F(H\alpha)_{fitted}}\). \({\rm F(H\alpha)_{fitted}}\) is evaluated at the shadow altitude of the data point using Table 2, \({\rm log_{10}(F(H\alpha)_{obs}/F(H\alpha)_{fitted})}\) is the residual in log space from the shadow altitude method, then \(\delta\) is simply given as,
\[\label{eq:delta95appendix} \begin{align} \delta &= {\rm rms}(\frac{\Delta {\rm F(H\alpha)}}{\rm F(H\alpha)_{fitted}}). \end{align}\tag{7}\]
The Gaussian assumption is only strictly true in linear space. However, since the residual is rather small in our case here, it is justifiable to approximate a Gaussian distribution even in log space. The approximated value of \(\delta\) is still useful in the comparison to the solar altitude method.
The solar altitude method evaluates the flux ratio between observed geocoronal \(\rm H\alpha\)flux and skyline \({\rm OH\lambda 6554}\), \(r_{\rm obs}\equiv {\rm F(H\alpha)_{obs}/F(OH\lambda 6554)_{obs}}\), as a function of the solar altitude. Again, the solar altitude is calculated from the time of observation, hence the residual of solar altitude is negligible.
The residual from the solar altitude method is in a flux ratio, \(r\). To convert this deviation back to linear space of geocoronal \(\rm H\alpha\)flux, we can first find the expected geocoronal \(\rm H\alpha\)flux, \({\rm F(H\alpha)_{fitted}}\), as given below,
\[\begin{align} {\rm F(H\alpha)_{fitted}} = r_{\rm fitted}*{\rm F(OH\lambda 6544)_{obs}}, \end{align}\] where \({r_{\rm fitted}}\) is evaluated at the solar altitude of the data point using Table 3.
Now similarly assuming a Gaussian distribution of uncertainty of \({\rm r_{obs}}\) centered around \({\rm r_{fitted}}\), the fractional residual of geocoronal \(\rm H\alpha\)flux is simply,
\[\label{eq:frac95uncer95Zhang}
\begin{align} \frac{\Delta {\rm F(H\alpha)}}{\rm F(H\alpha)_{fitted}} \approx\; \pm \sqrt{(\frac{\delta {\rm F(OH\lambda 6544)_{obs}}}{\rm F(OH\lambda 6544)_{obs}})^2 + (\frac{\Delta r_{\rm fitted}}{r_{\rm fitted}})^2},
\end{align}\tag{8}\] where \({\rm \delta F(OH\lambda 6544)_{obs}}\) is the uncertainty of \({\rm F(OH\lambda 6544)_{obs}}\) given by scipy.optimize.curve_fit when
fitting using Equations 1 and 3 , \(\Delta r_{\rm fitted} = r_{\rm obs} - r_{\rm fitted}\), and the sign of this equation depends on the sign of
\(\Delta r_{\rm fitted}\). \(\delta\) for the solar altitude method will be given by the same Eq. 7 as above, therefore, the sign of Eq. 8 is irrelevant when getting the root mean square. Although the Gaussian assumption is not strictly true for the flux ratio, the approximation is justified given the small residual.


Figure 21: The UVES skyline flux before (left) and after (right) combining the neighbouring skylines. The UVES skyline flux before combining the neighbouring skylines has a group of relatively weak skylines, as separated by flux level of \({\rm 2\times 10^{-17} erg~ cm^{-2}~ \AA^{-1}~ s^{-1}~ arcsec^{-2}}\) (red dashed line in both figures)..
Email: wkrickywong@link.cuhk.edu.hk↩︎
Email: rbyan@cuhk.edu.hk↩︎
https://www.eso.org/observing/dfo/quality/UVES/pipeline/sky_spectrum.html↩︎