Distance Estimate Method for Asymptotic Giant Branch Stars Using Infrared Spectral Energy Distributions

We present a method to estimate distances to asymptotic giant branch (AGB) stars in the Galaxy, using spectral energy distributions (SEDs) in the near- and mid-infrared. By assuming that a given set of source properties (initial mass, stellar temperature, composition, and evolutionary stage) will provide a typical SED shape and brightness, sources are color matched to a distance-calibrated template and thereafter scaled to extract the distance. The method is tested by comparing the distances obtained to those estimated from very long baseline interferometry or Gaia parallax measurements, yielding a strong correlation in both cases. Additional templates are formed by constructing a source sample likely to be close to the Galactic center, and thus with a common, typical distance for calibration of the templates. These first results provide statistical distance estimates to a set of almost 15,000 Milky Way AGB stars belonging to the Bulge Asymmetries and Dynamical Evolution (BAaDE) survey, with typical distance errors of ±35%. With these statistical distances, a map of the intermediate-age population of stars traced by AGBs is formed, and a clear bar structure can be discerned, consistent with the previously reported inclination angle of 30° to the GC–Sun direction vector. These results motivate deeper studies of the AGB population to tease out the intermediate-age stellar distribution throughout the Galaxy, as well as determining statistical properties of the AGB population luminosity and mass-loss-rate distributions.


INTRODUCTION
Elements of Galactic structure are largely derived from multi-wavelength observations and models of stars and gas in the Milky Way as well as through comparison with extragalactic systems.The Milky Way is often modeled as an asymmetric bulge observed in the near-infrared (near-IR) (Blitz & Spergel 1991), and a logarithmic structure of the spiral arms (Hou & Han 2014;Quiroga-Nuñez et al. 2017;Reid 2022).Individual stars in the bulge are difficult to map due to the very high extinction values hindering even near-IR observations.Specifically, in the Galactic Center (GC), A V can be as high as 90 magnitudes (Elmegreen et al. 2009).These regions are better studied using longer wavelengths as * Adjunct Astronomer at the National Radio Astronomy Observatory † Adjunct Professor at the Department of Physics and Astronomy, the University of New Mexico interstellar extinction is an inversely dependent function of the wavelength.The Bulge Asymmetries and Dynamical Evolution (BAaDE) radio-wavelength survey aims to present a comprehensive study of the inner regions of the Galaxy to improve our understanding of Galactic structure and dynamics, with a focus on the bulge stellar population distribution and age (Sjouwerman et al. 2017;Lewis et al. 2020;Sjouwerman et al. 2024).The BAaDE survey consists of 28,062 infrared color-selected red giant stars, the majority of which are of Mira-type and lie on the Asymptotic Giant Branch (AGB).Approximately 10,000 of these stars have measured line-ofsight velocities determined from SiO maser lines (Stroh et al. 2019;Lewis 2021).In order to optimize how these velocities are incorporated into dynamical models, and to allow any existing spatial separation between populations to be distinguished, a 6D phase-space (position-velocity) is ultimately desirable.With distance estimates, the single epoch BAaDE survey can provide 3D positions along with a line-of-sight velocity.Distance estimates further enable determination of intrinsic AGB stellar properties like luminosity, mass-loss rate and SiO maser luminosities.Mapping the intermediateage stellar population in the bulge complements the BeSSeL survey's delineation of young stars in the disk and spiral arms (Reid 2022).
In recent years, the Gaia satellite has been instrumental in determining parallaxes to a large number of stars in the Milky Way (Vallenari et al. 2023).However, the vast majority of the BAaDE AGB stars lack reliable Gaia parallaxes (e.g., Xu et al. 2019;Van Langevelde et al. 2018, and references therein), therefore alternative methods to determine distances to AGB stars must be explored (Quiroga-Nunez et al. 2022).Very Long Baseline Interferometry (VLBI) parallax measurements would provide another pathway, and successful VLBI parallax distances to AGB stars obtained using OH maser lines at 1.6 GHz have indeed been reported (Van Langevelde et al. 2003).However, measuring VLBI parallaxes for a couple of thousand sources would be excruciatingly time-consuming and is also technically challenging to perform at the frequencies of the SiO maser (43 and 86 GHz).
Distances to AGB stars have also been estimated through the phase-lag method, which relies on comparing the angular stellar size to the absolute one.The absolute AGB stellar size is obtained by considering the lag time between the variations in the stellar light and variations in either the dust-scattered light in the circumstellar envelope (CSE) or OH maser emission (Etoka et al. 2017).Similar to the VLBI parallaxes, applying this method to a large sample of AGB stars would be very time consuming, as this requires regular flux measurements as well as accurate determination of the angular sizes (Maercker et al. 2018).
Finally, a commonly explored method for variable stars is using a known Period-Luminosity (P-L) relation.For Milky Way Miras this is hampered by the lack of a well-defined P-L relation.P-L relations for the less metal-rich AGBs in the Large Magellanic Cloud (LMC) have been derived (Whitelock et al. 2008), but work is still ongoing to better define the relation within the Milky Way, including effects of the circumstellar envelope (Lewis et al. 2023).
Due to the sizeable AGB sample in our survey, we aim to explore a method which can be consistently applied to any AGB star within the full sample.In this paper, we discuss an approach using distance-calibrated IR Spectral Energy Distribution (SED) templates.Ancillary photometric data are obtained from sky surveys ranging from the optical to the far-IR.The proposed method is advantageous as it builds on utilization of existing infrared catalogs, and can be used for AGB stars throughout the Galaxy without necessitating new observations.The methodology is outlined in Sect.2, with the results from testing the method given in Sect.3. In Sect.
4 the results are applied to consider the 3D distribution of the BAaDE AGBs in the bulge region.

