Blueshifted Hydrogen Emission and Shock Wave of RR Lyrae Variables in SDSS and LAMOST

RR Lyrae stars are Population II stars commonly found in globular clusters. They are located at the intersection of the instability strip and the horizontal branch (HB), from A to F type (Smith 1995; Sesar 2012). RR Lyrae variables are classified, according to the number of oscillation modes, as fundamental mode (RRab), the first-overtone mode (RRc), or double-mode (RRd) stars (Soszyński et al. 2011). They obey the M_V –[Fe/H] relation (Muraveva et al. 2018) and the period–luminosity–metallicity relationship (Longmore et al. 1986; Catelan et al. 2004) in the infrared passbands, which makes them good “standard candles” for precise distance determinations for star clusters and nearby galaxies (Bhardwaj 2020). Some of the RR Lyrae variables, across all subtypes, show periodic amplitude and/or phase modulations, which is known as the “Blazhko effect” (Blažko 1907; Kolenberg et al. 2006).

Despite the successful application of RR Lyrae stars as a step of the distance ladder, the detailed dynamical evolutions of the interior are unclear and complicated. However, we can speculate on the internal structure of RR Lyrae variables by tracing the imprints of feature lines on spectra. Hydrogen emission lines, helium emission lines, line broadening and doubling phenomena, and neutral metallic line disappearance phenomena in RR Lyrae stars indicate shock waves propagating through their atmospheres (Chadid et al. 2008, 2017). The emissions are produced when atoms in the gas de-excite after being excited by a shock wave traveling into the atmosphere. The shock provides high enough energy to excite neutral hydrogen from the second quantum state upward. It compresses, heats, and accelerates the gas. After that, the atoms release energy in the cooler shock wake region and generate the emission (Gillet & Fokin 2014).

Hydrogen emission lines of RR Lyrae stars are called “apparitions” (Preston 2011). There are three famous “apparitions,” which got their sequence number based on the discovery time. Sanford (1949) first reported a hydrogen emission in RR Lyrae in 1949. The “first apparition” is an intense blueshifted emission, which mostly presents during ϕ ≈ 0.91 (in RRab), just before the maximum luminosity (Gillet et al. 2019). Duan et al. (2021) reported the first detection of blueshifted hydrogen emission (the “first apparition” at both Hα and Hβ) in the first-overtone and double-mode RR Lyrae variables. Now it is well accepted that the “first apparition” is generated from the hot emitting layer behind an outward-moving shock front (Schwarzschild 1952).

The “second apparition” is a weak hump at the blue side of the absorption-line wing, when ϕ ≈ 0.73, discovered in 1988 (Gillet & Crowe 1988). It is generated by the collision when layers of the upper atmosphere meet the photospheric layers during the infalling processes of the ballistic motion. It occurs during the small bump around ϕ ≈ 0.72, which was produced by the “early shock” (Hill 1972; Gillet & Fokin 2014). Moreover, the third one was discovered in 2011 (Preston 2011), which is a weak redshifted emission line that appears at around ϕ ≈ 0.30. There are two views on the origin of the “third apparition.” One is that it is possibly provided by a weakly supersonic and post-maximum infalling shock wave (called Sh_PM), which is generated from superimposing compression due to the hydrogen recombination front and the accumulation of weak compression waves (Chadid & Preston 2013). The other opinion is that this emission is a P-Cygni−like profile, which is the sign that the envelope surrounding the stellar surface is expanding, and the emission was produced when the shock wave is detached from the photosphere. To confirm which one is closer to the fact, further investigation is still needed.

