Possible next steps in XP spectra & the Milky Way

 These are literally only notes-to-self by HW

The Rich Heart of the Milky Way

[see also the Jupyter notebook: Metal Rich Stars from Andrae+2023]

I just made the plot using the Andrae+2023 giant table

We can then identify the most metal-rich stars in it [M/H]>0.47.

Their CMD's look "normal"; but iffy [M/H]-estimates and imprecise parallaxes

prevent clear age-dating.

We then look at the all-sky distribution, in different metallicity slices, which looks remarkable

 To do's:

• plot [M/H]_XP vs DR17 APOGEE - get new metallicities from Andy Casey (?)
• check whether there are any RV's (few?)
• get good distances (a la Hogg...) or from Xiangyu
• with good/assumed distances, say something about the ages
• build a model to constrain the level of rotation using only proper motions (that needs good distances). -- get proper motions!! (the notebook is set up for that)
• all this requires a clean-up/augmentation of the Jupyter notebook.

Giant Ages within 3 kpc

For giants with XP metallicities and very good parallaxes (5%?) we can get ages, as in the LMC, using HWR's sped-up code.   We should do that..  perhaps suggest to to Josh?

Calibration of the stellar mass scale on the lower MS

I have taken the Hwang et al 2023 code that calculates (and dynamically calibrates) M_*(absG,B-R) from the wide-binary dynamics and applied it to M-dwarfs. Compared to the what the Carmenes DR1 folks use, there is a serious offset.

Interpreting the Gaia astrometric orbit catalog: constraints on objects at 2 AU

 "Modelling" the Gaia DR3 Astrometric Orbits (a.k.a. learning astrophysics from complicated catalogs)

The Basic Idea

Gaia has produced many catalogs with constraints on interesting astrophysical phenomena. In particular, Gaia DR3 has produced a catalog of >130,000 "astrometric orbit solutions", where many-epoch observations of a sources have lead to an "orbit".  Specifically, the light-centroid position of the source on the sky can be (well) fit by a combination of parallax, proper motion and a Keplerian orbit.

This catalog (lovingly dubbed gaiadr3.nss_two_body_orbit) should be a treasure trove on questions on how stars deal with their angular momentum during formation (binarity, or planetary system), on stellar evolution (what if one of the stars in a binary dies to become a WD, NS or BH)?

But it turns out that it is not easy to go from "What's (not) in the catalog?" to astrophysical questions such as "How often do low-mass stars have even lower mass companions?"  "How common are giant gas planets at, say, 3 year periods, around stars of different mass?", "How common are brown dwarfs are companions?", or "How often are stars orbiting stellar-mass black holes?"

This  catalog<-->question  link is not easy for a variety of reasons. First, and foremost, no catalog contains everything. It usually (and sensibly) contains what the underlying experiment can detect or measure well. Second, acknowledging this incompleteness does not yet help. To learn from the catalog one must build a model that predicts the ensemble of catalog entries, and this model must incorporate the selection function (e.g. Rix et al 2022; and references therein for a recent exposition of this issue). The selection function specifies with what probability (often 0 or 1) an object of any properties (flux, temperature, distance, direction, etc..) would have ended up in the catalog or not. The selection function should be a function of variables that a) are "observables" and b) can be predicted by the model.

 For "literature" samples, the selection function is hopefully specified in terms that lend themselves to modelling; else it may need to be (painstakingly) (re-)constructed. By comparison, specifying a "model" is often easy. E.g. for the scienc e theme discussed above, the model could be given by a single number: what is the probability pJ that a solar-type star in the Galactic neighbourhood has (as dominant companion) a Jupiter-mass object in a 3-year orbit.

To get specific, let's make a plot of the catalog in which the Gaia mission speaks about (astrometric) binary orbits.   There are many things to be plotted about this sample, but the following Figure (showing all 134,000 sample members) cuts most to the chase.

The X-axis shows the inferred period, between 0.1 years and 10 years, with the prominent 'aliasing gap' at 1 year. The Y-axis shows the 'Astrometric Mass-Ratio Function' (AMRF), a clever quantify devised by Shahaf et al. ( 2019 ). It links the light-centroid motion (observable) to the motions of the primary and secondary. If the primary dominates, then AMRF becomes an approximate measure of the mass ratio, q. 

So, interesting companions are either at very high AMRF (black holes, neutron stars ??) or at very low values (planets, etc..).  And, the plot shows very little / nothing at extreme AMRF. What does this mean for the incidence of BH's and planets orbiting stars in the pertinent period range....?

