Every Data Point Counts: Stellar Flares as a Case Study of Atmosphere-aided Studies of Transients in the LSST Era

Due to their short timescale, stellar flares are a challenging target for the most modern synoptic sky surveys. The upcoming Vera C. Rubin Legacy Survey of Space and Time (LSST), a project designed to collect more data than any precursor survey, is unlikely to detect flares with more than one data point in its main survey. We developed a methodology to enable LSST studies of stellar flares, with a focus on flare temperature and temperature evolution, which remain poorly constrained compared to flare morphology. By leveraging the sensitivity expected from the Rubin system, differential chromatic refraction (DCR) can be used to constrain flare temperature from a single-epoch detection, which will enable statistical studies of flare temperatures and constrain models of the physical processes behind flare emission using the unprecedentedly high volume of data produced by Rubin over the 10 yr LSST. We model the refraction effect as a function of the atmospheric column density, photometric filter, and temperature of the flare, and show that flare temperatures at or above ∼4000 K can be constrained by a single g-band observation at air mass X ≳ 1.2, given the minimum specified requirement on the single-visit relative astrometric accuracy of LSST, and that a surprisingly large number of LSST observations are in fact likely be conducted at X ≳ 1.2, in spite of image quality requirements pushing the survey to preferentially low X. Having failed to measure flare DCR in LSST precursor surveys, we make recommendations on survey design and data products that enable these studies in LSST and other future surveys.


INTRODUCTION
Differential Chromatic Refraction (DCR) is the apparent shift of a celestial object on the sky due to atmospheric refraction.The magnitude of the shift is determined by the airmass, observing bandpass, and Spectral Energy Distribution (SED, which is dominated by the temperature and spectral composition) of the source.The direction of the shift is toward the zenith.DCR must be characterized effectively to ensure it does not compromise image quality, astrometric inference (Van De Kamp 1967), spectrophotometry (Filippenko 1982), and in surveys where time-domain phenomena are discovered and studied in template subtracted images (Lupton 2007).DCR, however, can also aid scientific inference in a few cases.For example, Kaczmarczik et al. (2009) used DCR to improve redshift mea-surements of quasars in the Sloan Digital Sky Survey (SDSS).While innovative techniques have been previously employed to study the chromatic evolution of stellar flares (e.g.Hedges et al. (2021)), the use of groundbased observations to infer information about flare SEDs via atmospheric refraction remains untested.We consider flares to be suitable candidates for inferential applications of DCR due to their stochastic and short-lived nature, which makes detections in multiple filters difficult, their nature as chromatic transients, and their heightened occurrence on cool stars like M dwarfs, the most common type of star in our galaxy.Currently, there are only a small number of measurements of flare temperatures, which our usage of DCR will enable.
Flares are stochastic, short-lived stellar transients caused by magnetic reconnection at the stellar photo-sphere (Pettersen 1989), most commonly on low-mass stars such as M-dwarfs (dM hereafter).Flares are tracers of stellar magnetic activity, and flare rates have been studied as a correlate with stellar age and rotation (Davenport 2016).Flares can also have a significant impact on planetary atmospheres, as high amounts of UV flux from flares can deplete their ozone layers (Tilley et al. 2019) and are thought to have a role in the formation of life, where they have been identified as both a potential trigger and inhibitor (Ramsay et al. 2021).
Photometrically, they appear as a highly chromatic sudden rise (typically ∼ 3 magnitudes in u band and between 0.5 and 1 magnitude in g band, see Davenport et al. 2012) followed by an exponential decay phase on a timescale ranging from a few minutes to 100-200 minutes in some cases (Yan et al. 2021).While the photometric evolution of flares is well-understood thanks to high-accuracy, rapid-cadence measurements by exoplanet-finding satellites like Kepler (Borucki et al. 2003) and TESS (Ricker et al. 2010), accurate color measurements remain scarce.Statistical inference on flares' origin and evolution mechanisms and their impact is limited by this lack of a large ensemble of information on flare temperatures and temperature evolution, prompting a turn to large next-generation astrophysical surveys, most notably the Vera C. Rubin Observatory Legacy Survey of Space (LSST).
LSST (Ivezić et al. 2019) is the premier ground-based photometric survey of the 2020s.It will continuously scan the whole of the southern sky over a 10-year period, producing over 3 million images with a sub-arcsecond resolution of the whole southern sky in six photometric filters: ugrizy.Based on flare rate measurements in SDSS Stripe 82 (Kowalski et al. 2009), the number of flares will range between 0.4 and 1.4 per image, depending on galactic latitude.While it is certain that Rubin will observe many flares, whether or not we can access scientifically valuable information from those detections is currently an open question.It is expected that the LSST primary survey, referred to as "Wide, Fast, Deep" (WFD), will take 2-3 exposures per pointing per night with a median internight gap of 34 minutes, repeating observations of a field in the same filter every few days (∼ 6 days in r, ∼ 21 days in g, median values, see Bianco & The Rubin Observatory Survey Cadence Optimization Committee 2023 for details on the most recent LSST observing plans).Given that Yan et al. 2021 reported median rise and decay times of flares on Sun-like stars in the Kepler catalog to be 5.9 minutes to 22.6 minutes respectively, LSST will observe most flares in the WFD with 1-2 data points per detected event, rendering the time-resolved photometric analysis enjoyed by e.g.Kepler impossible.The Deep Drilling Fields (DDF) will offer an opportunity for time-resolved flares and to identify even faint dM flare precursors.Roughly 5% of the total survey time will be devoted to continuous observations of five selected pointings.While the observational details for the DDFs are yet to be defined, in one hour over 100 consecutive exposures can be taken, over 20 in each of five filters of one 9.6 sq degree field of view.1When co-added with the main survey observations these images will extend the 5σ depth down to r ∼ 27.5, allowing time-resolved observations of even the faintest dM s in multiple filters.It must be noted that all DDFs are extragalactic fields, where the flare rate is expected to be suppressed.2Yet flares are common phenomena even in the extragalactic sky, and we expect hundreds to be detected in each DDF by the end of the 10-year survey.The LSST DDF program will enable traditional flare studies.However, by exploiting the high astrometric accuracy of Rubin, here we show that the Rubin WFD, with its unprecedented sample size, offers a valuable opportunity for flare science.
It is already expected that DCR will aid extragalactic studies of quasars (Yu et al. 2020;Richards et al. 2018) with Rubin, and we here propose that employing DCR as a scientific tool can be extended to variable and transient phenomena with Rubin.We use stellar flares as our case study and demonstrate in this paper that by taking advantage of Rubin's pristine image quality and astrometric precision to measure the difference in source location on the sky between quiescence and event due to DCR, we will be able to indirectly obtain information regarding the color, and therefore the temperature of the flare, from even a single point detection.We will show that an astrometric shift should be apparent in LSST astrometry for a star when flaring compared to when quiescent, with the effect being especially pronounced for a hot flare on a cool star observed at high airmass; but we further demonstrate that the effect will be measurable by Rubin even at its typically low airmass range.This means that, despite the sparse photometric sampling of Rubin, the DCR toolset can be used to leverage its high volume of flare detections in order to infer flare temperatures across an expansive sample.
In section 2 we describe how DCR is calculated and how we model the expected excess DCR during a flare with some temperature, at some airmass, in a given fil- ter.In section 3 we compare the expected DCR for flaring stars with the capabilities of the Rubin system and data analysis pipeline.In section 4 we describe our failed search for excess DCR produced by flares in surveys considered precursors to LSST, ZTF (Bellm et al. 2018) and Deep Drilling in the Time Domain with DE-Cam (Graham et al. 2023) and identify the bottlenecks that impaired a successful DCR measurement.In section 5 we discuss the technical requirements for measuring the ∆DCR from flares: the excess zenith-bound displacement of a flaring star compared to the apparent shift expected based on its quiescent SED, and provide recommendations to maximize the potential of this method for LSST and other surveys.In section 6 we summarize our conclusions and explain how this method can be embedded within Rubin's software ecosystem to support studies of flares and other chromatic phenomena throughout the survey lifespan.