Thanks to the rapid improvement of observing facilities, we can get an ever-expanding volume of observed spectroscopic data of RR Lyrae variables. In this work, we develop a pattern recognition searching algorithm and apply it to a large database of single-epoch spectra, the “first apparition.” We find spectra of RR Lyrae stars that show the most prominent observational characteristics of shock waves. We report hydrogen emissions in 103 RRab, 20 RRc, 3 RRd, and 1 Blazhko-type RR Lyrae variable. ⁷ The detection of blueshifted hydrogen emission in nonfundamental mode RR Lyrae variables has been discussed in our previous work (Duan et al. 2021). In this work, we present the results of RRab variables and the Blazhko-type RR Lyrae stars (RR-Bl, the Blazhko-type RR Lyrae star with Large Sky Area Multi-Object Fiber Spectroscopic Telescope survey (LAMOST) spectra and hydrogen emission). The Blazhko-type RR Lyrae star is characterized as an RRab pulsator using the pre-whitening sequence method. We build up the largest database of properties of the “first apparitions” in RR Lyrae stars based on the Sloan Digital Sky Survey (SDSS) and LAMOST, to investigate shock waves in RR Lyrae stars with a new view.

This paper is organized as follows. In Section 2, we describe the observations of both photometry and spectroscopy. The searching methods are interpreted in Section 3. We display the results of the search in Section 4. Moreover, we analyze the frequency components in RR-Bl, characterize it as an RRab star, and provide a possible Blazhko period in Section 5. The detection rate and mock spectra test are discussed in Section 6. Finally, the conclusions are covered in Section 7.

2.1. Photometric Observations

2.2. Spectroscopic Observations

Our spectroscopic data were collected from SDSS and LAMOST. In this work, we adopt low-resolution and single-epoch spectra. The periods of RR Lyrae stars are very short. Under the circumstances, the emission features can be overwhelmed by long-time exposures and the processes of combining spectra to get a higher signal-to-noise ratio (S/N). Therefore, we use low-resolution and single-epoch spectra to capture these features, which will quickly wear away otherwise.

SDSS (Eisenstein et al. 2011) is one of the most ambitious surveys in the history of astronomy and enjoys enormous influence. It saw its first light in 1998 and entered routine observations in 2000. Using the Sloan Foundation 2.5 m optical telescope (Apache Point Observatory, New Mexico), SDSS serves images and spectra of the northern sky, and it currently extends to the southern sky using the 2.5 m Du Pont optical telescope (Las Campanas Observatory, Chile). SDSS-III is a massive spectroscopic survey focusing on the distant universe, the Milky Way, and extrasolar planetary systems (Eisenstein et al. 2011). For the SDSS data set, some of the spectra come from the Baryon Oscillation Spectroscopic Survey (BOSS; Dawson et al. 2013). Spectra from SDSS cover wavelengths from 3800 to 9200 Å, and those from BOSS have wavelength coverage of 3650–10400 Å. They achieve R ∼ 1500 at 3800 Å and R ∼ 2500 at 9000 Å. In this work, we use SDSS DR9 because we utilize parameters from the SEGUE Stellar Parameter Pipeline (SSPP; Lee et al. 2008).

LAMOST (Cui et al. 2012; Deng et al. 2012; Zhao et al. 2012) is a Chinese national scientific research facility. This survey consists of two major parts: the LAMOST ExtraGAlactic Survey (LEGAS) and the LAMOST Experiment for Galactic Understanding and Exploration (LEGUE). It sets up a spectroscopic survey of millions of objects in much of the northern sky and has great potential to survey a large volume of space efficiently. We apply for low-resolution and single-exposure spectra in LAMOST DR6 (R ∼ 1800). Each plate was observed more than 3 times for SDSS and 1–3 times for LAMOST. For SDSS, the exposure time is 10–40 minutes each; for LAMOST, it is 10–30 minutes each. This provides the finest time resolution for RR Lyrae targets.

From this large database, we can screen out spectra of RR Lyrae stars. We cross-match the catalog of RR Lyrae stars and the data set of spectra with TOPCAT (Taylor 2005) and apply for single-epoch spectra from LAMOST. We get single-epoch spectra of 3526 RR Lyrae stars from SDSS DR9 and 5571 RR Lyrae stars from LAMOST DR1-6. The locations of RR Lyrae stars of this sample are shown in Figure 1.

Zoom In Zoom Out Reset image size

Based on the profiles of the “first apparitions,” we adapt two methods to search for this kind of feature. In this section, we will explain the technical details of the pipeline we use to catch the blueshifted hydrogen emission lines of RR Lyrae stars in a large sky survey database.