That's the basic issue to be addressed.  

The Basic Mathematics

Let's take as (interesting) toy case that we want to constrain is the rate at which solar-type stars have Jupiter-mass objects MJ in ~2AU orbits. 

This is a special case of the following general model that constrains the probability "p" that a stars of a c ertain mass and metallicity has a planet of Mp in an orbit size "a" and eccentricity "e".

This model is now to be compared to the "data" (which is in essence the set of catalog entries in the "orbit catalog" above).

In addition to the light-centroid parameters ("lc"), we need the parallax (to convert milli-arcseconds to AU), the apparent magnitudes \vec{m} (which set the S/N; see below), and an estimate of the primary star's mass (which comes from other information), as Kepler's law needs that.  The Gaia catalog contains such information for N_*=134,000 stars (see plot above).

The nest (and hard and crucial) step is to understand under which circumstances a system might have ended up in the catalog. Why this is crucial, we'll show below; let's for now just go with it.

For Gaia DR3 the following conditions have to be fulfilled (Halbwachs+2022).

These are S/N conditions on the parallax (varpi), and on the quality of the light-centroid semi-major axis (a_0). There is also a criterion on the eccentricity uncertainty. For operational reasons all these criteria depend on period.

The art, and bain, is to express the selection function S correctly. This function S is the probability that a system with certain physical characteristics -- stellar, and planetary mass, distance, period, etc.. -- would have passed the selection criteria. We have to be able to go to any (even just envisioned) star of certain mass distances etc., and answer: if that star had a planet of certain mass in a certain orbit, would it have resulted in a catalog entry. 

Initial version of such an approach was worked out in El-Badry+2023, for the regime of massive dark companions); there it seemed that the parallax/S/N criterion was the most stringent. If we assume that the parallax uncertainty is (mostly) a function of parallax, brightness of the stars and semi-major axis of the orbit, then we can make the selection function a relatively simple function of "model-predictable observables" [More to say here].

This then leads to the following prediction of how many systems we should see in the catalog:

How to interpret this? To predict the number of catalog entries that imply a planet of mass M_p in an orbit of size a -- for a given model p(M_p,a) -- we need to sum over all stars (yes, ALL stars) and sum up their selection function.  Of course the sum over all stars is silly, as for almost all stars in the Milky Way the selection function probability is S==0. So, in practice, it boils down to counting all stars in the catalog for which S==1 (in a given dM_P x da ) and comparing that to the number of actual catalog entries (using the binomial statistic). 

The delightful aspect is that -- once the selection function is specified -- the exercise becomes "trivial".

Selection Function Nitty-Gritty

<here we'll eventually spell out how to determine and then apply a suitable selection function>

Icochrone fitting (todo list)

 Next steps to explore the simplified isochrone fitting code

Code

Lives here: https://github.com/HWRix/stellarAges

To-Dos

Code changes / updates / tests

• use finer grids at <1 Gyrs
• understand whether the original Teff really needs to be modified, and why
• understand why the isochrones > 13 Gyrs give "funny" results
• the CMD distribution of stars identified as 8 < logAge < 8.5 looks funny
• I set a floor to the precision of the photometry; was that bad?
• conceptually, we still need to work out how to make all this "mass weighted"
• consider ingesting a prior on each star's extinction

Diagnostic plots

• look at the CMD distribution of stars in a grid of (FE_H and fitlogAge); what looks OK?
• look at the "tensions" between input and output: are there systematic differences between the input FeH, Teff, G and the best fit predictions of the isochrones?  
  Note: this requires that we also calculate the likelihood-weighted mean of the isochrone predictions, not only of the isochrone's "physical parameters"
• 
  

The Galactic bar, seen in stars of different [M/H]

 Which populations are bar-like in the Milky Way?

These are notes from a set of conversation between HWR and Madeline Lucey, during her visit to MPIA May 2023, augmented by a very few feasibility plots.  Worth exploring as a project?

The Milky Way has a prominent, presumably long-lived stellar bar, which is interesting in itself and whose dynamics affect the stellar dynamics well inside and outside of it. The exact length and pattern speed is still under debate, as is the question of how much it has slowed down.

Stellar populations of different age, abundances and kinematics (hot/cold) will contribute differently to forming and maintaining the bar, and will respond differently to the bar as an input perturbation.

Therefore it would be good to map, at first geometrically, the bar -- or the bar region in the MW -- in different stellar populations. Age might be the ideal population marker but is hard to determine consistently and precisely for stars that are 3-8kpc away from the Sun.  But metallicities, [M/H], are now available for vast samples in the inner Galaxy (say, RGC<5kpc). So, it appears feasible to make a metallicity-dependent map of the central galaxy.