METHODOLOGY
Following (Kaczmarczik et al. 2009), we can describe the formalism of DCR in the following four steps: defines λ eff as the effective wavelength of bandpass j, where S j is the transmission function of j, and f λ is the spectral flux of the source (or SED) in units of W m −2 µ −1 ; for an observed source then, the filterdependent refraction index n λ can be calculated from (2) R, the total angular deflection of the source due to DCR (in arcseconds), is then calculated as: where Z is the zenith angle, and where the airmass X is commonly approximated as X = sec(Z), which assumes a homogenous, plane-parallel atmosphere.This approximation is valid only for 60 • ≤ Z ≤ 75 • , so our calculations involving airmass will be restricted to X ≤ 2.1 (although we note that the LSST survey strategy includes observations at higher airmasses, see section 3 and subsection 5.2). Figure 1 illustrates the amplitude of DCR produced by a 10,000K blackbody source for a range of airmasses in the six LSST bandpasses.In section 4 we will demonstrate that even flares with moderately high temperatures observed at airmasses 1.05 ≤ X ≤ 1.2 will produce an astrometric shift that will be detectable at the LSST precision level.
2.1.The DCR effect induced by flares: ∆DCR Using the DCR formalism described in section 2, we created an expository model to illustrate the DCR produced by a flare.We assume the star is an M5 dwarf and model the quiescent SED of the source with a template spectrum spanning the optical and near-infrared wavelength ranges built by Davenport et al. (2012), with the optical component using observations from the Sloan Digital Sky Survey (SDSS).We clip this template spectrum at 12, 000 Å as this was sufficient for full coverage of the LSST ugrizy bands.The source in its quiescent state will already suffer from a DCR effect, but we expect that the temperature contrast between flare and quiescent states will determine a significant change in the apparent star position.The difference between the displacement in quiescent and event state is the excess displacement produced by the flare which we refer to as ∆DCR.A schematic illustration of ∆DCR is shown in Figure 2.
To simulate a flare, we add a blackbody spectrum at temperature T to the spectrum of the star (Figure 3).In the absence of a well-characterized spectral energy distribution (SED), flare spectra are canonically approximated as a 9,000-10,000K blackbody (Osten & Wolk 2015).Thus, in this simulation the spectra are normalized such that the blackbody component at the canonical flare temperature (T = 10, 000K) has the same total energy as the dM spectrum when both are integrated over the SDSS optical range (3, 850 − 9, 200 Å).We also explored the effect of changing the fraction of the stellar surface covered by the flaring region, referred to as a "filling factor".Assuming our energy-based calibration to be representative of a filling factor f f ∼ 0.05 and applying a flat scaling to the blackbody spectrum for other f f values, we tested filling factors ranging from 0.05 to 0.2 and found that we are not sensitive to this parameters for f f ≥ 0.05.The composite spectrum is convolved with the chosen LSST filter to calculate λ eff in Equation 1, allowing us to estimate the angular displacement from DCR during both quiescence and event, and their difference, or ∆DCR.

Measuring ∆DCR with the Parallactic Angle
Many astrophysical and observational effects may cause the apparent position of a star to change between two images of the same area of the sky.Although stars' proper motion is typically too slow to have a measurable effect in images collected within days, weeks, or months, images are however subject to various distortion effects such that a star may appear to be offset between two images, especially if observed in a different position of the CCD plane, with different telescope rotation, etc.These observational systematics should be corrected, and we will return to these spurious effects in section 4.These components of the star's displacement should have no preferential direction while the DCR-induced displacement should be strictly in the direction of the zenith.To determine the DCR-induced change in position of a given source between detections, we will measure the components of the apparent motion toward, tangential to, or away from the zenith direction.The parallactic angle is defined as the angle between two great circles, one passing through the source and the zenith, and the other passing through the source and the North Celestial Pole.It is calculated according to Meeus (1998) as: where δ is the object's declination, h is the object's hour angle, and ϕ is the geographic latitude of the observer's location on Earth.The component of the source's motion in the direction of the zenith (which we call d ∥ ) is then calculated as: ) where α is the object's right ascension, and P 2 is the parallactic angle of the source in its second position.