3.1. 1D Pattern Recognition Method

3.2. Cross-correlation Method

Given that the appearance of the emission lines introduces a sudden change on the spectra, the cross-correlation method is also useful for selecting blueshifted hydrogen emissions (Yang et al. 2014). If there is not a “first apparition,” there should only be one wide absorption line of Hα on the spectrum. On the contrary, there would be one additional emission line, which introduces a big difference in the shapes of two spectra. So if the difference between radial velocities is strangely large, there are possibly the “first apparitions.”

The cross-correlation method implemented here works as follows: We use the first observed spectrum in a pair as the template, while the second one is shifted by 1 km s⁻¹ steps between −250 and 250 km s⁻¹. The shifted template is linearly interpolated at the points of the wavelength of the observational spectrum to calculate the cross-correlation function. This shows the difference of radial velocities (ΔRV) between two spectra when the cross-correlation function reaches its maximum. Here we provide the results of differences of radial velocities with the Hα line in Figure 2.

Zoom In Zoom Out Reset image size

We locally normalize the spectra between 6540 Å < λ < 6590 Å by the fitting result of the continuum. The quantity σ of the spectra is calculated from the rms of the spectra from the first and last five points, within which the spectra are continuum dominated. We adopted results with both spectra satisfying S/N > 15, and the significance of the absorption line is over 3σ. The distributions of ΔRV for both the SDSS and LAMOST data sets are visualized in Figure 2. We check for the blueshifted hydrogen emission when ∣ΔRV∣ > 150 km s⁻¹.

All of the selected results are produced by the 1D pattern recognition pipeline. The cross-correlation method provides 13 RRab stars and 1 RRc star for SDSS and 44 RRab stars, 3 RRc stars, 1 RRd star, and 1 Blazhko-type RR Lyrae star for LAMOST. The distributions of ΔRV for selected RR Lyrae stars are visualized in Figure 2. The cross-correlation method does help, but not all of the spectra with emission show high ΔRV, and not all spectra with high ΔRV show emission. So our work is mainly based on the 1D pattern recognition pipeline. The cross-correlation method is used as a supplement.

As a result, we find 33 RRab stars and 10 RRc stars in SDSS and 70 RRab stars, 10 RRc stars, 3 RRd stars, and 1 Blazhko-type RR Lyrae star in LAMOST that show the “first apparition.” Visualization of the spatial distribution of the selected stars is shown in Figure 1. The detailed analysis of the results of RRc and RRd stars has been demonstrated in Duan et al. (2021). In this paper, we present the measurements of RRab and Blazhko-type RR Lyrae variables and compare the properties between different types of RR Lyrae stars.

The emission and absorption components are fitted with the scale width versus shape method for the hydrogen lines (Sersic 1968; Clewley et al. 2002). We use two functions (Xue et al. 2008; Yang et al. 2014; Duan et al. 2020) as

to fit the profile. Uncertainties are provided by error propagation with the covariance matrix and Monte Carlo method (Andrae 2010).

Figures 3 and 4 show the fitting examples of selected RRab stars in SDSS and LAMOST, respectively, which display the blueshifted emission in Hα. Figure 5 shows the fitting result of the Blazhko-type RR Lyrae star with the “first apparition” detected. The wavelength axis is in the stellar rest frame. Flux is normalized by continuum. The Hα emissions are shown as pink profiles. Hβ profiles are shown in the subplots. The significances of the blueshifted emissions for both Hα and Hβ are also shown. The candidates are adopted when the signal of Hα is over 2σ and the signal of Hβ is over 1σ.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The radial velocity of the blueshifted hydrogen emission in the stellar rest frame is calculated as

where λ_e1,α denotes the wavelength corresponding to the central wavelength of the emission line, λ_ab denotes the central wavelength of the absorption line, and λ₀ denotes the laboratory wavelength.