METHODOLOGY
Our method is based on scaling measured flux densities in multiple IR filters to match the flux densities of a template with a known distance.In simple terms, it is a variation of the standard candle technique.However, our sources are not typically considered standard candles since they are located over a broad range of luminosities on the Hertzsprung-Russell diagram.Moreover, they exhibit strong variability, which introduces uncertainties when using single-epoch photometry.
In order to address these limitations, we have expanded the standard candle technique by constructing template Spectral Energy Distributions (SEDs) and categorizing sources based on their SED shapes.This presupposes that sources with similar underlying properties (initial mass, stellar temperature, composition, evolutionary stage) will exhibit similar SED shapes1 and luminosities, and that the shapes will vary sufficiently for various values of those properties, allowing for a unique shape to be discerned for each value of the luminosity.This is supported by previous work on Mira-variables both in the Milky Way as well as in the LMC, demonstrating a dependency of the absolute magnitude and colors in both the near-IR and mid-IR for O-rich AGBs (Guandalini & Busso 2008;Glass et al. 2009;Lebzelter et al. 2018;Smith 2022).
Finally, we employ more than one photometry point (an entire 5-7 point SED) to minimize uncertainties introduced by variability.While this method may yield substantial uncertainties for individual sources, statistically we will be able to use the distances to discern the 3D distribution of the targets.We will also be able to infer (in a forthcoming paper), for example, the luminosity and mass-loss rate distributions within our sample.In this section we describe the methodology to derive a distance based on the SED shapes (Sect.2.1) and how variability uncertainties are folded in (Sect.2.2), followed by a description of interstellar extinction corrections (Sect.2.3) and the estimate of the distance error (Sect.2.4).