DETECTION POTENTIAL FOR LSST
Using the methodology described in section 2, we simulate the ∆DCR produced by a flare on an M5 star for a variety of airmasses and flare temperatures across all six Rubin bands.The ∆DCR for a fixed flare temperature at 10,000K for all 6 filters in the 1.1-2.1 airmass range is shown in Figure 4.The Rubin system specification requires a single-image absolute astrometric accuracy of Rubin of 0.1 arcsec (Ivezić & The LSST Science Collaboration 2018).This figure demonstrates that a 10,000K flare's ∆DCR is detectable in g-bands even at moderate airmasses; ∼ 8% of all WFD images will be in g band.
While the effect is also detectable in u-band, higher airmasses are required to generate a comparably large shift.This may be counterintuitive, as the DCR effect is more pronounced in bluer wavelengths.However, the ∆DCR is proportional to the change in λ eff during the flare compared to quiescent state.The quiescent SED has a prominent redward slope in the g band wavelengths which, integrated through the filter leads to a λ eff,quiescent (g) = 4, 988.4 Å which becomes λ eff,flare (g) = 4, 740.9 Å when the black body contribution dominates the SED during flare.Conversely, the quiescent SED in u band is more flat and only slopes redward significantly where the filter transmission is already very low.Thus the shift in λ eff is smaller (λ eff,quiescent (u) = 3, 728.0 Å compared to λ eff,flare (u) = 3, 655.8Å).
The g-band evolution of the ∆DCR with flare temperature is shown in Figure 5, showing that the ∆DCR grows rapidly between 5,000-10,000K, then flattens off at higher temperatures.The minimum system requirement and stretch goals for both absolute and relatve astrometric accuracy of Rubin is indicated on Figure 5, demonstrating that DCR could be used to probe flare temperatures as a low as 4,000K, depending on the type of data products used for the analysis (see section 5).
To maximize image quality, the automated LSST scheduling system will preferentially observe at low airmass by design (Ivezić & The LSST Science Collaboration 2018).The dependence of this effect on airmass then begs the question of whether or not LSST will perform observations at sufficiently high airmass for the effect to be detectable.However, the most recent fiducial simulation of the LSST (baseline v3.0 10yrs, Bianco & The Rubin Observatory Survey Cadence Optimization Committee 2023) suggests that nearly 10 6 WFD visits will be performed above airmass 1.25 (Figure 6), which, as we have shown in section 3, is sufficient to produce a detectable ∆DCR at typical flare temperatures (we will return to the LSST design in subsection 5.2).
These simple calculations indicate that Rubin should be capable of detecting even relatively cool flares at typical airmass on a typical M5 dwarf.However, contamination from image warping, chromatic aberation, DIA misalignments, etc. will complicate this simplified expectation.We tested our method on detected flares in precursor surveys selecting as flare candidates transient events with short duration (≲ 2-hour) and a change of magnitude of at least 0.4 mag in g band.

Zwicky Transient Facility
The Zwicky Transient Facility (ZTF, Bellm et al. 2018) was designed to detect transient objects across the entire northern hemisphere.It is considered a precursor survey for LSST and it will deliver images and lightcurves for 3,750 sqdeg/hour with alerts delivered in real-time and a significant fraction of its data made available without proprietary restrictions.The exposure time is similar to that of Rubin's LSST (30 seconds + 10 seconds readout) to a single image depth of ZTFr 20.5, and three filters are available (ZTFg, ZTFr, ZTFi) also with similar throughput to the Rubin filters.With a smaller footprint by a factor of four and three filters the revisit time is nominally nearly one order of magnitude shorter than LSST's, enabling multiple observations of the same flare to be collected.
Notable differences, however, for the purpose of our science, are the overall data throughput, about 10 times smaller than Rubin as measured in bytes of data, but which corresponds to a factor of nearly 100 fewer targets due to the decreased limiting magnitude, and a significant decrease in image quality, measured as Point Spread Function (PSF), leading to decreased astrometric accuracy.The instrumental pixel scale is 1.01 arcsec/pixel, compared to Rubin's ∼0.2 arcsec/pixel, and the median seeing-limited PSF is ∼2 arcsec, compared to ≲ 1 arcsec expected for Rubin.This forced us to limit the study to flares observed at high airmass to enable DCR-based temperature estimates, further decreasing the size of our flare sample.In our preliminary analysis, we inspected a sample of 17,000 bright dM in ZTF; 2.3% showed possible large flares in data collected airmass above 1.4,for a total of 414 candidates.To increase the confidence that the observed brightening was indeed a flare we required more than one observation within the event (i.e.within ∼ 2 hours).Of these 414 flares, only one was captured with more than one data point.ZTF's reconstructed astrometric accuracy per science image with respect to Gaia DR1 is ∼ 0.045 − 0.085 ′′ for sources extracted at a 10-σ limit, however, the seeing Full Width Half Maximum (FHWM) in this pair of images was 3.474" for the first flare epoch and 2.910" for the second.The seeing turned out to be the bottleneck in the application of our method; the astrometric solution generated by the ZTF pipeline measured a small ∼ 0.2 ′′ displacement between the first and second image compared to other stars in the field, too small compared to the seeing to confidently determine whether or not the star moved towards the zenith.