We provide measurements of redshift z_e1,α and radial velocity in the stellar rest frame V_e1,α , normalized flux Flux_e1,α , and FWHM of the emission and absorption FWHM_e1,α in the stellar rest frame of the blueshifted hydrogen emission. Parameters of the selected RR Lyrae variables and measurements of the “first apparitions” are summarized in Tables 1 and 2. RR Lyrae stars with SDSS spectra that exhibit the first apparition are marked as RRabs (s for SDSS) and those with LAMOST spectra and hydrogen emission as RRabl (l for LAMOST). The star “RR-Bl” denotes the Blazhko-type RR Lyrae star with LAMOST spectra and hydrogen emission (l for LAMOST). It was temporarily classified as a “Blazhko-type RR Lyrae variable” because it was not accurately classified owing to the significant modulation in the Catalina Sky Survey. In this work, RR-Bl is characterized as an RRab star using the pre-whitening sequence method in Section 5.

The distribution of selected stars on the period–amplitude diagram is shown in Figure 6. A clear gap divides RRab and RRc stars into two groups. RRab stars have longer periods and larger mean amplitudes than RRc stars. Box plots of the measurements are displayed in Figure 7. The mean value of normalized flux and the FWHM of emission or absorption of RRc stars are higher than those of RRab stars. The mean value of V_e1,α of RRab stars is higher than that of RRc stars.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 8 shows the log T_eff–log g diagram. We adopt measurements of T_eff and log g generated by SSPP if available. Despite that the pipeline is not optimized for RR Lyrae stars and these values are derived from co-added spectra, the parameters from the pipeline can provide an overall description of log g–log T_eff distribution of various subgroups. Overall, RRc stars are hotter and have larger log g than RRab stars. We also display simulated HB evolutionary tracks for stars in the background, which was calculated with expansion and semiconvection of the core and enhanced oxygen composition from Dorman (1992), to show regions of different masses. The line consisting of red, green, and blue points in the background indicates evolutionary tracks for stars with M = 0.54, 0.58, and 0.68 M_⊙, respectively.

Zoom In Zoom Out Reset image size

The presence of the “first apparition” indicates a certain range of phase, just before the maximum luminosity. Figure 9 shows that points during a certain phase interval from RR Lyrae variables with similar properties are concentrating on a linear path on the log T_eff–log g diagram when the effective temperature changes, which is demonstrated by a series of simulations using Modules for Experiments in Stellar Astrophysics (MESA; Paxton et al. 2019). We simulated nine RRab and nine RRc stars with various initial effective temperatures in identical duration of evolution, while other initial parameters remain the same. The initial parameters are M_init = 0.55 M_⊙, L_init = 41.687 L_⊙, X_init = 0.736, and Z_init = 0.008, which indicates [Fe/H] = −0.33 (Das et al. 2018). T_eff ranges from 6300 to 7100 K in step width of 100 K. In Figure 9, variation of color indicates the changing of luminosity, which can stand for the changing phase. Points at around ϕ = 0.91 are located on a linear path, highly characterized by Pearson correlation coefficients (PCCs) of 1.000 and 0.999 for RRab and RRc, respectively. The observation points on the log T_eff–log g diagram have a diffuse distribution because our selected stars have various physical properties such as mass, luminosity, and metallicity. The different parameters make them deviate from one linear path.

Zoom In Zoom Out Reset image size

We analyze the r-band light curve of RR-Bl with a time span of 974.66 days with a standard successive pre-whitening method (Moskalik & Kołaczkowski 2009). Observations in poor quality (catflags ≠ 0) have been removed. We fit the data using a nonlinear least-squares procedure with the sine series of the following form at each step:

where f_k denotes independent frequencies detected in discrete Fourier transform of the data and their possible linear combinations. The pre-whitening sequence for RR-Bl is shown in Figure 10. We find a full light curve solution for RR-Bl and produce the fitting results of the light curve in Figure 11.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