SED shape categorization and distance extraction
The dependence of AGB SEDs on stellar and dust parameters has been addressed previously by modeling (e.g., Groenewegen 2006;Ventura et al. 2013;Dell'Agli et al. 2015;Jiménez-Esteban & Engels 2015).AGB stars generally possess dust-containing CSEs due to the significant mass-loss taking place in the form of stellar winds.Excluding interstellar extinction, an AGB SED can be approximated to have shorter-wavelength (optical and near-IR) radiation from the central star (which can be significantly absorbed by dust in the CSE) and mid-IR emission from the dust in the CSE.In order to categorize the SEDs we use three different colors: Figure 1 shows the SEDs for two well-known AGB stars, S Crt and OZ Gem, illustrating that the color differences between two AGB objects can be significant.For targets and/ or templates where MSX data is not available, AKARI data at 9µm and 18µm are used instead, forming the By constructing distance-calibrated SED templates for a wide range of the three colors, a source falling within a certain SED shape category can then have its distance extracted at any wavelength: where d tmpl is the template distance, F λ,tmpl and F λ,tgt are the flux density of the template and target, respectively.A K s is the extinction at K s -band and Z λ describes the extinction curve.A final distance estimate is taken as the median of the distance estimates calculated at each wavelength, to minimize variability uncertainties.If the source had an even number of data points, the median was defined as the mean of the two middle points.We note that we experimented with using an average, applying various weighting functions based on, for instance, Z λ as extinction effects are more significant at the shorter wavelengths.However, due to the limited set of wavebands < 10, the average is susceptible to offset data points significantly skewing the distance value, and we therefore opted to use the median.
The color-matching method used in this work relies on mass-loss rate as the determinant factor of the CSE mid-IR colors used ( ).This assumption was leveraged in choosing the initial BAaDE sample, enabling the formation of a set of AGBs with probable SiO emission in their CSEs (Sjouwerman et al. 2009).These objects are also inclined towards being Mira-type evolved stars, with expected mass-loss rates from a few 10 −7 to a few 10 −5 M ⊙ yr −1 (Whitelock et al. 1994;Le Bertre & Winters 1998;Höfner & Olofsson 2018).In order to confirm the assumption of mass loss dominating the mid-IR colors, we calculated the mid-IR [9] − [18] colors resulting from a stellar photosphere with T eff = 3, 300 K as a function of the mass-loss rate.The photospheric models were adopted from Gustafsson et al. (2008).These calculations showed that already at small mass-loss rates of ∼ 5 × 10 −7 M ⊙ yr −1 the mid-IR colors are dominated by the mass loss.Increasing the mass-loss rate by a factor of 10 causes a redder color of 0.5 − 1 mag (not a linear relation), measurable in the color range of our data set where the mid-IR colors vary with approximately 1.5 mag.Varying T eff does not significantly affect the resulting mid-IR colors in the model as mid-IR is within the Rayleigh-Jeans limit of a blackbody spectrum.

Variability-induced uncertainties
For each of the three colors used, the color range is split up into "bins" within which we consider the colors to be comparable.In order to determine the width of the bins representative of a given template, we note we are limited by the strong source variability and the associated uncertainty from using primarily single-epoch archival data from MSX, AKARI, and limited-epoch data from 2MASS catalogs.Using color bin widths significantly narrower than the variability in a certain band will not provide any benefit.
To estimate the typical source variability and how it changes with wavelength, we utilize existing IR catalogs from instruments equipped with comparable filters: 2MASS J and DENIS J, 2MASS K s and DENIS K s , and MSX A and AKARI 9 bands.The various surveys are not observed simultaneously, so comparing the source magnitudes between the surveys provide a statistical estimate of typical variabilities.Figure 2 shows the distribution of the magnitude differences for these three comparable filter sets, and the resulting 1σ widths (derived from the width containing 68% of the sources under the assumption of a symmetric distribution) are 0.61, 0.52, and 0.27 mag, respectively.For this paper, we therefore adopt a typical bin width of ±0.25 mag around the bin center value.Table 1 lists the uncertainty estimates for the three wavelengths plus the values interpolated to other wavebands used in the survey, which will represent the 1σ source variability as a function of wavelength for our sample.

Interstellar extinction corrections
As most of our sources are in the Galactic plane and we use near-IR photometric magnitudes, an interstellar extinction correction is required before a target can be categorized and matched with a template.There are multiple interstellar extinction maps available in the literature and most of them are limited to certain regions of the sky and/or have limitations in the assumed distance to the target.Our objective was to use an extinction correction method that could be applied consistently to the whole sample across the Galactic plane without any distance assumption.

SED color-excess extinction estimates
One of the traditional methods to determine the extinction A λ involves taking the ratio of the spectral flux density of an obscured reddened target, F λ , with that of a de-reddened, extinction free star roughly of the same spectral class, F λ,0 (e.g., De Marchi et al. 2014).This is also known as the 'pair method', and requires the target and the de-reddened source to have the same spectral type, distance, and absolute luminosity.Satisfying all these conditions, specifically for our set of sources without distances, is difficult.Hence, instead of determining the extinction directly we calculate the color excess, which is a quantity independent of the distance.The color excess is defined as: where F ref,λ is the flux density of the reference, dereddened object, and F λ is the flux density of the target.Using 2MASS J and K s as λ i and λ j , respectively, A K s can subsequently be calculated through: where C JK s is a constant which depends on the slope of the power-law relating extinction and wavelength in the near-IR (A λ ∝ λ −α ).The value used for C JK s is 0.537, consistent with an α = 1.9 following the Cardelli et al. (1989) extinction law.
In our case, the F ref,λ will not be the underlying, extinctionfree star but instead refer to the SEDs provided from wellknown AGB sources with very little interstellar extinction.The F ref,λ may thus still contain reddening from the CSE.This matters little for the method (hereafter denoted 'SED color-excess extinction method') given that we compare an observed SED to to its reference SED (with minimal interstellar extinction) based on the two having the same SED color characteristics including the effects of both CSEs.We note that the result is an estimate of the interstellar extinction only, excluding the CSE extinction.
To select targets which have comparable SED shapes as that of the reference, mid-IR MSX Using S Crt and U Her as references, extinction estimates to 7,259 and 11,594 unique sources, respectively, were calculated.In addition, 1,878 sources were accessible to both references; for those we averaged the two extinction estimates.Negative extinction values were derived for 190 sources, and a closer inspection of those sources revealed that most are very bright in the near-IR, implying they are likely foreground sources and hence with little interstellar extinction.We opted to simply exclude these objects in the following calculations.In summary, we obtained positive extinctions to 20,541 sources.
The error in A K s includes the propagation of the source flux density uncertainty due to the source variability (see Sect. 2.1) and the reference SED flux density measurement errors.An additional error term needs to be included due to the lack of our knowledge of T eff for stars in our sample.Following Aringer et al. (2016), for a typical AGB star mass and radius with log(g) = 0 and solar abundance, changes in T eff between 2, 800 − 3, 800 K produces changes in the (intrinsic) [

Comparison to other extinction estimates
The resulting A K s values extracted through the SED colorexcess method can be compared to those obtained by other methods, for example the 2D extinction maps published by Nidever et al. (2012), and Gonzalez et al. (2012).The BEAM calculator is based on color excess using Red Clump (RC) stars along a given line-of-sight (Gonzalez et al. 2012), and illustrate the typical differences found by the color-matched SED method and 2D extinction maps.Figure 3 shows the color-matched SED A K s values versus the BEAM calculator values for the BAaDE sample.The data is divided up in two brightness sets, where the top panel are foreground sources selected by considering the uncorrected K s band magnitude where K s < 6 mag (for l < 0) and K s < 5.5 (for l > 0) represents all the foreground sources, using the results from Trapp et al. (2018).The remaining, fainter sources constitute the bulge sample, plotted in the bottom panel.For the foreground sample our SED color-excess method produces less extinction compared to BEAM, which is reasonable since BEAM is focused on the extinction in the bulge where RC stars are numerous.In the bulge sample, A K s values for sources close to the plane agree well between the two methods, but the SED color-excess method derives additional extinction for sources at latitudes |b| ≳ 3 • .This may be due to the RC stars being biased toward the near side of the bulge, and/or several BAaDE targets originating from the far side of the bulge.BEAM has a quoted uncertainty on the extinction values which ranges from 0.1 mag when |b| ≈ 4 • up to 0.35 mag closer to the plane.
We further compared our A K s values to those achieved by following Lewis et al. (2023) and Messineo et al. (2005) who applied a 2MASS Red Giant Branch (RGB) matching method.The SED color-excess method produces extinction values which are slightly larger compared to the RGB method for sources likely in the bulge and close to the plane, and less extinction for more nearby targets.Messineo et al. (2005) quotes an uncertainty on their extinction values which is equal to 0.1 mag when A K s ≈ 0.6 mag and can range up to 0.7 mag for higher extinction values.These comparisons illustrate the effectiveness of the SED color-excess method for our AGB sample as we do not need to apply an a-priori known distance.

Distance error estimates
The error in the distance estimates obtained through our method are likely large due to the non-standard candle nature of the sources in addition to their large variability in the infrared.A minimum error at each wavelength can be estimated using error propagation with the error coming both from the template distance error, ∆d tmpl , and from the distance scaling in Eq. 1: ∆F tmpl and ∆F tgt are the flux density errors of the template and target, respectively.The last term includes uncertainties propagated from the extinction correction (Eq.1): (5) where ∆F tgt,raw denotes the uncertainty in the raw target flux density before extinction correction, ∆Z λ is the extinction curve uncertainty, and ∆A K s the assigned extinction For the bulge sample in the bottom panel the two methods agree well for sources in the plane, and deviate for sources at slightly higher absolute latitudes.
value uncertainty.For ∆A K s a typical value of 0.5 mag for S Crt and 0.47 mag for U Her is assumed, derived from the extinction estimates described in Sect.2.3.The ∆Z λ uncertainty is harder to assess, and we therefore will set this value to 0 for simplicity.∆F tgt,raw contains the observational photometric error and the source variability error.Note that the photometric error is much smaller (for example, a few percent for 2MASS) than the error due to the source variability which can be ≈ 1 − 2 magnitudes in the near-IR, falling off to < 1 magnitude in the mid-IR.While we do not have lightcurves for all targets, we instead apply the typical variability error to each source depending on the waveband (Table 1).The total distance error is then reduced by  error achieved by the above, this provides us with a rough typical error for the source distances at a given wavelength.

DISTANCE ESTIMATES
With the interstellar extinction corrections performed, the methodology for distance estimates outlined in Sect.2.1 can now be applied.Two different types of SED templates are used to cover a broad range of target colors.The first type is based on nearby, well-studied AGBs with VLBI parallax measurements available for the distance calibration (Sect.3.1) and to which the SED distances can be compared.In order to extend the accessible color range, the second type of templates is formed by using a group of targets assumed to have a mean distance of that of the GC, thereby providing the template distance calibration (Sect.3.2).

VLBI parallax calibrated templates
To validate the method, a sample of AGB stars with VLBI parallax measurements were identified in the literature.We constrained the sample to 14 stars for which IR data could be consistently be collected from Vizier via a name search, creating an IR data set with flux densities between 1 and 18µm from 2MASS, WISE and AKARI (Skrutskie, M. F. 2003; AKARI Team 2020; Wright, Edward L. 2019).The sample of stars, along with their parallax measurements, is listed in Table 2.These VLBI sources are located well above the plane, all with |b| > 20 • , and their VLBI parallax measure- ments show they are all located closer than 1.3 kpc.Interstellar extinction was therefore ignored for this sample.All the 14 objects have J and K s band photometry as well as the AKARI 9µm and 18µm (only 1 of the 14 sources have MSX data, thus AKARI 9µm and 18µm was used instead).These 14 objects cover a [J] − [K s ] color range of 1.90 mag, an AKARI [9] − [18] color range of 0.98 mag, and a [K s ] − [A] color range of 3.09 mag.
Various methods of forming a template were tested, including a median of the values across sources with similar SEDs, and individual object SEDs.It was found that working with an individual object's SED as a template worked well as long as the template object was not an outlier in terms of its colors compared to the other sources.This is in agreement with our assumption that the SED shape matching is of essence.Figure 4 shows the comparison of VLBI parallax to SED distance estimates using T Lep and R Uma as templates.For this plot, three VLBI sources fell within ±0.25 mag of the colors of T Lep and four VLBI sources fell within ±0.25 mag for R Uma.A Pearson correlation coefficient of 0.96 and 0.98 is achieved for the T Lep and R Uma sources, respectively.
In the next step we apply the method to the BAaDE sources with colors matching the template, and then use Gaia parallaxes to test the derived distances.T Lep was selected as the template due to the largest possible overlap in the BAaDE sample color regimes, allowing for a large number of color-matched BAaDE sources to which distances can be estimated.4,208 BAaDE AGB sources with A K s values (Sect.2.3) fall within the color range with ±0.25 mag of those of T Lep, and for which we calculated SED distances.Distance error estimates follow Sect.2.4 with a few modifications.First, as T Lep is a single source which is well studied, the template flux density variability applied was taken from the NASA/IPAC Infrared Science Archive (IRSA), where the uncertainty value has been estimated over a number of observations at each wavelength and the near-IR photometric errors ranged from 0.26-0.35mag.Second, T Lep is a nearby source (∼327 pc; Table 2) and no interstellar extinction was applied.With these assumptions the relative distance uncertainties ranged between ±33% − ±52%.The error variation is primarily driven by the number of wavelength points used, which varies between four and seven for the sources matched to the VLBI template.
A comparison distance data set was constructed through a cross-match to the Gaia DR3 database.Out of the 4,208 targets, 541 have associated Gaia DR3 parallaxes with parallax errors < 20 % and for which we can expect to derive parallax distances straightforwardly (Bailer-Jones 2015).However, AGB stars are prone to large errors in Gaia astrometry measurements, due to their obscuration, large variability, and extended size (Xu et al. 2019;Van Langevelde et al. 2018).Andriantsaralaza et al. (2022) points out difficulties for Gaia in providing reliable parallaxes for AGB stars, even if parallax errors are limited to < 20%.They further note that the relative parallax errors must be corrected with a factor depending on the G magnitude, with the largest error inflation factor for the brightest objects (G < 8 mag).Following their work to correct the parallax errors for our objects, we plot the SED distances versus the Gaia parallaxes in the right hand side panel of Fig. 4. Considering the Pearson correlation coefficient excluding 32 sources outside 2σ of the 1-1 correlation, (since the Pearson correlation test is sensitive to outliers) we obtain a moderately negative correlation with the Gaia parallaxes, which increases to a strong negative correla- tion for the brighter objects.For the brightest sources (G < 8) the Pearson correlation coefficient is -0.83, for sources with 8 < G < 12 it is -0.65, and for the faintest sources with G > 12 the correlation is the lowest at -0.35.A systematic difference between the SED distances and the Gaia parallax distance is apparent once the SED distances exceed 2000 pc, where the SED distances are consistently larger than the Gaia parallax distance.This portion of the plot consists of an increasingly larger number of fainter sources (G > 12), and the sources with corrected parallax errors larger than 20%.In this portion, the average deviation between the SED distance and the Gaia distance is approximately 1 kpc.For distances beyond about 1-2 kpc Andriantsaralaza et al. (2022) observe a similar offset between the inverse of the Gaia parallaxes and distances derived using priors once Gaia parallax errors exceed 18%, indicating the inverse of the Gaia parallaxes are systematically underestimating the distances.
From the comparison with VLBI and Gaia parallaxes, we conclude that the SED template method of estimating distances to AGB stars works well statistically when using a VLBI parallax calibrated template.

Inner Galaxy distance-calibrated templates
Using AGB sources with VLBI parallaxes as SED templates provided promising results.However, the VLBI tar- gets are all nearby and relatively blue compared to the BAaDE sample as a whole (Fig. 5), which necessitates constructing additional templates.We use the sub-sample defined in Lewis et al. (2023), which is likely to have a mean distance of the distance to the Galactic center of 8.277 kpc (Abuter et al. 2022).Their sample is based on selecting BAaDE sources with Galactic longitude |l| < 3 • and Galactic latitude |b| < 4 • which have absolute SiO line-of-sight velocities > 100 km s −1 .The coordinate selection is implemented to select the sources within the central galactic region and the velocity cut is applied to separate out foreground from background stars (Lewis et al. 2023).This GC sample of 518 sources, being a subset of the BAaDE sample, has MSX data which allows using the [A] − [D] color instead of the AKARI colors used in Sect.3.1.With the MSX photometry an additional four filters are included (MSX A , MSX C , MSX D, MSX E), further driving down the uncertainty due to variability (MSX Team 2019).
The newly formed GC sample was first divided into nine equally sized [K s ] − [A] color bins (GC1-GC9) as the [K s ] − [A] color spanned the largest color range of the three colors (0.74−7.45 mag).From each bin a template was constructed by using the median color and median flux density.Using the median instead of the mean is appropriate as our flux density distributions are not Gaussian.Fig. 5 shows the position of the GC templates as well as that of T Lep in the color-color regimes, respectively.Table 3 lists the template mean colors, the number of sources used to construct a given template, and how many target sources were color-matched to each template using ±0.25 magnitudes around each template color center.By using the GC sample for templates, we thus obtained distances to 10,446 color-matched BAaDE sources (Table 4).Errors were estimated for the templates and targets following Sect.2.4.In contrast to the VLBI template case which used a single source (Sect.3.1), the GC templates were formed using N data points at a given wavelength, and the template flux density errors are driven down according to

√
N. This results in template errors which are small compared to the target flux density uncertainties now dominating the error.Consequently, typical distance errors range between ±27% to ±41%.The variation of the errors depends primarily on the number of wavebands that are being used, ranging between 4 and 11.The smaller errors obtained using GC templates compared to the VLBI template reflects the improvement gained by using the median from a large set of sources rather than using a single source SED.The VLBI sample, however, is too small to allow for this statistical approach to forming a template.
For the GC color-matched sample there is no independent set of parallax distances to compare to, and instead the SED color-matched distances are compared to a set of distances derived from Period-Luminosity (P-L) relations.First, we compared our distances with the P-L distance estimates from OGLE periods derived using LMC Miras (Iwanek et al. 2023).3,553 sources with distances are cross-matched between their sample and ours.The top panel of Fig. 6 shows the LMC P-L distances versus the SED distances, and the SED distances are systematically larger than the P-L distances.A couple of reasons causing this can be considered; Iwanek et al. (2023) obtained a distance to the GC of 7.66 kpc (implying the P-L distances are short), and also applied a P-L luminosity relation from the LMC which has not yet been proven to hold for the more metal rich Miras in the Milky Way.We further note that the P-L relations used were derived from near-IR photometry, which is affected by CSE extinction.Our sources, being chosen to be likely maserbearing, probably have more substantial envelopes than the LMC sources from which the LMC relation was made.
P-L relations for very red Mira variables are indeed uncertain with a large scatter, but Lewis et al. (2023) have shown that P-L relations in the mid-IR show a tighter correlation than in the near-IR for maser-bearing Galactic Miras.For our sample P-L distances were calculated for MSX A, C,  2023) sample derived from P-L relations for Miras in the LMC.The distances systematically deviate from the SED method which yields larger distances.Bottom: Comparison of the SED distances for the 3,722 cross-matched sources between our sample and an OGLE sample for which a P-L relation in the mid-IR for Milky Way Miras has been applied (Lewis et al. 2023), showing an improved correlation.and D bands separately using the relations given in Lewis et al. (2023) and the interstellar extinction methods described in this work.The mean distance was then applied to each source, and compared to the SED color-matched distance (bottom panel in Fig. 6).In comparing the MSX P-L distance from Lewis et al. (2023) to the SED-derived distance, we find a stronger correlation than for the distances from Iwanek et al. (2023).The systematic deviation is largely removed when using a mid-IR P-L relation derived for Milky Way AGBs.On average, the LMC P-L distances are 10% shorter than the SED distances, whereas the MSX P-L distances are only 2% larger than the SED ones.This is not completely surprising, given that the near-IR relations were derived from low-metallicity LMC objects.There is also scatter due to  3. The Milky Way disk with a radius of 14 kpc centered at the GC at (0, 0) is indicated with a solid circle, and the position of the Sun is indicated with a star at (0, 8.227).The finger-like structures are regions of deeper scans in the MSX survey.An outline of the bulge is shown with the ellipse centered on (0,0) with a semi-major axis of 4 kpc and a semi-minor axis of 2.2 kpc oriented at an angle of 30 • clockwise from the GC-Sun direction vector.
variability in the different bands (including the OGLE data), but variability amplitude and extinction is smaller in the mid-IR compared to the near-IR.

DISCUSSION
Using the SED color-matching technique, distance estimates to a set of 14,654 AGB sources have been derived.In the following discussion related to bulge structure, we largely focus on the 10,446 reddest targets color-matched to the GC templates, as this sample contains most of the bulge objects.

Overall distance distribution
The resulting distances can be assessed by constructing a 2D Cartesian Milky Way face-on plot which provides a visualization of the resulting source distribution for each of the nine GC templates separately (Fig. 7).For GC4-9, there is no obvious mean distance shift of redder sources toward farther distances, indicating interstellar extinction corrections applied were reasonable.For the three bluer templates, GC1-3, the mean distances are shifted to smaller distances, which may indicate a problem with interstellar extinction corrections, However, this shift can be understood if considering the source K s versus [J] − [K s ] magnitude-color diagram, where most of the bluer targets are also the brighter ones and thus likely to be foreground disk sources (Trapp et al. 2018).

Discerning the bar structure
If the bulge contains a bar structure whose long axis has an inclination of 30 • to the GC-Sun direction vector, we expect that the distance distributions toward the near and far sides of the bar structure should show an offset between their peaks.In order to probe this expectation, we constructed separate histograms of the distances for sources between 8 • < l < 16 • (near side) and −12 • < l < −4 • (far side).Note that we are not selecting the two cuts symmetrically in longitude, but such that we are likely to intersect the bar at at a similar radius from the GC.We then removed foreground sources in each sample by considering the uncorrected K s band magnitude where K s < 6 mag (for l < 0) and K s < 5.5 (for l > 0) represents all the foreground sources, using the results from Trapp et al. (2018).After applying these magnitude cuts, Fig. 8 demonstrates a clear distinction between the distributions, indicative of the presence of a stellar bar.These distributions can be compared to simple stellar density distribution models of the bulge/bar.The first model is an elliptical bar-like distribution with an inclination angle of 30 • and with an exponential fall off in the number of sources in the radial direction away from the GC.This model shows two distinct peaks similar to our data.In the second model we use a spherical distribution with an exponential fall off, demonstrating that no peak offset would be observed for a symmetric stellar distribution.We note that the aim of the models is not to match the data perfectly, but rather to illustrate that a spherically symmetric stellar distribution does not result in distance distribution offset between the positive and negative longitude cuts.More detailed modeling work will be needed to fully represent our data.By separating sources likely to be in the disk from those in the bulge via the K s magnitude cut, Fig. 9 presents the resulting BAaDE AGB source density distributions for the bulge and foreground disk, respectively.The bar structure is clearly discernible and consistent with the inclination angle of 30 • to the GC-Sun direction vector reported on by other groups (e.g., Wegg et al. 2015).A proper fit of the bar inclination angle from our data at this point is difficult due to two issues.First, the far side of the bar contains fewer sources, demonstrating that the far side may be affected by observational selection effects such as sensitivity and source confusion arising from the MSX catalog from which the BAaDE sample originated.The near side of the bar is more pronounced and may be less affected by such observational effects.Second, the extension to the GC direction beyond the GC is the result of a deeper scan MSX performed in the GC area (|l| ≤ 0.5 • ), representing an uneven sampling depth.We refrain from making conclusions about the resulting bar angle until these limitations in our sample have been resolved.
The foreground sources belonging to the disk also contains a selection bias, as the BAaDE survey was constructed from MSX sources at latitudes |b| < 6 • .With a disk scale height of 300 pc, the survey does not reach the full scale height until a distance of around 3.4 kpc, reflected in the disk source density plot.This is consistent with the findings of Quiroga-Nuñez et al. (2020), finding a similar distribution for a selection of BAaDE targets within 2 kpc distance from the Sun.

CONCLUSIONS AND FUTURE WORK
By using infrared survey data, we have tested a method of estimating distances to AGB stars.The method relies on color-matching targets to distance-calibrated SED templates.The method allows for interstellar extinction corrections prior to extracting the distances.The results show distances which are consistent with VLBI and Gaia parallaxderived distances, as well as with distances determined from P-L relations derived in the mid-IR.Typical distance errors are estimated to ±35%, and we note that using a large set of sources to form the templates aids in driving down the total errors in the method.
The method was applied to the BAaDE sample and provided distances to almost 15,000 AGB stars.By mapping the sources we find that the intermediate-age AGB population traces the bar structure.However, to more reliably model the full bar structure including the far side a more uniform star coverage throughout the bulge region is needed.Future work therefore includes employing Machine-Learning techniques to fold in additional observed properties of the stars, including periods, to access targets which may not have all IR filters required for the color-matching applied in this paper.Other IR catalogs with deeper photometry than MSX (e.g., the AKARI and GLIMPSE catalogs), will also be folded into our methods in the future.The distance estimates will aid in determining luminosity and mass-loss rate distributions for Galactic AGB stars, and will be presented in a forthcoming paper.
Y.P. and R.B. acknowledge support from the the National Aeronautics and Space Administration (NASA) under grant number 80NSSC22K0482 issued through the NNH21ZDA001N Astrophysics Data Analysis Program (ADAP).R.S.'s contribution to the research described here was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with NASA, and funded in part by NASA via ADAP award number 80NM0018F0610.The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.
235µm and 2.159µm), [A] − [D] (MSX 8.28µm and 14.65µm), and [K s ] − [A].The [A] − [D] color thus informs primarily about the CSE properties, while the shorter wavelength [J] − [K s ] informs about the central star emission modified by the CSE.The [K s ] − [A] color represents the boundary region between the two components.
[9] − [18] color as a substitute for [A] − [D].This is justified by comparing the two colors for sources in our sample containing both AKARI and MSX data.The AKARI [9] − [18] color correlates with with the MSX [A] − [D] color with a Pearson correlation coefficient close to 0.7.