DECam Deep Drilling Fields
The Dark Energy Camera (DECam; Flaugher et al. 2015) at the Blanco 4-meter telescope is the instrument that enabled the Dark Energy Survey (Dark Energy Survey Collaboration et al. 2016).The camera itself is a single chip of the same kind as those that will constitute the LSST camera mosaic, making DECam data naturally comparable with LSST's, with high image throughput, similar image quality (0.263 arcsecond/pixel resolution, with a telescope located near the site of LSST leading to similar sky properties), and similar system wavelength coverage (grizY filters).Graham et al. (2023) used DECam to survey two of the LSST DDFs: COSMOS (Scoville et al. 2007) and ELAIS (Oliver et al. 2000).This led to a precursor survey of the LSST DDFs with 5-sigma limiting magnitudes r ∼ 23.5 mag (single exposure).However, in addition to the shallower depth, there are other significant differences.This DECam DDF program was one of the co-founding members of the DECam Alliance for Transients (DECAT), within which multiple PIs of DECam programs pooled their time to enable dynamic queuelike scheduling and time-domain science.With a field of view of 9.6 deg 2 , the LSST DOE camera will cover each field with a single pointing.Conversely, the DECam field of view (3 deg 2 ) is smaller and DECam DDF program used three adjacent pointings in COSMOS and two in ELAIS.Every night of observations cycled through each pointing five times, obtaining a sequence of gri images per pointing.However, because the fields are covered with multiple pointings, 10 or 15 observations per night are obtained in the area where the pointings overlap (see Figure 7).In contrast, LSST will obtain tens of consecutive images in five filters on each night when a DDF will be observed.Taking these differences into consideration, the DECam DDF program is an effective test-bed for our DCR flare studies, but not for temporally resolved flare investigations.Importantly, the data processing pipelines for DECam DDF and LSST also differ: LSST will process its data with dedicated pipelines (Bosch et al. 2018), while the DECam DDF fields data are processed with existing software to detect transient events (Graham et al. 2023).Briefly, the difference-image analysis (DIA) pipeline implemented in this survey ingests raw images directly from the NOIRLab data archive and performs standard data reduction procedures.Source Extractor (Bertin & Arnouts 1996) is used to detect all sources in the image, SCAMP (Bertin 2006) is used to calculate the astrometry for each chip and match each source with stars drawn from the Gaia DR2 catalog, and then SWARP (Bertin 2010) is used to solve for the world coordinate system (WCS) using these objects.Object catalogs are generated from the reference images with Source Extractor and aligned with SCAMP and SWARP.The image subtraction is done with HOTPANTS (Becker 2015).Notably, there is no DCR correction in the difference image analysis as implemented for the DECam DDF fields.Source Extractor is used on the resultant difference image to identify residual signals and extract fluxes via forced photometry.Lastly, using the algorithm described in Goldstein et al. (2015), each detec-tion is assigned a "real/bogus" score3 , ranging from 0 (bogus) to 1 (real).Objects (i.e.detections) within 2" of a previously detected object are associated with the same candidate ID.
We elected to focus our analysis on the COSMOS g band sample (recall our preference for g band from section 3).We selected a subset of the publicly available DECam DDF transient catalogs4 in the COSMOS field, requiring the following: • at least 1 g-band detection, • all detections within 0.5 days (single-night), • a mean "real/bogus" score RB ≥ 0.6 across all detections.
These cuts resulted in a sample of 1230 candidates, and the coordinates of this g band detections' sample are shown in Figure 7.To identify flares in the data and measure the change in magnitude during the event we applied the following additional cut, resulting in a final sample of 1015 COSMOS flare candidates: • at least two g-band detections, each with real/bogus score RB ≥ 0.55 We calculated d ∥ , the component of the candidates' motion along the parallactic angle (toward or away from the zenith), for all pairs of g-band detections with the same on-sky association.d ∥ is plotted against the measured change in g magnitude in Figure 8 in order to visualize the sample in an informative phase space and select interesting candidates for additional inspection.We expected flares to be found in the upper left quadrant of the plot, i.e. pairs of detections where the brightness of the source increased (relative to its initial difference image detection) and also moved on the sky in the direction of the zenith at the time of observation, or lower right quadrant for dimming flares moving away from zenith.In addition, we favor candidates observed at higher airmasses, but with little to no change in airmass between observations, as a change in airmass dominates the chromatic contribution to DCR (see Figure 1).However, as shown in the left panel of Figure 8, few observations of the DECam DDF fields occur at airmass above X = 1.2, and none of these high-airmass detections displayed sufficient change in brightness to be considered for additional inspection.We selected three objects for further investigation (marked with circles in Figure 8): the two most extreme outliers in the upper-left quadrant that are ostensible rising flare candidates, and an object with large zenithbound motion but no magnitude change, as a control.All of them show a small change in airmass between the two observations in a pair, as desired.While a few DE-Cam DDF transients that passed our cuts are detected at high airmass X > 2 (Figure 8, left panel), none of them have significant enough magnitude changes to be considered as flare candidates.All of our three objects of interest are observed at airmass X ∼ 1.2.Thus, in order to show a significant ∆DCR, if they were flares, these first two candidate events would have to be extremely hot.The coordinates, detection times, g band photometry, and d ∥ for these candidates are shown in Table 1.To ascertain the stellar origin of these transients we look for a quiescent-state source in the tem-plates. 5The DIA triplets of the three candidates circled in Figure 8 are shown in Figure 9. From left to right, each triplet contains the "new" image (also called the "direct" or "science" image, obtained by the DECam DDF program), the "reference" (a template built from archival DECam images obtained in previous years), and the "difference" image (PSF-matched subtraction of the reference image from the new image).For none of the three candidates were we able to find a precursor source in the reference image.It is possible that the quiescent sources were not present in the reference image because they were beyond the survey magnitude limit.However, the templates reach about one magnitude deeper than the search images (Graham et al. 2023), and a flare with 1% surface coverage on an M3 dwarf, would generate a brightening of ∼0.62 magnitudes compared to the qui- escent source in g band (Davenport et al. 2012), which should enable the detection of the quiescent source even for transients detected close to the single image 5 − σ limit (see Table 1).Since dM are typically bright in the infrared (IR) wavelengths, rather than the optical wavelengths, we also searched for a quiescent source in the IR by examining the WISE data (Wright et al. 2010) at the candidate's coordinates (Figure 10), but likewise did not see a source in these images.This is not surprising: given the effective magnitude limit of WISE W1 (16.6 mag; Cutri et al. 2012), a faint transient in the DECam DDF data, such as our three candidates, is expected to be at the detection limit of the WISE data in quiescent state.This raises the possibility that the candidates may not be stellar sources, but rather Solar  1 is plotted (where only one circle is visible, as for DC21engi, the detections are overlapping).No source can unambiguously be identified at the location of the transients, except for DC21engi.
System Objects or extragalactic transients.We checked additional DECam alerts to see if the transients were detected at later epochs: extragalactic transients have typically slower evolutions, so they should remain detectable in successive nights.Inspecting data with the DECAT LBL Viewer, we found no subsequent alerts detections for any of the four candidates.Finally, we checked if any known Solar System Object (SSO) is expected at the coordinates of the detection by inspecting the IAU Minor Planet Center (MPC) catalogs using the Minor Planet Checker6 and we found SSOs within one arcminute (the minimum search radius in the MPC) of the candidate coordinates, at the time of detection for all three candidates.This leaves asteroids as the most likely source of contamination given the brightness of the transient and the fact that each only appeared in a single night of observations.
We then took a more traditional approach and began with ascertaining if any of the DECam DDF detections selected by our cuts were indeed flares, then  (Brown et al. 2021).We selected Gaia sources with G bp − G rp > 2 and M G > 5 to choose K − M dwarfs and performed an additional cut on parallax error (≤ 1 ′′ ) over three areas corresponding to the three DECam pointings that compose the DDF's COSMOS field.This led to the selection of 6993 Gaia DR3 sources.
We then crossmatched the positions of the Gaia sources to within 1.0" of the initial 1230 DECam DDF candidate positions using the match to catalog 3D method of the astropy.coordinates.SkyCoord object.This crossmatch left us with a single candidate, DC21engi, whose DIA triplet for two g-band detections is shown in Figure 9 and whose WISE W1 image is shown in Figure 10.We analyzed the astrometric changes of DC21engi between quiescence and event by taking the Gaia coordinates as the bonafide quiescent position.In Figure 11 we examine the relative positions of the first and second DIA g band detections, using the parallactic angle to indicate what the motion of the source would be if it moved directly toward zenith as expected for DCR.While the initial detection (g ∼ 22.6) of the event and the associated motion from its Gaia position is suggestive of DCR (i.e. in the zenith direction), the second detection (g ∼ 23.0) defies the expectation that the source would move back towards its quiescent location as the photosphere cools and the object position is measured in a direction perpendicular to the parallactic angle.The bottleneck here is likely the precision of the DECam astrometric solution, which dominates the uncertainty in the position: astrometric match requirements for images in Graham et al. (2023) are set by median astrometric residual across the entire field which must not exceed 0.15".This margin is shown by the errorbars on the detection positions, indicating that this candidate does not conclusively demonstrate motion (produced by DCR or otherwise) during the transient event.