According to f₀ = 1.84226, we derive the fundamental period P₀ = 0.54281 days and A₀ = 0.242 mag. Apart from f₀ and its linear combinations nf₀, we detect significant signals that can be expressed as nf₀ + λ_c f_c , where f_c ≈ 1, which can be explained as daily cadence. According to the pre-whitening sequence for RR-Bl, only one independent dominant frequency is detected. So RR-Bl is an RRab star with a strong Blazhko effect.

The additional close side peaks at the fundamental frequency denote long-term modulation (Smolec et al. 2015). In Figure 12, the amplifying pattern is shown. The weak signatures beside nf₀ indicate a possible Blazhko-type modulation, with f_B ≈ 0.009 and P_B ≈ 111.111 days. The signals other than f₀ + f_B are not significant enough compared to others. Larger data sets are needed for a more precise analysis.

Zoom In Zoom Out Reset image size

The observational detection rate of blueshifted emission features in hydrogen lines of RR Lyrae stars is defined as

where N_obs represents the number of RRab stars with featured emission components in LAMOST survey spectra. N_sat denotes the input number of RRab variables in the survey plan that are observed and with the S/Ns of spectra greater than 15. Such a definition of detection rate has an implicit assumption, that all RRab stars possess shock-triggered emission and may or may not be detected in random sampling spectroscopic observations, depending on the observational timing and the strength of the emission feature.

The results of η are as follows: 1.84%—RRab in SDSS; 1.06%—RRc in SDSS; 2.87%—RRab in LAMOST; 1.45%—RRc in LAMOST; 2.97%—RRd; and 2.56% for Blazhko-type RR Lyrae stars in LAMOST. The detection rate of the “first apparition” of RRc stars is significantly lower than that of RRab stars. We suggest that the detection rate in LAMOST is higher than that in SDSS because the sources in SDSS are fainter and the pixel resolution of LAMOST is higher than that of SDSS. For SDSS, we choose the “spCFrame spectra,” and the pixel resolutions are about 4150 at Hα and 4180 at Hβ. For LAMOST, the pixel resolutions are about 7800 at Hα and 8400 at Hβ. The pixel resolution influences our results because one selecting criterion is that both emission and absorption components should contain at least two observational points. We set up two simulations for RRab variables in LAMOST: one on mock spectra to test the performance of our pipeline, the other on the observational strategy to investigate the occurrence of the “first apparition.”

First, we test the performance of our pipeline on mock spectra generated based on the distributions of real observational parameters. We propose that the recognition rate of the shock wave signatures can be described as

where A = “the first apparition is observed on the spectra” and B = “the signal meets the criterion.” We recognize P(A) as the real detection rate, while P(AB) is the detection rate from the search. P(AB) is lower than P(A) owing to finite S/N and our strict criterion of selecting the stars with the “first apparition,” which is described in Section 3.

We generate the mock spectra that have the “first apparition” and apply our 1D pattern recognition pipeline on the mock sample to estimate P(B∣A). The parameterizations of the simulated “first apparitions” are shown in Figure 13. The mock spectrum is regenerated as two Gaussian-like profiles between 6540 and 6590 Å. The central line of the absorption component is fixed at the rest wavelength of Hα. The blueshift of the emission, the ratio between the flux of the emission and absorption component, and the ratio between the FWHM of the emission and absorption component determine the shape of the simulated spectrum. The parameters of the Gaussian profiles are set as random values under Gaussian distribution based on the fitted distributions of the parameters of the selected RRab stars in LAMOST. Random noises are added up to the profiles to reproduce real situations, under Gaussian distribution, where σ is the error of normalized flux. The resolutions and S/Ns of mock spectra are taken from real observations.

Zoom In Zoom Out Reset image size

We apply our pipeline on the mock sample and set up this simulation for 50 times. We only adopt the result of RRab stars in LAMOST, because the pixel resolution of SDSS is low and the search in SDSS relies more on eye-checking. Moreover, the size of the sample of RRc, RRd, and Blazhko-type RR Lyrae stars is not large enough to provide reliable parameter distribution. Figure 14 displays the result that P(B∣A) = 73.48% ± 0.43%. This means that, for a spectrum exhibiting blueshifted hydrogen emission with random flux, FWHM, and blueshift, there should be a 73.48% ± 0.43% chance that our algorithm would recognize.