However, naive maps n(X,Y | [M/H]) are likely to fail (== look confusing) for two reasons:

• (patchy) dust extinction and complex selection function
• imprecise distances that blur the image

A [M/H]-dependent bar map with RC stars

Does it make sense to try to construct a [M/H] dependent bar map?

Red clump (RC) stars of known metallicity are ~10% standard candles in their infrared fluxes (say, W1).  This implied distance precision of 5% is better than the expected parallax-based distance precision for heliocentric distances >3kpc. So, with good (perfect) identification of a set of stars as RCs of a given [M/H] a 5% distance precision map could conceivably be constructed.

The basic idea -- floated by HWR -- is to take advantage of the fact that the large majority of stars with luminosities within 0.5mag of M_W1(RC) ~ -1.7 and an appropriate temperature or color range are RCs anyway; once we know only that, we know their distance to 5-10%.  So one could consider: pick a broad RC candidate box in M_W1 (needs only poor parallax), Teff, logg space; identify likely RCs; then assume we know their photometric distance to 5-10%; and make a map.

Why might this work?  Rene Andrae's catalog contains 1.7M RC candidates beyond 4kpc (parallax<0.25), of which only 8000 have parallax S/N > 20.

Issues

• How to best identify RC stars, using XP (and Teff, logg,M/H), parallax, B,G,R,JHK,W1 etc..? There are simple ways (make cuts) and fancy ways (train and classify)
• How to deal with dust etc..

On Metal-Poor Stars at the Heart of the Milky Way

 Metal-Poor Stars at the Heart of the Milky Way

Motivation and Goals

When exploring the oldest stars in the Milky Way, there are two major components that are well-studied:

• the old disk (aka the thick disk, or alpha-enhanced disk). Our recent exploration has found that we can trace stars back to ages of 13 Gyrs, but that the [M/H] distribution extends from [M/H] ~ -1 for the oldest stars, to [M/H] ~ 0 for the youngest.  Which stars are responsible to enhance the gas of the old disk to [M/H] = -1 is not known:
• the stellar halo, which has been parsed by recent Gaia- and/or spectroscopy enabled analyses into a set of components of components that reflect the ancient merger history.  Their age - [M/H] distribution and abundance patterns generally imply an ex situ origin. The most prominent component dominating the inner halo (<15 kpc) is the GSE ("Gaia-Sausage-Enceladus"), stars on radial orbits with apo-centers ranging from 10kpc to 100 kpc, GSE stars reflect a merger with the proto-Milky Way about 11 kpc ago.
• 
  

  

  Age-[Fe/H] distribution from Xiang&Rix2022, showing that in samples at the soral radius, the old disk stats at 13 Gyrs and [Fe/H] ~ -1; the (ex-situ?) halo stars have a different enrichment track.

  
  

However, we know little about ancient, metal poor stars (say, [M/H] < -1) at the heart of the Milky Way, say, within 1-5 kpc of the Galactic Center. We should expect stars in the heart of the MW to be a mix of

• in-situ stars that reflect this earliest star formation in the MW's main potential well. Why "earliest"? Because it presumably did not take long at the heart of the MW to enrich birth gas to [M/H] >= -1.  Such stars are very much expected to exist from cosmological simulations. Whether they would be best dubbed as "innermost halo", or as "innermost old/thick disk" would depend on their angular momentum and terminology preferences.
• stars that spend much of their life at far greater Galactocentric radii, but are currently passing through the center. I.e. we should expect GSE stars near their peri-centers.

Selection Approach

Much of past research has focussed on finding halo stars away from the Galactic disk, as this mitigates the detrimental effects of dust extinction and of the sample contamination by more metal-rich stars. [Notable exceptions: early efforts with the sky mapper.]

Here we want to explore whether the Gaia DR3 BP/RP spectra afford a clean and efficient selection of metal-poor stars in the inner galaxy.  We do this by selecting RGB/RC stars towards the Galactic center (+-30 degrees) and deriving data-driven revised [M/H] estimates by

• deriving narrow-band filter fluxes (Stromgren and similar; a la "the Galaxy in your favourite colors")  -- Vedant Chandra
• deriving [M/H] estimates via xgboost in training with the APOGEE DR17 data set .

The initial query was:

SELECT  source_id,ra,dec,phot_g_mean_mag,phot_bp_mean_mag,bp_rp,parallax,parallax_error