TECHNICAL RECOMMENDATIONS FOR ∆DCR DETECTION
The difficulties encountered in the process of validating this method on precursor surveys emphasize the extremely high image quality and astrometric fidelity requirements demanded by such a method.We propose a set of astrometric accuracy benchmarks to maximize atmosphere-aided studies, continuing to use flares as our case study.Figure 12 shows the ∆DCR induced as a function of flare temperature and airmass in LSST g-band, assuming a distribution of flare T ef f as measured by Howard et al. (2020) using a combination of Evryscope (Law et al. 2015) and TESS data, modeling flares as blackbodies, and a per-visit airmass as in the baselinev3.010yrs simulation of LSST.

Astrometric requirements
. By combining these distributions, we can estimate the probability of measuring ∆DCR in LSST as a function of flare temperature.Figure 13 shows the conditional probability of producing a ∆DCR-induced astrometric shift detectable in LSST during a flare given the flare temperature distribution in Howard et al. (2020): P (X > X crit |T ef f > T ), where X crit (T, q) is the minimum airmass necessary for a flare of temperature T to produce a ∆DCR greater than the astrometric accuracy q.The probability of measuring ∆DCR from a flare with peak temperature at least 10,000K is 56.8% for a survey with astrometry accurate to within 10 mas, but this probability falls to 27.1% at the 100 mas level and 4.1% at the 200 mas level.
New methodologies for increasing astrometric accuracy have been proposed that could improve upon the 100 mas requirement set in Ivezić & The LSST Science Collaboration (2018).Fortino et al. (2021) leveraged Gaia astrometry on stars visible in both surveys and Gaussian process regression (GPR) to reduce the astrometric variance caused by atmospheric turbulence.Testing on the orbit of trans-Neptunian object Eris (r ≈ 18.5) they found that GPR corrections reduced the root-mean-square (RMS) residuals in riz band from 10 to 5 mas.The observations used to perform the correction to the Eris orbit were done at median airmass X = 1.29 which exceeds the median airmass of baseline v3.0 10yrs setting the expectation that similar improvements could be observed in ∆DCR relevant LSST observations.Furthermore, they state explicitly that the GPR technique should be applicable to LSST, and that Rubin's larger aperture compared to DECam means that more stars will be available for the fit in the turbulence-dominated leading to improved results.We note that the RMS reduction is proportional to the number of Gaia stars available in the image, indicating that better results (up to a factor of 5) can be expected in Galactic fields, where the rate of flares is enhanced (see Fortino et al. 2021 fig. 5) Full characterization of the uncertainty in both the quiescent and flaring position of the source is necessary for DCR-based inference of flare temperature.If prompt follow up is not required, or even not possible for very short duration events such as flares, the Rubin annual Data Releases can be used for the analysis, leveraging the high Rubin's internal relative astrometric accuracy (Figure 5).This analysis can also be performed using the world-public alert packets sent to brokers within 60 seconds of observation and the Prompt data products released 24 hours after observation (Jurić et al. 2022).However, these data products do not include positional uncertainties, thus direct analysis of the images is required.