Zoom In Zoom Out Reset image size

Second, we set up a following simulation on the real observational strategy to investigate the occurrence of “first apparition.” The simulation is based on the real parameters, including detection mode, exposure time, observation interval, observing frequency, and the period of the corresponding RR Lyrae star for each spectrum satisfying S/N > 15.

In the simulation we preset that every star in the testing sample would have a “first apparition.” The exposure time for a spectrum with S/N > 15 is then used as an observing window. Exposure times and observation intervals are obtained from real spectra. In the simulation, we set that the “first apparition” shows up for a short period of time before maximum luminosity in each pulsation cycle. Chadid (2011) reported that the “first apparition” appears over ∼5% of the whole pulsation cycle. Gillet et al. (2019) reported that the “first apparition” presents during ϕ = 0.892–0.929 in RR Lyr, accounting for about 3.7% of the whole period. Their researches are based on bright RR Lyrae stars. Our sample is much larger and contains a lot of faint stars. In the simulation, we assume that the “first apparition” accounts for about 4% of the whole period.

As a matter of fact, too long-time exposures can smooth out the emission features in the hydrogen lines. Therefore, we assume that the survival time of the emission feature and the exposure time should overlap with at least 70% of the exposure time for a clear signal. Each case in which the emission overlaps with more than 70% of the exposure time yields one valid catch for a star. The number of valid catches of one star changes with observation time. The possibility of detecting a star with the “first apparition” is P(B∣A), [1 − (1 − P(B∣A))²], or [1 − (1 − P(B∣A))³] during the time when the star has one, two, or three valid catches, respectively.

For the result, the number of detected RRab stars with blueshifted hydrogen emission in LAMOST from the simulation is about 371.79, while our observational result is 70. The result of the simulation is larger than our observational result. In our simulation we assume that every star in our sample exhibits the “first apparition.” However, not every RRab star has a visible “first apparition” in the real survey owing to variant shock strength. The occurrence of “first apparition” can be estimated as

where O = “the star shows the first apparition.” We define P(O) as the occurrence of “first apparition.” N_theo is the number of detected RRab stars with blueshifted hydrogen emission in LAMOST from simulation. N_obs denotes the number of detected RRab stars with blueshifted hydrogen emission in LAMOST in our survey. Here N_theo = 371.79, N_obs = 70, and P(O) ≈ 18.83%.

We suggest that P(O) is underestimated. First, our simulation does not contain the selecting of Hβ emissions owing to the low resolution of the spectra, so P(B∣A) is overestimated. Another reason for the overestimation of P(B∣A) is that we used the observed distribution but not the intrinsic distribution for emission flux. Moreover, in the progress of fitting we abandon some spectra that may show shock signatures but are hard to identify or provide no valid measurements, or contain no valid Hβ emissions. That is to say, more than 18.83% of the RRab stars show relatively strong “first apparition.”

In this work, we develop a handcrafted 1D pattern recognition searching algorithm and apply it to a large data set of single-epoch spectra of RR Lyrae stars, in order to fetch out the “first apparitions.” Through this survey, we found 33 RRab stars and 10 RRc stars in SDSS and 70 RRab stars, 10 RRc stars, 3 RRd stars, and 1 Blazhko-type RR Lyrae star in LAMOST. Based on the searching results, we set up the first population study of the RR Lyrae variables showing hypersonic shock waves.

We build up the largest database of blueshifted hydrogen emission in RR Lyrae stars. The features of the “first apparitions” are fitted by two Sérsic profiles. We provide measurements of redshift z_e1,α and radial velocity in the stellar rest frame V_e1,α , normalized flux Flux_e1,α , and FWHM of the emission and absorption FWHM_e1,α in the stellar rest frame of the blueshifted hydrogen emission. The distributions of measurements for different types are compared. We provide detailed analysis for the light curve of the Blazhko-type RR Lyrae star with ZTF DR5. We characterize this Blazhko-type RR Lyrae star as an RRab star with strong Blazhko modulations with a possible Blazhko period P_B ≈ 111.111 days.