FROM gaiadr3.gaia_source 

WHERE

parallax < 100.*power(10.,0.2*(0.9 - (phot_g_mean_mag - 1.5*(bp_rp-1.))))

and

parallax < 1.

and

abs(b)<30 and (l<30 or l > 330.)

and 

bp_rp between 1.0 and 3.5

and

phot_bp_mean_mag < 15.5

where the 1st condition selects for stars with M(de-reddened) < 0.9, i.e. eliminates the MS; the color cut is designed to eliminate stars bluer than the blue edge of the unreddened, metal-poor RC; and -- importantly -- phot_bp_mean_mag < 15.5 implies that the BP/Rp spectra should have good S/N in the blue.   This query returns 2.1M stars.

The features used in the [M/H] estimate by xgboost are <...> and the training set is the full APOGEE DR17 RC/RGB set.

The distribution of these stars in the B-R and Stromgren-m1 space is shown below, where the "streak" in the bottom-left corner is reddened hot stars, that we excise before training (they are not well represented in the training set).

Verification

The upshot is that it looks like we can select low-[M/H] samples cleanly and effectively. Here are a few plots that illustrate this:

Cross-/Self-Validation of [M/H] estimates (X-axis) against presumed ground-truth (APOGEE; Y-axis). The figure shows that a) for [M/H] > -0.9 one gets robust and precise [M/H] estimates (whose precision could presumably be improved). More pertinent for the questions at hand is that the selection below [M/H] = -0.9 works robustly, albeit with larger scatter. But stars stars selected to have [M/H] < -1 are at [M/H] < -1.

Resulting [M/H]-distribution of the overall 2.1M star sample of the "inner galaxy". This deserves more analysis, but the metal-poor end looks very plausible, and not swamped by false positives. In particular there are 10x fewer stars at [M/H] = -2 than at [M/H] = -1, as expected for naive closed-box models.

The next 3 panels show the on-sky distribution of the metal-rich sample ([M/H] > -0.4), an intermediate sample ( -0.9 < [M/H] < -0.4 ), and the metal-poor sample ( [M/H] < -0.9). The latter sample is (in projection) very much concentrated towards the Galactic center.