Figure 2 .
Figure2.Distribution of the magnitude differences observed between 2MASS J and DENIS J, 2MASS K s and DENIS K s , and MSX A and AKARI 9 bands, providing an estimate of typical source variability in these bands.The vertical lines show the regions within ±1σ assuming a symmetric distribution.
[A] − [D] colors were used as they are least affected by interstellar extinction 2 .For sources lacking MSX photometry, AKARI [9] − [18] was substituted for [A] − [D].Two reference sources were selected, S Crt and U Her, both with negligible interstellar extinction and with mid-IR colors falling within the color of the full BAaDE sample.S Crt and U Her have distinct AKARI [9] − [18] colors of 0.74 and 1.19, respectively, and by selecting sources with colors of ±0.25 magnitudes within those of the reference objects, approximately 74% of the BAaDE sample could be covered.
J] − [K] values ranging between 1.1 − 1.35.The corresponding changes in the [A] − [D] or [9] − [18] colors are significantly smaller.The resulting total A K s error is 0.5 mag and 0.48 mag for estimates based on the S Crt and U Her templates, respectively.

Figure 3 .
Figure 3. Extinction values for the BAaDE sample as estimated by the SED color excess method compared to the Gonzalez et al. (2012) BEAM calculator.The color scale represents the absolute value of the latitude of the target.The top panel shows a sample of likely foreground sources, for which the SED color excess method finds lower values consistent with them being nearby disk sources.For the bulge sample in the bottom panel the two methods agree well for sources in the plane, and deviate for sources at slightly higher absolute latitudes.