Finally, we set up two simulations for RRab variables in LAMOST. For the first one, we apply our algorithm on mock spectra to test the performance of our pipeline and to check the influence of the selecting criteria. The other simulation is based on the real observational strategy to investigate the occurrence of the blueshifted hydrogen emission in RRab variables in LAMOST. The result suggests that more than 18.83% of the RRab stars exhibit relatively strong “first apparition.” The nature of RR Lyrae variables will be more and more clear with the enormous volume of upcoming observational data.

The suggestions and comments by the anonymous referee are gratefully acknowledged. We thank the help from Dr. Hao-Tong Zhang, Dr. Zhong-Rui Bai, and Dr. Jian-Jun Chen for getting the single-epoch spectra from LAMOST. We acknowledge discussion with Dr. Anupam Bhardwaj. X.-W. D. acknowledges research support from the Cultivation Project for LAMOST Scientific Payoff and Research Achievement of CAMS-CAS and the Peking University President Scholarship. L.-C.D. acknowledges research support from the National Science Foundation of China through grants 11633005. X.-D.C. also acknowledges support from the National Natural Science Foundation of China through grant 11903045. H.-W.Z. acknowledges research support from the National Natural Science Foundation of China (NSFC) under No. 11973001 and National Key R&D Program of China No. 2019YFA0405504. We thank Mark Taylor for the TOPCAT software. The CSS survey is funded by the National Aeronautics and Space Administration under grant No. NNG05GF22G issued through the Science Mission Directorate Near-Earth Objects Observations Program. The CRTS survey is supported by the U.S. National Science Foundation under grants AST-0909182 and AST-1313422. This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. The Gaia archive website is https://archives.esac.esa.int/gaia. This publication makes use of data products from the Wide-field Infrared Survey Explorer, which is funded by the National Aeronautics and Space Administration. The All-Sky Automated Survey for Supernovae (ASAS-SN) group is partially funded by Gordon and Betty Moore Foundation 5 yr grant GBMF5490. They are also supported by NSF Grants AST-151592 and AST-1908570. ATLAS construction and operations are funded by grants 80NSSC18K0284 and 80NSSC18K1575 under the NASA Planetary Defense Office and Near-Earth Objects Observations program (NEOO). This work also made use of data from the Zwicky Transient Facility project (ZTF), based on observations obtained with the Samuel Oschin 48-inch Telescope at the Palomar Observatory as part of the Zwicky Transient Facility project. ZTF is supported by the National Science Foundation under grant No. AST-1440341. Funding for SDSS-III has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, and the U.S. Department of Energy Office of Science. The SDSS-III website is http://www.sdss3.org/. Participating Institutions of the SDSS-III Collaboration include the University of Arizona, the Brazilian Participation Group, Brookhaven National Laboratory, Carnegie Mellon University, University of Florida, the French Participation Group, the German Participation Group, Harvard University, the Instituto de Astrofisica de Canarias, the Michigan State/Notre Dame/JINA Participation Group, Johns Hopkins University, Lawrence Berkeley National Laboratory, Max Planck Institute for Astrophysics, Max Planck Institute for Extraterrestrial Physics, New Mexico State University, New York University, Ohio State University, Pennsylvania State University, University of Portsmouth, Princeton University, the Spanish Participation Group, University of Tokyo, University of Utah, Vanderbilt University, University of Virginia, University of Washington, and Yale University. The Guoshoujing Telescope (the Large Sky Area Multi-Object Fiber Spectroscopic Telescope LAMOST) is a National Major Scientific Project built by the Chinese Academy of Sciences. Funding for the project has been provided by the National Development and Reform Commission. LAMOST is operated and managed by the National Astronomical Observatories, Chinese Academy of Sciences.

Software: NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), AstroPy (Astropy Collaboration et al. 2018), Matplotlib (Hunter 2007), Scikit-learn (Pedregosa et al. 2011), PyAstronomy (Czesla et al. 2019), TOPCAT (Taylor 2005).