The impact of observing strategy choices
The observing strategy of Rubin LSST, which is still being finalized at the time of writing of this manuscript (Bianco et al. 2021), has led to the current survey strategy recommendation (Bianco & The Rubin Observatory Survey Cadence Optimization Committee 2023).The the survey strategy is simulated via the LSST Operation Simulator (Delgado et al. 2014;Delgado & Reuter 2016;Naghib et al. 2019); the median and maximum airmass of the LSST current observing strategy proposal (baseline v3.0 10yrs) is shown in Figure 14.Highairmass visits are favorably distributed across the WFD footprint, with visits above airmass 2.0 occurring at least once anywhere in the WFD footprint.
Four overarching science goals set the survey requirements and the high level design of the survey strategy (Ivezić et al. 2019): probing dark energy and dark matter, building an unprecedentedly complete orbital catalog of Solar System objects, exploring the transient and variable sky, and advancing our understanding of the Milky Way via the resolved stellar population.The system is designed to simultaneously maximize field of view and depth, i.e. to maximize the volume of space-time that can be surveyed, leading to a fast survey that can scan the entire southern sky in ∼ 3 days, which maximizes the discovery potential in time domain.Many of the science deliverables of LSST, however, require highresolution imaging and exquisite image quality, which pushes the survey design to prefer low airmass observations.Nonetheless, in the complex optimization exercise that weighs in the needs of many science deliverables, it is inevitable that a subset of the observations will be performed at medium, and even high airmass.Over the course of several years, as the observing strategy for LSST has evolved under the guidance of many community contributions (Bianco et al. 2021), we observe that the power in the high-airmass tail of the distribution of airmasses has increased.As the optimization of the observing strategy is refined under a complex net of science goals, additional constraints on the pointing com-pete with optimizing the airmass choice for image quality.This tail of high airmass observations, however, is welcomed by our science case and all atmospheric-aided studies (Richards et al. 2018;Yu et al. 2020) as it increases the probability of measuring (∆)DCR.Figure 15 shows the evolution of the distribution of airmasses over five recent opsims (simulations of the LSST 10-year survey).The most recent recommendation for the Rubin LSST observing strategy (Bianco & The Rubin Observatory Survey Cadence Optimization Committee 2023) has led to the most favorable baseline among those considered in this work for DCR-assisted temperature and SED measurements.