Figure 4 .
Figure 4. Left: The SED method distance estimates compared to independently measured VLBI parallaxes.Right: The SED distance estimates compared to Gaia parallax measurements.The blue symbols indicate the brightest Gaia sources with G < 8 and the red ones with 8 < G < 12.Both VLBI and Gaia parallax distances show a strong correlation with the SED distances.In both panels the black line denotes a 1-1 correlation.

Figure 5 .
Figure 5. Color-color diagram using [K s ] − [A], [J] − [K s ], and [A] − [D] colors.The distribution of the 20,541 BAaDE sources is indicated with the blue density contours, and the template bins are outlined with rectangles.The sources with VLBI parallaxes used for testing the method are plotted with plus-symbols.Note that the VLBI source template selection is confined to the bluest colors, while the GC templates access a broader range of the color-color space.

Figure 6 .
Figure 6.Top: Comparison of SED distances with 3,553 sources cross-matched with the Iwanek et al. (2023) sample derived from P-L relations for Miras in the LMC.The distances systematically deviate from the SED method which yields larger distances.Bottom: Comparison of the SED distances for the 3,722 cross-matched sources between our sample and an OGLE sample for which a P-L relation in the mid-IR for Milky Way Miras has been applied(Lewis et al. 2023), showing an improved correlation.