What gives credence to the fact that our metal-poor selection picks up stars in the innermost galaxy is the parallax distribution: the metal-poor stars (green) are peaked at 8.2 kpc (blue line0; the metal-rich stars (blue) are seen throughout the disk; intermediate [M/H]  stars are in orange.

HWR note to self; the flags were set wrong for the plots; the sample is too small.

Results: what do we see?

The first result is indeed the previous figure: there is a population of (~20k) stars that are metal-poor ([M/H] < 0.9) and very centrally concentrated: there is a metal-poor (and presumably ancient) heart of the MW.

That immediately brings up the question of "on what orbits do these low-[M/H] stars move"? Are they confined to the inner MW, or just passing through. 

In general about 3/4 of the sample have RVS velocities. The low-[M/H] sample is in the most crowded part and has only 1/3 RVS velocities. Fron those (Price-Whelan) we calculate orbits.

[Note distance uncertainties are tricky, and deserve attention]

This is what the apocenter -- eccentricity plane looks like:

The plot above excites HWR: the vast majority of stars have a broad eccentricity distribution from 0.2 to 0., and remain confined to the inner ~5 kpc. They are not just passing through. At high eccentricities there is a tail of stars with apocenters of 10kpc to >100 kpc, as expected e.g. for members of GSE, "just passing through" the center. [Need to check ho much of that is in the literature]

{Selection effects matter}

This plot overlays in orange the [M/H]-distribution of all stars on highly eccentric orbits (ecc>0.9). One can think of them as kinematically-selected GSE members.  If they were GSE members, the metal rich tail would seem puzzling. What's going on?

It gets potentially more interesting if we plot for the highly eccentric orbits, ecc>0.85, the distribution [M/H] vs r_{apo}:

 

One sees:

One interpretation is: at r_apo > 10kpc the [M/H]-distribution is quite GSE-like; for tightly bound orbits the [M/H] is higher: somet6hing different, or the more-metal-rich ex-core of GSE?

Implications and Next Steps

• tidy up analysis and write paper(s)
• think about getting ages

• make sure all the metal-poor stars here (and throughout the plane) are in SDSS-V targetting!!

on BP/RP metallicities in the disk

 Science goals:

Approach:

First implementation:

Results:

Exploration: where does the [M/H] information come from?

classic Stromgren filter combination

finding low M/H stars form Stromgren alone

RC vs not-clump stars

Are the simple ways to do an all-sky QSO selection with Gaia and WISE?

 What are good ways to select QSO samples before GDR3 ?

Goal

• Gaia will provide a QSO category with redshifts in DR3. There will be ~6M candidates; but, of course,  a significant fraction of them may not be QSOs. I.e. sample purity will be an  issue. And, redshift aliasing (from identification of the emission lines) may also be an issue.
• Are there ways to
1. reduce the non-QSO contamination?
2. help break redshift ambiguities?
• Here we focus on 1. ,  by asking how well (completeness and purity) can one select QSO samples with existing all-sky information.
• specifically we want to use Gaia and WISE only. 
• for the moment we'll stick to G<20, else Gaia data constraints become weak.
• we'll also stick to | b | >20 deg for now, to reduce near-plane contaminants.
• We want to exploit the three simple conditions
• quasars have 0 parallax  ("consistent with")
• quasars have 0 proper motion ("consistent with")
• most quasars z< 4 have W1-W2 > 0.5  (while most stars etc have W1-W2~0).

           Are the first two conditions sensible? The two plots below show the statistics of these quantities for the sample of 150k SDSS DR16 quasars G<20, selected independent of Gaia.

How far do we get with this?  Proposed approach: try it out on Stripe 82 and then apply to whole sky

Stripe 82 Experiment

• Basic numbers
• 8967 QSOs in S82, of which 5471 are G<20; of those 5242 also have a WISE match (<2").
• if we query:   

SELECT *

FROM gaiaedr3.gaia_source

WHERE  phot_g_mean_mag < 20.0

and parallax_over_error between -3 and 3

and  (pm*pm)/(pmra_error*pmra_error+pmdec_error*pmdec_error) < 9

and dec between -1.2 and 1.2 and (ra < 58 or ra > 309)

we get 20516 sources, of which 13040 have WISE matches.

         Their color distribution looks like this (plot below is for all sky with same Gaia Query)

         

This is the color-color distribution of non-moving, no-parallax Gaia sources G<20: top left lump: QSOs; sources with W1-W2~0 stars of all kinds of colors (presumably distant Giants, and -- se later LMC/SMC stars). For orientation, here the color distribution of spec. confirmed S82 QSO with G<20

This suggests to make a color cut as indicated below. Note that it turns out that S82 likely has overlooked quite a number of red/reddened QSO, which are in a plume towards the bottom-right.

What does that selection yield in S82?

Completeness: From the initial 20k sources Gaia selected, 13k have a WISE match, and 6500 are in the (blue) color-cut region above (compared to 5242 know S82 QSOs with WISE and G<20), after applying basic Gaia quality flags.  Of the 5242 QSOs with spectra, 4827 are picked up bi this selection: 92% completeness.

Purity: But there are 24% more sources in the Gaia-WISE-selected sample in S82 [TBD: need to check boundaries..] than spectroscopically conformed S82 QSOs. Conjecture: this is mostly (but not only), because SDSS is incomplete in red (reddened?) QSO, as the following plot suggests:

Here blue is the S82-spec sample, and gray the Gaia-WISE selection.

Last but not least, what if one applies the G<20, no-proper-motion, no-pls, red-in-(W1-W2) selection fo the entire sky with |b|>20?   

• the initial Gaia query yields 4.88M sources
• of these, 1.95M sources have a WISE match within 2"
• if one then makes the (B-R) -- (W1-W2) color cut (blue region a few plots up), one gets 747,000 QSO candidates with a sky distribution like this:

First addendum:

Are there smarter Gaia-WISE colors to single out QSOs?

Conjecture: yes, e.g the "dumbest" color, G-W1 vs W1-W2

if you then select the red area

You get a sky-distribution that is astoundingly uniform

What do we see:

• basic uniformity across sky. Yay!
• vestiges of the bulge, the magellanic clouds and one cluster with tidal tails in the North???
• imprint of dust-dimming
• a few articfacts that smell of Gaia-sky-scan

Second addendum:

What if we push to G=21?  Same procedure as above, yields another 600k objects with a rather uniform sky distribution..

Addendum  3:

What about photometric redshifts? Good enough to help with aliasing?

A qualitative look at the color-color plane looks promising:

And, using nearest neighbour, or NN (2nd plot) on z = f(G,BP,RP,W1,W2,A_G) yields

which looks "OK".

Questions:

• is that sample (700k after further cleaning) interesting, if we get phot-z + Gaia-spec-z for most?
• what happens if one pushed the same procedure fainter?