Data Products Considerations
The astrometric requirements of LSST are sufficiently demanding that DCR corrections may be required to achieve them (Ivezić & The LSST Science Collabora-Figure 13.Probability of detecting ∆DCR induced by a flare at or above peak effective temperature T for four different astrometric accuracy limits.The probability of DCR detection is measured as the conditional probability P (X > Xcrit|T ef f > T ), where Xcrit(T, q) is the minimum airmass necessary for a flare of temperature T to produce a ∆DCR greater than the astrometric accuracy q.The temperature distribution follows the measured flare temperatures in Howard et al. (2020).The small number of flares at large T causes the discontinuities in the probability at T ef f > 20, 000K.
tion 2018, Swinbank et al. 2020), and understanding how and when these corrections are applied is necessary to develop a schema for measuring ∆DCR with the Rubin data products.While the position of a flaring star will be available in an alert prompted by a flare, the centroid associated with the alert will be measured on the difference image.Difference imaging is at the core of Rubin's transient detection, and as Rubin is not equipped with an atmospheric dispersion corrector, DCR presents a challenge for image subtraction methods, most notably in the form of "dipoles" in the subtracted image caused by mis-subtraction of sources observed at different zenith angles.At the Difference Image Analysis (DIA) stage, DCR-matched templates can be used to mitigate these dipoles.This method (described in detail in Sullivan 2018) iteratively forward-models the unrefracted sky and uses the result as a template for image subtraction.However, in practice, LSST will not produce enough data in the bluer bands in the first year of operations to produce DCR-matched templates (Swinbank et al. 2020).
The LSST Science Pipeline does not currently implement wavelength-dependent PSF modeling (e.g.Meyers & Burchat 2015) to account for DCR effects, although this will probably be added early in the survey.DCR produces both a bulk shift of the PSF centroid (first or-der effect) as well as a second order effect on the shape of the PSF, increasing the dispersion along the zenith direction.Kolmogorov turbulence in the atmosphere leads to a linear isotropic contraction or dialation of the kernel.
In Figure 16 we explore the impact on DCR correction approaches on the resulting DIA -and thus our ability to infer flare properties -for a simulated flaring M dwarf.The initial model (top row) shows the most ideal scenario, where the T=10,000K flare occurs at the same airmass as a template.The PSF is modeled as a Gaussian with F W HM = 0.7 arcsec seeing.As expected, the resulting DIA shows a bulk shift toward zenith from the blue flare of ∼ 0.10 arcsec or ∼ 0.5 pixels at the LSST system plate scale at an airmass of X = 1.22.This model represents perfect DCR correction, since the template and flare occur at the same airmass, and would allow us to correctly infer the flare's temperature in an ideal (i.e.noise free) scenario.
The second scenario in Figure 16 demonstrates the flare occurring at a higher airmass (X = 1.74) than the quiescent template.If appropriate DCR correction for the red dM SED is applied, then the resulting DIA is again due to the blue flare flux alone.Note that the DIA shift is larger and more broadened compared to the first scenario, due to the increased airmass, which naturally results in a more precise constraint on the flare temperature as expected (e.g.see Figure 4 and Figure 5).This scenario represents the pipeline applying perfect DCR corrections for all known source SEDs at every airmass, making ∆DCR measurable from DIA.
The third scenario shown in Figure 16 demonstrates an incorrect or ad-hoc correction for DCR.Here we have assumed that the flare observation was at a high airmass, and the field was corrected for some average DCR effect without any knowledge of the quiescent source SED.For this demonstration, we assumed the dM had a DCR correction applied for a solar-type G2V star.The impact here is an over-correction for DCR for the quiescent red dM source.The blue flare flux still appears in the resulting DIA in Figure 16f.However, the resulting DIA signal is distorted from the expectation (Figure 16e) due to the incorrect DCR correction, which would negatively impact the accuracy of the inferred flare temperature.
This result of these DCR simulations show that 1) the ∆DCR signal should be a good indicator in DIA products of flare activity for low-mass stars at moderate airmasses, even with imperfect DCR corrections, and 2) careful and accurate DCR corrections based on the quiescent stellar SED are required to correctly infer flare temperatures using DIA.Plots produced within the Rubin Metric Analysis Framework.(Jones et al. 2014).To produce these maps, a sky segmentation into healpixels with resolution 64 (Górski et al. 2005) 6. CONCLUSIONS We have proposed a method to extract spectral information from chromatic transients by taking advantage of the change in differential chromatic refraction (DCR) between quiescence and event states, or ∆DCR, using flares as a case study.We used a composite spectrum derived from SDSS dM spectra to demonstrate the capability for LSST to measure the ∆DCR induced by flares with effective blackbody temperature 10,000K at or above X=1.2 in g band, enabling ∆DCR temperature measurement of flares over a significant fraction of total survey operations.Depending on the type of data products available for the analysis, flare temperatures as low as 4000K could also probed by Rubin using DCR.In order to prepare to deploy this methodology on the Rubin pipeline, we tested its efficacy on two precursors to LSST, the Zwicky Transient Facility (ZTF, (Bellm et al. 2018)) and the "Deep Drilling in the Time Domain with DECam" program (Graham et al. 2023).Our initial tests of the ZTF data quickly revealed a critical issue with using ZTF as a validation testbed, as the ZTF pixel scale (1.01"/pixel) is an order of magnitude larger than the astrometric change we expect to be induced by the ∆DCR produced by a typical flare and the typical seeing exceeds 2".While the system properties and observing conditions of the DECam DDF observations are more favorable, and indeed comparable to the LSST's, the processing pipeline did not lead to accurate enough astrometric solutions to enable ∆DCR measurements.We estimated how accurate an astrometric solution must be in order to be able to measure ∆DCR, given an observed flare's peak temperature and found that a survey must be accurate to within 0.01" to measure 56.8% of ∆DCR events produced by flares above 10,000K.This initial exploration of studying stellar flare temperatures using DCR in the context of Rubin LSST has demonstrated the potential for next-generation surveys to use incredible advances in image quality, astrometric solutions, and data throughput to break new ground using novel, albeit unconventional methods.However, the challenges encountered in our search for ∆DCR in precursor surveys emphasize the strict technical requirements for such a method to deliver reliable results.Future work on this project will involve testing this method by injecting artificial flares into simulated LSST-like images and testing the method's ability to precisely recover positional offsets and accurately infer flare temperatures using Rubin data products.

Figure 1 .
Figure1.Angular deflection from true position for a source with a 10,000K blackbody spectrum for airmass 1.1 ≤ X ≤ 2.1 in all six LSST bands.A zoomed section is shown to reveal the separation between curves corresponding to separate bandpasses.

Figure 2 .
Figure 2. Light incident from a star is deflected by the atmosphere.The amount of deflection depends on the color (i.e.temperature) of the source and the amount of atmosphere the light passes through.The chromatic change during a flare event should produce an excess in the normal DCR at quiescence, labeled as ∆DCR in the figure.

Figure 3 .
Figure 3.In pink, the spectrum of an M5 dwarf (composite spectrum built partially from SDSS observations, Davenport et al. 2012).In blue, the spectrum of a 10,000K blackbody to simulate the flare.In purple, the sum of the two aforementioned spectra, representing the spectrum of the dM during the flare event.The blackbody and dM spectra contain the same total energy over the SDSS optical range (3, 850 − 9, 200 Å).The transmission functions for the LSST ugrizy photometric system are shown in grey.

Figure 4 .
Figure 4. Expected magnitude of the ∆DCR effect for a flare SED approximated by a 10,000K blackbody as a function of airmass and filter in the Rubin Observatory ugrizy observing system(Olivier et al. 2008).Blue coloring corresponds to a ∆DCR shift detectable by Rubin, and red coloring corresponds to an undetectable shift, given the absolute astrometric accuracy goal of 0.1 arcsec.

Figure 5 .
Figure 5.The magnitude of a star image displacement during a flare compared to the quiescent star position as a function of flare temperature in LSST g-band for four different airmass values.The minimum requirement and stretch goal for absolute astrometric accuracy of Rubin are shown by the the black and grey dashed lines, respectively.The minimum requirement and stretch goal for relative astrometric accuracy of Rubin are shown by the black and grey dash-dot lines, respectively.

Figure 6 .
Figure 6.Airmass distribution of the current LSST survey strategy proposal (baseline v3.0 10yrs, Bianco & The Rubin Observatory Survey Cadence Optimization Committee 2023).Cumulative distribution showing the number of visits at or below a given airmass.Plot produced within the Rubin Metric Analysis Framework (Jones et al. 2014).

Figure 7 .
Figure 7. Location of DECam DDF flare candidate g-band objects in the COSMOS field.

Figure 8 .
Figure 8. Scatter plots of 1015 flare candidates in the DECam DDF survey of the COSMOS field (for selection criteria discussed in subsection 4.2).The x axis is the change in g magnitude between the objects and the y-axis is the change in the component of the source's movement in the direction of the zenith.Each point represents one unique candidate ID with at least two DIA "objects" (two distinct observations).The change in magnitude and position are calculated as the difference between the two objects with the largest absolute change in magnitude.Left panel: For observations with airmass X < 1.2, points are shown as grey circles.Points are colored by the airmass of the initial detection for airmasses X > 1.2.Right panel: Each point is colored by the change in airmass between the two objects and histograms show the marginalized distributions of g magnitude change (top) and zenith-bound displacement (right).Three candidates of interest, as described subsection 4.2, are circled.The single candidate whose quiescent counterpart can be identified in the Gaia DR3 dataset, DC21engi, is marked by an arrow.

Figure 9 .
Figure 9. Image triplets for DECam DDF candidates, as shown by the DECAT LBL Pipeline Candidate Viewer.From top to bottom, the candidates are: DC21fjnb, DC21fygb, DC21flly, and DC21engi."New" denotes the search image, "Ref" the reference image, and "Sub" the difference image.A quiescent source consistent with a dM is seen in the reference image for DC21engi.

Figure 10 .
Figure10.WISE W1 band images; the coordinates of candidate flares DC21flly, DC21fjnb, DC21fygb, and DC21engi are marked as labeled.In all panels, North is up and East is left.Each panel is 3.8' on the side.A 12" radius circle centered on each detection reported in Table1is plotted (where only one circle is visible, as for DC21engi, the detections are overlapping).No source can unambiguously be identified at the location of the transients, except for DC21engi.

Figure 11 .
Figure 11.Schematic of the position change of candidate DC21engi during a DECam transient event captured in two detections, at times t1 (2021-04-09 02:23:03.576)and t2 (2021-04-09 02:41:18.079).The Gaia DR3 coordinates of the source are shown in black, and the coordinate axes are shifted such that the Gaia position is at the origin.The source is assumed to be quiescent at tquiescent (2021-04-09 02:00:00), and the direction to zenith at tquiescent and the position of the source in the DECam DIA subtraction are shown in blue.The maximum astrometric residuals (0.15") are indicated by the error bars on the DECam event positions.Arrows point to the zenith-bound direction at the time of the observation, showing the axis along which the next observation should be found if its motion were dominated by DCR.

Figure 12 .
Figure 12.Center: ∆DCR induced as a function of flare temperature and airmass in LSST g-band.Right: Box-andwhiskers plot of peak effective flare temperatures measured by Howard et al. (2020).Top: histogram of the per-visit airmass in the current LSST baseline observing strategy (baseline v3.0 10yrs, Bianco & The Rubin Observatory Survey Cadence Optimization Committee 2023).The median airmass and temperature are indicated by orange lines.

Figure 14 .
Figure 14.Airmass sky distribution of the current LSST survey strategy proposal (baseline v3.0 10yrs, Bianco & The Rubin Observatory Survey Cadence Optimization Committee 2023): skymaps are shown for the median airmass (left) and maximum airmass (right) in g-band.Plots produced within the Rubin Metric Analysis Framework.(Jones et al. 2014).To produce these maps, a sky segmentation into healpixels with resolution 64(Górski et al. 2005)

Figure 15 .
Figure 15.Top: Distribution of observation airmasses over five LSST recent observing strategy proposals (as labeled).While it remains small, the planned fraction of high-airmass observations in LSST has grown over time.Bottom: the corresponding probability of measuring a ∆DCR: as Figure13but fixing q = 0.05 and for the same five simulations of the LSST as in the top panel.

Figure 16 .
Figure 16.Simulation of the impact of three possible DCR correction strategies applied in g-band DIA.The pixel scale of the LSST system (0.2 arcsec/pixel) is indicated by the grid and the zenith direction by an arrow.(a): a quiescent dM source observed at airmass X = 1.22.(b): The source is observed during a T = 10, 000 K flare at the same airmass.In both panels b and d, the flare flux is normalized such that the total flux at flare is 50% greater than the total flux in quiescence.The difference of the flare and quiescent images is shown in the (c); contours help the reader quantify the PSF width and centroid in the difference image.The quiescent component is subtracted cleanly in DIA and the residuals associated with the flare directly quantify ∆DCR.Matching the DCR templates with correct assumptions on the quiescent SED would replicate this result.In (d ), the same flare is observed at airmass X = 1.74.DCR shifts the centroid of the quiescent source in the zenith direction by ∼ 42 arcsec.This shift is corrected in the visualization.The flaring PSF center is offset from the center of the image by 0.21 arcsec due to ∆DCR, and the PSF is significantly elongated in the zenith direction (with PSF variance increase ∆V = 0.14 arcsec 2 ).(e): The quiescent and flare images are differenced after modeling the DCR effects on the quiescent source using the true dM SED.The quiescent components is removed cleanly and the residuals can be used to measure ∆DCR directly as in (c).(f ): For the flare observed at airmass X = 1.74 (d ), the same DCR correction as for (e) are applied, but based on a G2V (solar type) quiescent SED: the flare and quiescent images are differenced resulting in non-zero residuals for the quiescent component and distorted residuals for the flare component.

Table 1 .
Candidate ID, Coordinates, timestamps, and g magnitudes with errors, as reported by the DIA pipeline described in subsection 4.2, for each DIA detection of each candidate circled in Figure8and DC21engi.The coordinates, detection time, magnitude, and component of the motion in the direction of zenith as measured between the first and each successive observation are indicated.|| -vs-∆g and see if we can detect a ∆DCR.To do this, we crossmatched the DECam DDF source location with a sample of known cool stars in the COSMOS field from the Gaia DR3 catalog