Figure 7
Figure 7. 2-D Cartesian plot to show the distance distribution for targets color-matched to the GC templates.Starting from the top left, each template corresponds to a redder [K s ]-[A] color according to Table3.The Milky Way disk with a radius of 14 kpc centered at the GC at (0, 0) is indicated with a solid circle, and the position of the Sun is indicated with a star at (0, 8.227).The finger-like structures are regions of deeper scans in the MSX survey.An outline of the bulge is shown with the ellipse centered on (0,0) with a semi-major axis of 4 kpc and a semi-minor axis of 2.2 kpc oriented at an angle of 30 • clockwise from the GC-Sun direction vector.

Figure 8 .
Figure 8. Top: Histogram showing SED distance distributions for the near (red) and far (blue) side of the bar with a clear shift between peaks of the two distributions.Middle: Distance distributions for the near and far side of an elliptical bar model, inclined at 30 • to the Sun-GC direction vector.Bottom: Distance distributions for a modeled spherical bulge.Note that the data (top panel) contains a contribution of foreground sources, lacking in the models used for the middle and bottom panels, explaining the much broader distributions present in the SED data.

Figure 9 .
Figure9.Top: Bulge AGB density distribution after removing sources most likely to belong to the disk through the magnitude cuts reported on inTrapp et al. (2018), see also the text.The Milky Way disk with a radius of 14 kpc centered at the GC at (0, 0) is indicated with a solid circle.An outline of the bulge is shown with the ellipse centered on (0,0) with a semi-major axis of 4 kpc and a semi-minor axis of 2.2 kpc oriented at an angle of 30 • clockwise from the GC-Sun direction vector.The near side of the bar is clearly discernible, while the far side lacks a comparable source coverage.Bottom: The distribution for the targets removed from the top panel, and which most likely are foreground disk sources.GC at (0, 0) is indicated with a solid circle and the position of the Sun is indicated with a star at (0, 8.277).The colorbar represents number of sources.
made use of the NASA/IPAC Infrared Science Archive, which is funded by the National Aeronautics and Space Administration and operated by the California Institute of Technology.This research made use of data products from the Midcourse Space Experiment.Processing of the data was funded by the Ballistic Missile Defense Organization with additional support from NASA Office of Space Science.This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation.This research is based on observations with AKARI, a JAXA project with the participation of ESA.

Table 1 .
Derived variability errors for the BAaDE sources as a function of wavelength.

Table 2 .
√N when N wavelength points are used.Although an individual source may have an error much larger (or smaller) than the total The 14 sources constituting our comparison sample with VLBI parallax measurements (ϖ VLBI

Table 3 .
The median color of the constructed templates for each of the three color regimes, the number of sources that went into each template (#Tmpl), and how many target sources could be colormatched with the template to estimate a distance.Note that the VLBI template used a single source, T Lep.

Table 4 .
SED distances to the 10,446 BAaDE sources matched to the GC templates including the estimated 1σ error.The full table is available electronically.