# Electric field-dependent dynamic polarizability and state-insensitive
conditions for optical trapping of diatomic polar molecules

Svetlana Kotochigova1∗ and David DeMille2 1Department of Physics, Temple
University, Philadelphia, PA 19122-6082, USA
2Department of Physics, Yale University, New Haven, CT 06520, USA

###### Abstract

Selection of state-insensitive or “magic” trapping conditions with ultracold
atoms or molecules, where pairs of internal states experience identical
trapping potentials, brings substantial benefits to precision measurements and
quantum computing schemes. Working at such conditions could ensure that
detrimental effects of inevitable inhomogeneities across an ultracold sample
are significantly reduced. However, this aspect of confinement remains
unexplored for ultracold polar molecules. Here, we present means to control
the AC Stark shift of rotational states of ultracold diatomic polar molecules,
when subjected to both trapping laser light and an external electric field. We
show that both the strength and relative orientation of the two fields
influence the trapping potential. In particular, we predict “magic electric
field strengths” and a “magic angle”, where the Stark shift is independent of
the DC external field for certain rotational states of the molecule.

The advantage of using state-independent light traps for precision frequency
measurements with ultracold atoms has been demonstrated in several experiments
Katori ; TrapScience08 . Applications of this approach were analyzed in the
context of optical atomic clocks and coherent control of atoms and photons
within an optical cavity. For these applications the frequency of the laser
beam that creates a far-off resonant optical dipole trap is chosen such that
the AC Stark shift of the ground and one excited electronic atomic level are
the same. In this way any optical atomic transition between these levels is
unaffected by the trapping light. In Refs. Flambaum ; Beloy ; Derevianko ;
Lundblad  state-insensitive trapping conditions have also been found for
microwave transitions in the atomic ground state by using a combination of the
vector and tensor components of the AC Stark shift and an external magnetic
field.

Ultracold molecular systems possess unique proporties that are considered to
make them potentially useful as tools for precision measurements special_issue
; NJP_review  and quantum computing DeMille02 ; Zhao . Therefore it is
desirable to extend the zero-differential AC Stark shift technique to these
more complex systems. The idea of “magic” frequencies for vibrational Raman
transitions in homo-nuclear Sr2 molecules was first explored in Refs.
Zelevinsky ; Kotochigova09 . This molecule is proposed for a search of
possible time variation of the electron-to-proton mass ratio. The AC Stark
shift of a molecule is determined by the dynamic molecular polarizability
$\alpha(\nu)$, which is a function of radiation frequency $\nu$ and its
polarization.

Polar molecules have a permanent dipole moment and their levels can be shifted
and mixed with one another by applying an external electric field. This opens
up a new way to create “magic” trapping conditions for two rotational levels
of the molecule. In the presence of an external electric field, $J$ is not a
good quantum number and all states, even the “rotationless” ground state, have
an anisotropic polarizability Ospelkaus . The anisotropy of the dynamic
polarizability of these levels manifests itself as a dependence on the
relative orientation of the polarization of the trapping laser and the DC
electric field. The combined action of these two fields can be a powerful tool
to manipulate and control ultracold molecules trapped in an optical potential.

The behavior described here has potential applications to several experiments
that have been envisioned for diatomic polar molecules held in optical traps.
The implications are particularly striking for the use of polar molecules in
an optical lattice as quantum bits. As was described in Ref. DeMille02 , a
pair of rotational states forms a suitable quantum bit. However, it was
pointed out that this system is susceptible to decoherence due to intensity
fluctuations in the optical trapping lasers, if the dynamic polarizabilities
of these levels differ as is generally the case. This in turn leads to very
stringent requirements on the laser intensity stability for such a system to
be practical. As we will show, this limitation can be removed by adjusting the
experimental parameters to guarantee that the dynamic polarizabilities of
these states are equal. In particular, for this proposed system (where a
spatial gradient of the electric field $\vec{\mathcal{E}}$ is required, so
working at a magic electric field value is impossible) it should be possible
to use light polarized at the “magic angle” relative to the static field in
order to eliminate this potentially dangerous source of decoherence.

The existence of magic electric field values is also of possible use for
envisioned applications where polar molecules in optical lattices are employed
in novel types of many-body systems. This includes, for example, cases where
the properties of long-range molecular interactions are tailored by a
combination of a static electric field and resonant microwave fields coupling
different rotational states Micheli06 ; Buchler06 ; Buchler07 . Working at a
magic electric field value in such systems could ensure that the inevitable
inhomogeneities in the intensity of the trapping light across a large sample
would not change the resonant condition for the microwave drive fields. Hence,
working under such “magic” conditions might be necessary to implement
proposals of this type.

Motivated by these ideas, we calculate the near-infrared dynamic
polarizability of various rotational levels of the $v=0$ vibrational state of
the X${}^{1}\Sigma^{+}$ potential of the KRb and RbCs molecules, under the
simultaneous influence of trapping electromagnetic and static electric fields.
We calculate dynamic polarizability using the computational techniques
developed in previous publications Kotochigova1 ; Kotochigova2 ; Kotochigova3
. Our calculations with an external electric field are predominantly performed
at a laser frequency of 9174 cm-1 (or wavelength of 1090 nm), which
corresponds to an often-used frequency to trap atoms and molecules in
ultracold experiments. In the near infrared this laser frequency is
sufficiently far away from molecular resonances of the excited electronic
states that heating due to photon scattering is negligible. We focus on
external electric field strengths up to 15 kV/cm, a value which should be
experimentally accessible.

Figure 1: Dynamic polarizability in the absence of an electric field and in
atomic units for the $v$=0, $J$=0 and 1 levels of the X${}^{1}\Sigma^{+}$
ground state of KRb and RbCs as a function of trapping frequency in the near
infrared domain. The dashed lines show the polarizability for the $J$=0 and
$M$ = 0 level, which is independent of the light polarization. The solid lines
correspond to the $J$= 1 state with magnetic sublevels $M$ = 0 and $\pm$ 1
illuminated by linear polarized light along the $\hat{z}$ direction. One
atomic unit of polarizability corresponds to 4.68645$\times 10^{-8}$
MHz/(W/cm2).

We begin by studying the dynamic polarizability or AC Stark shift for the
$J=0$ and $J=1$ rotational levels of the $v=0$ X${}^{1}\Sigma^{+}$ state of
KRb and RbCs in the near infrared, without an external electric field. The
procedure used to determine the complex molecular polarizability has already
been described in our previous paper Kotochigova1 . Figure 1 illustrates the
results of this study in the absence of an electric field by showing the
dynamic polarizability $\alpha$ as a function of laser frequency $h\nu$. The
range of laser frequencies spans a technologically relevant near-infrared
optical domain and includes the lasing frequency near 9174 cm-1 used in Ref.
Ospelkaus . The figure also shows that the dynamic polarizability of the
molecular state depends on the rotational quantum number $J$ and its
projection $M$. The curves for the different states do not cross as a function
of the laser frequency. In other words, at zero electric field we can not find
a “magic” frequency for J=0 and J=1 states using trapping light in the near
infrared.

Field-dependent dynamic polarizability

Here we extend the idea of the AC Stark shift for the rotational levels of the
X${}^{1}\Sigma^{+}$ ground state to the mixing of these levels in a static
electric field $\vec{\mathcal{E}}=\mathcal{E}\hat{z}$ along the space-fixed
$\hat{z}$ direction. Our analyses show that the experimentally accessible
electric fields mix only a few low-lying rotational states. The ultimate goal
is to introduce the dynamic polarizability of the mixed rotational levels as a
function of an external DC electric field and the polarization of the AC
trapping field. In a heteronuclear molecule this Stark mixing is primarily due
to the permanent electronic dipole moment $\vec{d}$; we neglect the much
smaller effects due to Stark mixing with other electronic states. For the
Hund’s case (a) X${}^{1}\Sigma^{+}$ state of an alkali-metal dimer, the total
electron spin and orbital angular momentum are not coupled to the molecular
rotation. The molecular wavefunction in the lab frame, $|vJM\rangle_{z}$, is
then given by

$|vJM\rangle_{z}\equiv|vJ\rangle\times
Y_{JM}(\hat{R})=\left\\{\frac{\psi_{vJ}(R)}{R}|{\rm
X}^{1}\Sigma^{+}\rangle\right\\}\times Y_{JM}(\hat{R})\,\,,$ (1)

where $v$ is a vibrational quantum number, $J$ and $M$ are the molecular
rotational angular momentum and its projection along the $z$ axis,
$\psi_{vJ}(R)$ is the radial rovibrational wavefunction, $Y_{JM}(\hat{R})$ is
a spherical harmonic, $\hat{R}$ is the orientation of the molecule relative to
the electric field direction $z$, $|{\rm X}^{1}\Sigma^{+}\rangle$ is the
electronic wavefunction with projections defined along the internuclear axis.
For the X${}^{1}\Sigma^{+}$ state, in even modest electric fields the nuclear
spins are well-decoupled from the other spins and angular momenta; hence we
ignore the nuclear spins here. Then the Hamiltonian for such a system becomes

$H=\sum_{vJM}E_{vJ}|vJM\rangle_{z}\times{}_{z}\langle
vJM|-\vec{d}\cdot\vec{\mathcal{E}}\,,$ (2)

where $E_{vJ}=G_{v}+B_{v}J(J+1)$ and $G_{v}$ and $B_{v}$ are the vibrational
energy and rotational constant of vibrational level $v$, respectively. Higher
order rotational corrections are negligible. For the alkali-metal dimers KRb
and RbCs we have that $\Delta G_{v}=G_{v+1}-G_{v}$ is on the order of 50-100
cm-1 for small $v$, while $B_{v}$ is on the order of 0.017-0.037 cm-1.

We evaluate the matrix elements of the operator
$-\vec{d}\cdot\vec{\mathcal{E}}$ by noting that after averaging over the
electronic wavefunction $|{\rm X}^{1}\Sigma^{+}\rangle$ it reduces to
$-d(R)C_{10}(\hat{R})\mathcal{E}_{z}$, where $d(R)=\langle{\rm
X}^{1}\Sigma^{+}|d|{\rm X}^{1}\Sigma^{+}\rangle$ is the $R$-dependent
permanent electric dipole moment,
$C_{lm}(\hat{R})=\sqrt{4\pi/(2l+1)}Y_{lm}(\hat{R})$ are tensors of rank $l$,
and $\mathcal{E}_{z}$ is the electric field strength. Consequently, this
operator conserves the projection quantum number $M$. The matrix element
between two rovibrational states is

${}_{z}\langle
vJM|-\vec{d}\cdot\vec{\mathcal{E}}|v^{\prime}J^{\prime}M^{\prime}\rangle_{z}=-\delta_{MM^{\prime}}\,d_{vJ,v^{\prime}J^{\prime}}\,\mathcal{E}_{z}\int
d\hat{R}\,Y_{JM}^{*}(\hat{R})C_{10}(\hat{R})Y_{J^{\prime}M^{\prime}}(\hat{R})$
(8)
$\displaystyle\quad\quad=-\delta_{MM^{\prime}}d_{vJ,v^{\prime}J^{\prime}}\mathcal{E}_{z}(-1)^{M}\sqrt{(2J+1)(2J^{\prime}+1)}\left(\begin{array}[]{ccc}J&1&J^{\prime}\\\
-M&0&M^{\prime}\end{array}\right)\left(\begin{array}[]{ccc}J&1&J^{\prime}\\\
0&0&0\end{array}\right)\,,$

where $\delta_{MM^{\prime}}$ is the Kronecker delta function,
$d_{vJ,v^{\prime}J^{\prime}}=\int_{0}^{\infty}dR\,\psi_{vJ}(R)d(R)\psi_{v^{\prime}J^{\prime}}(R)$,
and $(\cdots)$ are 3-$j$ symbols (see e.g. Brink ). This matrix element is
nonzero when $J+1+J^{\prime}$ is even, according to the parity selection
rules, and is independent of the sign of $M$ and $M^{\prime}$. For the small
$J$ values of interest here we can assume that the $J$ dependence of
$\psi_{vJ}(R)$ is negligible. Moreover, coupling between vibrational levels
can also be ignored, as the dipole moment $d(R)$ is a slowly varying function
with $R$ and the spacing between vibrational levels is large compared to the
rotational splitting. For tensor operators $C_{lm}(\hat{R})$ of rank 0, 1, and
2 the relationship between Cartesian $x$, $y$, and $z$ components and
spherical $m=-1$, $0$, and $+1$ components can be found in Refs. Bonin ;
Varsholovich .

For each projection $M$ and vibrational level $v$, eigenvalues and
eigenvectors of the Hamiltonian, Eq. (2), are obtained by the direct
diagonalization of the Hamiltonian matrix including rotational states $J$ from
$|M|$ up to some value $J_{\rm max}$. For the external electric field strength
and $J$ values of interest, $J_{\rm max}=10$ is sufficient for convergence. We
label the eigenenergies by $E_{v\tilde{J}M}$ with corresponding eigenvectors
$|v\tilde{J}M\rangle=\sum_{J}U^{vM}_{\tilde{J},J}|vJM\rangle_{z}$. Here
$\tilde{J}$ is an integer index with values $\tilde{J}=|M|,|M|+1,...$, such
that the eigenstate $|v\tilde{J}M\rangle$ adiabatically connects to the
electric field-free eigenstate $|vJM\rangle_{z}$ with $J=\tilde{J}$. For
$|M|>0$ the levels with projection quantum number $-M$ and $M$ remain
degenerate. The dipole matrix element between states $|v\tilde{J}M\rangle$ and
$|v^{\prime}\tilde{J}^{\prime}M^{\prime}\rangle$ of the X potential is given
by $\langle
v\tilde{J}M|d_{\sigma}|v^{\prime}\tilde{J}^{\prime}M^{\prime}\rangle=\sum_{J,J^{\prime}}U^{vM}_{\tilde{J},J}\,U^{v^{\prime}M^{\prime}}_{\tilde{J}^{\prime},J^{\prime}}\times{}_{z}\langle
vJM|d_{\sigma}|v^{\prime}J^{\prime}M^{\prime}\rangle_{z}$, where $\sigma$ is a
spatial index that can be expressed either in Cartesian coordinates $x,y,z$,
or spherical coodinates $q=0,\pm 1$.

Now we are able to calculate the dynamic polarizability of the mixed
rotational eigenstates at laser frequency $\nu$. It is determined by the
properties of the operator $\alpha_{\sigma\sigma^{\prime}}(\nu)$ defined by
Bonin ; Stone

$\displaystyle\alpha_{\sigma\sigma^{\prime}}(\nu)$ $\displaystyle=$
$\displaystyle\sum_{\gamma}\left\\{\frac{1}{E_{\gamma}-E_{v\tilde{J}M}+h\nu}+\frac{1}{E_{\gamma}-E_{v\tilde{J}M}-h\nu}\right\\}d_{\sigma}|\gamma\rangle\langle\gamma|d_{\sigma^{\prime}}$
(9)

with $\gamma$ enumerating eigenstates of the ground as well as excited
electronic potentials in the presence of the electric field and
$\sigma,\sigma^{\prime}=x,y,$ or $z$. We will focus on the $J=0$ and $J=1$
levels and without loss of generality assume that the $x$ axis of our
coordinate system lies in the plane spanned by the electric-field direction
and the orientation of the linearly-polarized laser light. Moreover, we are
interested in the situation where the level shifts due to the laser are small
compared to those induced by the electric field. For the isolated levels with
M=0, the polarizability is determined by the diagonal matrix element
$\alpha^{v\tilde{J}M=0}_{\sigma\sigma^{\prime}}(\nu)=\langle
v\tilde{J}M=0|\alpha_{\sigma\sigma^{\prime}}(\nu)|v\tilde{J}M=0\rangle$. The
$|M|>0$ dynamic polarizability needs to be treated by degenerate perturbation
theory within the two-dimensional subspace $|v\tilde{J}M\rangle$ and
$|v\tilde{J}-\\!M\rangle$. In fact, our choice of the $x$ direction and the
symmetry properties of the dynamic polarizability ensure that the linear
combinations
$|v\tilde{J}M,\pm\rangle=\\{|v\tilde{J}M\rangle\pm|v\tilde{J}-\\!M\rangle\\}/\sqrt{2}$
with $M>0$ are the correct eigenstates. For these states the dynamic
polarizability $\alpha_{\sigma\sigma^{\prime}}^{v\tilde{J}M,\pm}(\nu)=\langle
v\tilde{J}M,\pm|\alpha_{\sigma\sigma^{\prime}}(\nu)|v\tilde{J}M,\pm\rangle$.
For diatomic species, only the diagonal elements
$\alpha^{v\tilde{J}M,\pm}_{xx}$, $\alpha^{v\tilde{J}M,\pm}_{yy}$, and
$\alpha^{v\tilde{J}M,\pm}_{zz}$ are nonzero. Hence in an oscillating electric
field
$\vec{\mathcal{E}}_{o}(t)=\mathcal{E}_{o}(0)\textrm{Re}\\{\vec{\epsilon}e^{i2\pi\nu
t}\\}$ (where $\vec{\epsilon}$ is the complex unit vector indicating the
polarization), the state $|v\tilde{J}M,\pm\rangle$ shifts in energy by an
amount $\Delta E$ given by $\Delta
E=-\sum_{\sigma,\sigma^{\prime}}|\mathcal{E}_{o}(0)|^{2}\alpha^{v\tilde{J}M,\pm}_{\sigma\sigma^{\prime}}(\nu)\epsilon_{\sigma}\epsilon^{*}_{\sigma^{\prime}}/4$.
The sum is over the spatial indices $\sigma,\sigma^{\prime}=x,y,z$.

We are interested in near infrared laser frequencies, which are detuned away
from resonances with rovibrational levels of the electronically excited
potentials. In particular, we focus on wavelengths between 1000 nm and 1100
nm. Starting from the $v=0$ vibrational level of the X${}^{1}\Sigma^{+}$
states of KRb and RbCs, photons of this wavelength do not possess sufficient
energy to reach rovibrational levels of the electronically excited singlet
${}^{1}\Lambda$ potentials. The photons do have enough energy to reach
vibrational level of the triplet b${}^{3}\Pi$ potential considering only
single-photon excitations. However, such transitions require relativistic
spin-orbit coupling to the A${}^{1}\Sigma^{+}$ potential to acquire a non-zero
dipole matrix element; such spin-orbit induced couplings are small enough to
neglect under the conditions of interest here.

This allows us to make several approximations. Firstly, we can use non-
relativistic potentials and transition dipole moments and in the calculation
of the dynamic polarizability only consider singlet ${}^{1}\Sigma^{+}$ and
${}^{1}\Pi$ potentials. Secondly, assuming a large detuning (such that
$|E_{\gamma}-E_{v\tilde{J}M}-h\nu|\gg E_{vJM}$ for all states of interest) we
can neglect the electric field and rotational dependence of the energy
denominators in Eq. (9). Finally, we find that in the near infrared the
contribution to the polarizability from intermediate states $\gamma$ in the
ground X${}^{1}\Sigma^{+}$ state is small.

With this in mind the polarizability becomes

$\displaystyle\alpha^{v\tilde{J}M,\pm}_{\sigma\sigma^{\prime}}(\nu)$
$\displaystyle\cong$
$\displaystyle\sum_{J,J^{\prime}}U^{v|M|}_{\tilde{J},J}U^{v|M|}_{\tilde{J},J^{\prime}}\sum_{ev_{e}\Lambda}\left\\{\frac{1}{E_{ev_{e}\Lambda}-E_{v}+h\nu}+\frac{1}{E_{ev_{e}\Lambda}-E_{v}-h\nu}\right\\}$
$\displaystyle\quad\quad\quad\quad\times\sum_{J_{e}M_{e}}{}_{z}\langle
vJM,\pm|d_{\sigma}|ev_{e}J_{e}M_{e}\Lambda\rangle_{z}\,{}_{z}\langle
ev_{e}J_{e}M_{e}\Lambda|d_{\sigma^{\prime}}|vJ^{\prime}M,\pm\rangle_{z}\,,$

where the energy $E_{v}$ is a typical vibrational energy in the ground state
potential, and the rovibrational wavefunctions
$|ev_{e}J_{e}M_{e}\Lambda\rangle$ of the electronically excited states with
approximate energy $E_{ev_{e}\Lambda}$ are given by

$|ev_{e}J_{e}M_{e}\Lambda\rangle_{z}\equiv|ev_{e}\Lambda\rangle\times|J_{e}M_{e}\Lambda\rangle_{z}\equiv\left\\{\frac{\phi_{v_{e}}(R)}{R}|e~{}^{1}\Lambda\rangle\right\\}\,\times\left\\{\sqrt{\frac{2J+1}{4\pi}}D^{J_{e}*}_{M_{e}\Lambda}(\hat{R})\right\\}\,,$

where the vibrational and electronic dependence has been isolated in the ket
$|ev_{e}\Lambda\rangle$ and the rotational dependence in
$|J_{e}M_{e}\Lambda\rangle_{z}$, respectively. The wavefunction
$\phi_{v_{e}}(R)$ is the radial rovibrational wavefunction, the ket
$|e~{}^{1}\Lambda\rangle$ is the electronic state with projection quantum
number $\Lambda$ defined along the internuclear axis, and
$D^{J}_{MM^{\prime}}(\hat{R})$ is a Wigner rotation matrix that describes a
symmetric top rotational wavefunction. The sum over $\Lambda$ includes both
positive and negative values, where $\Lambda=0$ corresponds to excited
${}^{1}\Sigma^{+}$ electronic states and $\Lambda=\pm 1$ to ${}^{1}\Pi$
states.

As stated before the energy of the excited state and $\phi_{v_{e}}(R)$ depend
on the electronic state $|e~{}^{1}\Lambda\rangle$ and vibrational level
$v_{e}$, but not $J_{e}$ and $M_{e}$. Consequently, the transition dipole
moments separate into ${}_{z}\langle
vJM,\pm|d_{\sigma}|ev_{e}J_{e}M_{e}\Lambda\rangle_{z}=\langle
vJ|d(R)|ev_{e}\Lambda\rangle F^{J_{e}M_{e}\Lambda}_{JM\pm,\sigma}$, where
$F^{J_{e}M_{e}\Lambda}_{JM\pm,\sigma}$ is an integral of the product of three
Wigner rotation matrices $D^{J}_{MM^{\prime}}(\hat{R})$ over the orientation
of the molecule $\hat{R}$ that can be evaluated using angular momentum algebra
Brink .

Moreover, we have verified that $\langle vJ|d(R)|ev_{e}\Lambda\rangle$ is
nearly independent of $J$. The sums over $J_{e}$ and $M_{e}$ in Eq. (Electric
field-dependent dynamic polarizability and state-insensitive conditions for
optical trapping of diatomic polar molecules) can now be performed and we
finally find

$\displaystyle\langle
v\tilde{J}M,\pm|\alpha_{\sigma\sigma^{\prime}}(\nu)|v\tilde{J}M,\pm\rangle=$
$\displaystyle\sum_{\Lambda}\alpha^{\Lambda}(\nu)\sum_{J,J^{\prime}}U^{v|M|}_{\tilde{J},J}U^{v|M|}_{\tilde{J},J^{\prime}}\times\sum_{J_{e}M_{e}}F^{J_{e}M_{e}\Lambda}_{JM\pm,\sigma}F^{J_{e}M_{e}\Lambda}_{J^{\prime}M\pm,\sigma^{\prime}}$

with the $\Lambda$-dependent

$\displaystyle\alpha^{\Lambda}(\nu)=\sum_{ev_{e}}\langle
vJ|d|ev_{e}\Lambda\rangle\langle
ev_{e}\Lambda|d|vJ\rangle\left\\{\frac{1}{E_{ev_{e}\Lambda}-E_{v}+h\nu}+\frac{1}{E_{ev_{e}\Lambda}-E_{v}-h\nu}\right\\}\,.$
(12)

The parallel $\alpha^{0}(\nu)$ and perpendicular $\alpha^{1}(\nu)$
contributions to the polarizability are due to transitions to the $\Sigma$ and
$\Pi$ states, respectively.

Magic DC Electric Field

Figure 2: Dynamic polarizability at a wavenumber of 9174 cm-1 of the $v$=0,
$\tilde{J}$=0 and 1 levels of the X${}^{1}\Sigma^{+}$ ground state of KRb (top
row) and RbCs (bottom row) as a function of the field strength of an external
electric field. The polarization of the trapping field is parallel (left
panel) and perpendicular (right panel) to the direction of the electric field.
The circles indicate a crossing point between $\tilde{J}$=0 and 1
polarizabilities. The degeneracy of the states $\tilde{J}=1$ and $M=\pm 1$ is
lifted by perpendicularly polarized laser light. The new states are labeled
$|+1\rangle\pm|-1\rangle$, corresponding to $|\tilde{J}=1,M\pm\rangle$
respectively. The $M=\pm 1$ states remain degenerate for parallel
polarization.

The polarizability depends on the stength of an external electric field
through the use of the unitary matrices $U^{vM}_{\tilde{J},J}$ of Eq.
(Electric field-dependent dynamic polarizability and state-insensitive
conditions for optical trapping of diatomic polar molecules). Figure 2 shows
the dynamic polarizability of the ground states of the KRb and RbCs molecules
as a function of the external electric field strength. The left panels of
these figures correspond to the parallel $(\vec{\epsilon}=\hat{z})$ and the
right panels to the perpendicular $(\vec{\epsilon}=\hat{x})$ polarization of
the trapping light relative to the electic field direction
$\vec{\mathcal{E}}$. In all cases the polarizability depends on both the state
index $\tilde{J}$ and the angular momentum projection $M$. For KRb (Fig. 2,
top row) a “magic” electric field strength exists at $\mathcal{E}=10$ kV/cm,
where the polarizability of $\tilde{J}$=0, $M$=0 and $\tilde{J}$=1, $M$=0
states coincide. This is possible due to the polar character of the molecules.
For RbCs (Fig. 2, bottom row) two “magic” electric field strengths exist. For
the pair of states $\tilde{J}=0,M=0$ and $\tilde{J}=1,M=0$ a crossing occurs
near 2 kV/cm. Another crossing appears for the states $\tilde{J}=1,M=0$ and
$\tilde{J}=1,M=\pm 1$ at 4.7 kV/cm. They both are at a much smaller field
strength than for KRb, since RbCs has a larger permanent dipole moment and
smaller rotational splittings in the ground state. Note that the $M=0$ “magic”
electric field occurs at the same field strength for both parallel and
perpendicular polarization of the trapping light. In fact, they are the same
for any polarization.

Magic Angle

Figure 3: Dynamic polarizability of the $v$=0, $\tilde{J}$=0, $M$= 0 level of
the X${}^{1}\Sigma^{+}$ ground state of RbCs as a function of external
electric field strength and the angle $\theta$ between the direction of the
electric field $\vec{\mathcal{E}}$ and the polarization $\vec{\epsilon}$ of
the trapping light at a wavenumber of 9174 cm-1. The contour lines are marked
by the polarizability value in atomic units.

Figure 4: Dependence of the dynamic polarizability on the angle between the
polarization $\vec{\epsilon}$ of the trapping light at a wavenumber of 9174
cm-1 and the direction of the external electric field $\vec{\mathcal{E}}$, for
different values of $\mathcal{E}$. The polarizability of states with
projections $M=0$ and $M=\pm 1$ of the X${}^{1}\Sigma^{+}$ ground state of KRb
and RbCs are shown in the left and right panels, respectively. For clarity
$\mathcal{E}$ ranges from 0.6 kV/cm to 6 kV/cm in steps of 0.6 kV/cm for KRb
and from 0.3 kV/cm to 3 kV/cm in steps of 0.3 kV/cm for RbCs. At the “magic
angle” $\theta=\theta_{0}=54^{o}$, the polarizability of the states
$\tilde{J}=0,M=0$ and $\tilde{J}=1,M=0$ are identical and independent of the
strength of the static electric field. For the $\tilde{J}=1,M=\pm 1$ states no
magic angle exists.

In future experiments one would expect to be able to change the angle between
the static and dynamic electric fields. Figure 3 shows a surface plot of the
dynamic polarizability as a function of $\mathcal{E}$ and the angle $\theta$
between $\vec{\mathcal{E}}$ and polarization of the trapping light
$\vec{\epsilon}$, for linearly polarized light. The polarizability depends
smoothly on both $\theta$ and $\mathcal{E}$.

Figure 4 compares the dynamic polarizability of the states
$|\tilde{J}=0\rangle$ and $|\tilde{J}=1,M=0,\pm 1\rangle$ as a function of the
angle $\theta$ of the linear polarization of the optical field relative to the
static field (such that
$\hat{\epsilon}=\cos{\theta}\hat{z}+\sin{\theta}\hat{x}$), for several static
electric field strengths within the range from 0 to 6 kV/cm for KRb and 0 to 3
kV/cm for RbCs. All curves with $M=0$ cross at the angle $\theta=\theta_{0}$
such that cos${}^{2}\theta_{0}$ = 1/3, or $\theta_{0}\approx 54$ degrees. We
refer to this as the ”magic angle” since here the AC Stark shift is
independent of the internal state of the molecule. This behavior occurs in
many contexts and is a simple consequence of the rank-2 tensor structure of
the polarizability BKDbook . The angular dependence of the $\tilde{J}=0$ level
is smaller than that of the $\tilde{J}=1$ levels as expected. It inherited the
zero-electric field properties of the scalar $J=0$ state. As seen in Fig. 2,
the RbCs polarizability of the $\tilde{J}=0,M=0$ and $\tilde{J}=1,M=0$ levels
cross at a magic electric field strength of 2 kV/cm. For Fig. 4 this implies
that the curves start to “overlap” when $\cal E$ is equal or larger than this
magic field value.

The polarizability operator $\alpha_{\sigma\sigma^{\prime}}(\nu)$, defined by
Eq. (9), is a reducible rank-two tensor operator. Hence it can be expressed as
a sum of irreducible tensor operators $\alpha^{(k)}(\nu)$ of rank $k=0,1,2$.
In terms of these irreducible tensor operators, the AC Stark shift $\Delta E$
is proportional to the diagonal matrix element of the operator
$\sum_{\sigma\sigma^{\prime}}\alpha_{\sigma\sigma^{\prime}}(\nu)\epsilon_{\sigma}\epsilon^{*}_{\sigma^{\prime}}$,
which can be written in the general form

$\sum_{k=0}^{2}\sum_{q=-k}^{k}(-1)^{q}\alpha^{k}_{q}(\nu)\epsilon\epsilon^{k}_{-q},$
(13)

where $q$ is a spherical tensor projection index and the explicit forms of the
irreducible spherical tensors $\alpha^{k}_{q}(\nu)$ and
$\epsilon\epsilon^{k}_{q}$ are given in Refs Bonin ; Varsholovich .

In order to derive the “magic angle” condition for $M=0$ states, we consider
the effect of each term in the expansion (13). The term with $k=0$ corresponds
to the scalar polarizability; the operator
$\alpha^{0}_{0}(\nu)=\sum_{\sigma\sigma^{\prime}}\alpha_{\sigma\sigma^{\prime}}(\nu)\delta_{\sigma\sigma^{\prime}}$
has, under our approximations, diagonal matrix elements that are independent
of $J$ (or $\tilde{J}$) and $M$ for all states of interest, and similarly the
quantity
$\epsilon\epsilon^{0}_{0}\propto\vec{\epsilon}\cdot\vec{\epsilon}\,{}^{*}=1$
is independent of the polarization of the optical field. The term with $k=1$,
corresponding to the vector polarizability, in general is significant.
However, for the special case of linearly polarized light where
$\vec{\epsilon}$ is real, the quantities
$\epsilon\epsilon^{1}_{q}=(\vec{\epsilon}\times\vec{\epsilon}\,{}^{*})_{q}$
vanish and hence the effect of the vector polarizability is zero. For the
tensor polarizability terms (with $k=2$), from the Wigner-Eckhart theorem only
the operator component with $q=0$ gives rise to a non-zero diagonal matrix
element for $M=0$ states. Hence the contribution of this term is proportional
to $\epsilon\epsilon^{2}_{0}$. Without loss of generality we can define the
linear polarization as
$\vec{\epsilon}=\epsilon_{x}\hat{x}+\epsilon_{z}\hat{z}=\cos{\theta}\hat{z}+\sin{\theta}\hat{x}$.
In this case
$\epsilon\epsilon^{2}_{0}\propto\epsilon_{z}\epsilon_{z}^{*}-\vec{\epsilon}\cdot\vec{\epsilon}\,{}^{*}/3=\cos^{2}{\theta}-1/3$.
Hence the contribution to $\Delta E$ due to the tensor polarizability also
vanishes for all states with $M=0$, when the optical field is linearly
polarized at the “magic angle” $\theta=\theta_{0}$. Under this condition the
only contribution to the dynamic polarizability is from the $k=0$ scalar term,
which is the same for all states of interest.

The property of the “magic” electric field discussed in Fig. 2 can also be
understood in terms of the tensor structure of the polarizability. At certain
values of the applied DC electric field $\mathcal{E}$, the rank-2 components
of the polarizability, $\alpha^{2}_{0}$, of the $\tilde{J}=0,M=0$ and
$\tilde{J}=1,M=0$ levels becomes the same. When this condition is met, the AC
Stark shift becomes independent of the direction of the trapping light’s
linear polarization.

Here we finish with an example of a specific implementation of a “magic angle”
3-D lattice. Let the lattice be formed by three orthogonal retroreflected
laser beams $a,b$, and $c$, with initial propagation directions $\hat{k}$
given by $\hat{k}_{a}=\hat{y}$, $\hat{k}_{b}=(\hat{x}+\hat{z})/\sqrt{2}$, and
$\hat{k}_{c}=(\hat{x}-\hat{z})/\sqrt{2}$. These three beams each have a
different frequency $\nu$, such that
$\nu_{a}=\nu_{b}-\delta_{b}=\nu_{c}-\delta_{c}$; here the offset frequencies
$\delta_{b,c}$ must satisfy $\nu_{a,b,c}\gg\delta_{b,c}\gg f_{mot}$, where
$f_{mot}$ is the motional frequency of the molecules in the optical trapping
potential. The use of different frequencies (which can be generated from a
single laser by using e.g. acousto-optic modulators) in this manner eliminates
the effect of interference terms between the different laser beams: such terms
average to zero rapidly over the time of motion of the atom, and hence can be
neglected. The resulting average trap potential is then simply the sum of the
potentials due to each individual laser beam. Finally, the polarizations of
the three beams can be chosen as
$\hat{\epsilon}_{a}=\sqrt{2/3}\hat{x}+\sqrt{1/3}\hat{z}$;
$\hat{\epsilon}_{b}=\sqrt{2/3}\hat{k}_{c}+\sqrt{1/3}\hat{y}$; and
$\hat{\epsilon}_{c}=\sqrt{2/3}\hat{k}_{b}+\sqrt{1/3}\hat{y}$. In each case,
$|\hat{\epsilon}\cdot\hat{z}|=\cos{\theta_{0}}$.

## References

  * (1) J. Ye, H. J. Kimble, and H. Katori, Science 320, 1734 (2008).
  * (2) H. Katori, M. Takamoto, V. G. Pal’chikov, V. D. Ovsiannikov, Phys. Rev. Lett. 91, 173005 (2003).
  * (3) V. V. Flambaum, V. A. Dzuba, and A. Derevianko, Phys. Rev. Lett. 101, 220801 (2008).
  * (4) K. Beloy, A. Derevianko, V. A. Dzuba, and V. V. Flambaum, Phys. Rev. Lett. 102, 120801 (2009).
  * (5) N. Lundblad, M. Schlosser, and J. V. Porto, Phys. Rev. A 81, 031611(R) (2010).
  * (6) A. Derevianko, Phys. Rev. A 81, 051606(R) (2010).
  * (7) J. Doyle, B. Friedrich, R. V. Krems, and F. Masnow-Seeuws, Special issue on ulatracold polar molecules, Eur. Phys. J. D 31, 149 (2004).
  * (8) L.D. Carr, D. DeMille, R.V. Krems, Jun Ye, New J. Phys. 11, 055049 (2009).
  * (9) D. DeMille, Phys. Rev. Lett. 88, 067901 (2002).
  * (10) R. Zhao, el al. Nat. Phys. 5, 100 (2009).
  * (11) T. Zelevinsky, S. Kotochigova, and Jun Ye, Phys. Rev. Lett. 100, 043201 (2008).
  * (12) S. Kotochigova, T. Zelevinsky, and Jun Ye, Phys. Rev. A 79, 012504 (2009).
  * (13) S. Ospelkaus, K.-K. Ni, M. H. G. de Miranda, B. Neyenhuis, D. Wang, S. Kotochigova, P. S. Julienne, D. S. Jin, and J. Ye, Faraday Discuss. 142, 351 (2009).
  * (14) A. Micheli, G.K. Brennen, and P. Zoller, Nature Phys. 2, 341 (2006).
  * (15) H. P. Büchler, E. Demler, M. Lukin, A. Micheli, N. Prokof’ev, G. Pupillo, and P. Zoller, Phys. Rev. Lett. 98, 060404 (2007).
  * (16) H. P. Büchler, A. Micheli, and P. Zoller. Nature Phys. 3, 726 (2007).
  * (17) S. Kotochigova and E. Tiesinga, Phys. Rev. A 73, 041405(R) (2006).
  * (18) S. Kotochigova, E. Tiesinga, P. S. Julienne, New J. Phys. 11, 055043 (2009).
  * (19) S. Kotochigova, New J. Phys. 12, 073041 (2010).
  * (20) D. M. Brink and G. R. Satchler, Angular momentum, (Clarendon Press, London, 1993).
  * (21) K.D. Bonin and V.V. Kresin, Electric-dipole polarizabilities of atoms, molecules, and clusters, (World Scientific, Singapore, 1997).
  * (22) D.A. Varsholovich, A.N. Moskalev, and V.K. Khersonskii, Quantum theory of angular momentum: irreducible tensors, spherical harmonics, vectors coupling coefficients, 3nj symbols, (World Scientific, Singapore, 1988).
  * (23) A. J. Stone, The theory of intermolecular forces, (Clarendon Press, London, 1996).
  * (24) D. Budker, D.F. Kimball, and D.P. DeMille, Atomic Physics: an exploration through problems and solutions, 2nd ed., (Oxford Univ. Press, Oxford, 2008).

##  I Acknowledgments

This work is supported by a MURI grant of the Air Force Office of Scientific
Research; by NSF; and (for DD) by DOE. SK acknowledges helpful discussions
with J. Ye, D. Jin, and B. Neyenhuis.


11institutetext: Max-Planck-Institut für Astronomie, Königstuhl 17, 69117
Heidelberg, Germany
11email<EMAIL_ADDRESS>22institutetext: Department of Astronomy, University
of Michigan, Ann Arbor, MI 48109, USA

# Benchmarking Gaia DR3 Apsis with the Hyades and Pleiades open clusters

Wolfgang Brandner 11 Per Calissendorff 22 Taisiya Kopytova 11

(Received 2 May 2023 / Accepted 5 June 2023)

###### Abstract

Context. The Gaia astrophysical parameters inference system (Apsis) provides
astrophysical parameter estimates for several to 100s of millions of stars.

Aims. We aim to benchmark Gaia DR3 Apsis.

Methods. We have compiled $\approx$1500 bona fide single stars in the Hyades
and Pleiades open clusters for validation of PARSEC isochrones, and for
comparison with Apsis estimates. PARSEC stellar isochrones in the Gaia
photometric system enable us to assign average ages and metallicities to the
clusters, and mass, effective temperature, luminosity, and surface gravity to
the individual stars.

Results. Apsis does not recover the single-age, single-metallicity
characteristic of the cluster populations. Ages assigned to cluster members
seemingly follow the input template for Galactic populations, with earlier-
type stars systematically being assigned younger ages than later-type stars.
Cluster metallicities are underestimated by 0.1 to 0.2 dex. Effective
temperature estimates are in general reliable. Surface gravity estimates
reveal strong systematics for specific ranges of Gaia BP-RP colours.

Conclusions. We caution that Gaia DR3 Apsis estimates can be subject to
significant systematics. Some of the Apsis estimates, like metallicity, might
only be meaningful for statistical studies of the time-averaged Galactic
stellar population, but are not recommended to be used for individual stars.

###### Key Words.:

open clusters and associations: individual: Hyades, Pleiades – Stars:
abundances – Stars: evolution – Stars: fundamental parameters – Galaxy:
stellar content – Hertzsprung-Russell and C-M diagrams

##  1 Introduction

Precisely determined stellar parameters are the key to many astrophysical
investigations. The Gaia astrophysical parameters inference system (Apsis,
Bailer-Jones et al. 2013) was devised to extend the scientific content and
potential of Gaia beyond its core astrometric mission of providing precise
measurements of positions, parallaxes, and proper motions of more than 109
sources. Gaia Data Release 3 (DR3) includes astrophysical parameters of up to
470 million sources. The Apsis estimates are based on averaged low-resolution
spectra in the Gaia BP and RP bands with a resolution of R=20 to 60 combined
with Gaia photometry and astrometry, and on data from the Radial Velocity
Spectrograph (RVS), which provides R$\approx$11500 spectra in the wavelength
range 846 to 870 nm for stars with G$\lesssim$15.2 mag. From this stellar
astrophysical parameters like effective temperature, age, mass, surface
gravity log g, global metallicity [M/H] and elemental abundances, etc. are
inferred (Creevey et al. 2022; Fouesneau et al. 2022; Delchambre et al. 2022).

Gaia Apsis promises to unlock vastly improved insights into a variety of
research topics, like the Solar neighbourhood and young nearby moving groups,
stellar structure and evolution, open and globular clusters, and the overall
structure of the Milky Way with its distinct streams, bar, and spiral arms,
and more diffuse disk and halo components. Precise stellar astrophysical
parameters are also essential for the study of exoplanets, whose properties
are in general determined relative to the properties of the stellar host (e.g.
Magrini et al. 2022; Berger et al. 2023).

As pointed out by Bailer-Jones et al. (2013), the veracity of the Apsis
methods can only be assessed on real data. This provides our motivation to use
independently determined properties of bona fide single stars in the Hyades
and Pleiades open clusters as external validations of — and benchmarks for —
Gaia DR3 Apsis estimates.

The structure of the paper is as follow. In Sect. 2 we review the compilation
of the single star sequence of the Hyades open cluster, and apply the same
methodology to the Pleiades to validate PARSEC isochrones. In Sect. 3 we
compare cluster ensemble and individual stellar astrophysical properties with
Apsis estimates. In Sect. 4 we discuss the findings and potential explanations
for some of the issues identified. We conclude in Sect. 5 with recommendations
on the utilisation of DR3 Apsis estimates and an outlook.

##  2 Single star sequences in the Hyades and Pleiades

Figure 1: Colour-absolute magnitude diagram of bona fide single stars in the
Pleiades open cluster based on GAIA DR3. Overplotted is the family of the best
fitting PARSEC isochrones.  Figure 2: Minimum distance in colour-magnitude
space of each bona fide single star in the Pleiades from the best fitting
isochrone, plotted against the absolute magnitude (lower abscissa) and colour
(upper abscissa) on the isochrone. A positive residual corresponds to stars,
which are brighter (redder) than predicted by the isochrone. A negative
residual corresponds to stars, which are fainter (bluer) than the isochrone.
The uncertainty for each source in Gabs corresponds to the Gaia uncertainty in
absolute G magnitude, and the uncertainty in the residual corresponds to the
quadratically added uncertainties in absolute G magnitude and BP-RP colour.
Table 1: Compilation of abundance, $\alpha_{\rm ML}$, and age estimates for
the Pleiades open cluster based on isochrone fitting

[Fe/H] | X | Y | Z | $\alpha_{\rm ML}$ | mass range | age | PD1 | reference
---|---|---|---|---|---|---|---|---
|  |  |  |  | [M⊙] | [Myr] |  |
| 0.68 | 0.30 | 0.02 | 1.6 | [0.8,4.8] | 100 | BD | Meynet et al. (1993)
|  |  |  |  | [0.03,4.8] | 120 | DA |  Bouy et al. (2015)2
$+0.03\pm 0.05$ |  |  |  |  | [$\approx$2.5] | $95\pm 35$ | TY | Brandt & Huang (2015)
$>+0.29$ |  |  |  | 1.82 | [0.95,4.8] | $125\pm 35$ | TY | Gossage et al. (2018)
$+0.18$${}^{+0.29}_{-0.15}$ |  |  |  | 1.82 | [0.6,4.8] | 160${}^{+180}_{-60}$ | 2M | Gossage et al. (2018)
|  |  |  |  | [0.2,4.8] | 120$\pm 20$ | G3 | Heyl et al. (2022)

> 1key to photometric data set (PD): 2M - based on 2MASS photometry (Cutri et
> al. 2003); BD - based on BDA (Mermilliod 1992); DA - based on DANCe (Bouy et
> al. 2013); G3 - based on GAIA EDR3 photometry (Gaia Collaboration et al.
> 2021a); TY - based on TYCHO photometry (Høg et al. 2000); 2aimed at deriving
> empirical isochrones, thus only a comparison to theoretical isochrones, but
> no formal isochrone fitting.

Table 2: Astrophysical parameters of bona fide single stars in the Pleiades
open cluster

Gaia DR3 ID | RA | DEC | d | G | mass | $\sigma_{\rm mass}$ | log Teff | $\sigma_{\rm logTeff}$ | log L | $\sigma_{\rm logL}$ | log g | $\sigma_{\rm logg}$
---|---|---|---|---|---|---|---|---|---|---|---|---
| (deg) | (deg) | (pc) | (mag) | (M⊙) | (M⊙) | (K) | (K) | (L⊙) | (L⊙) | (cm/s2) | (cm/s2)
68828225308044288 | 52.797083 | 24.232163 | 129.717 | 17.3285 | 0.2011 | 0.0199 | 3.4712 | 0.0063 | -2.2212 | 0.0546 | 4.8196 | 0.0305
69335619861034752 | 52.816578 | 25.255258 | 130.680 | 8.0884 | 1.6837 | 0.0389 | 3.9026 | 0.0076 | 0.9233 | 0.0421 | 4.3043 | 0.0040
69585140280997760 | 52.824851 | 26.028613 | 138.364 | 12.3085 | 0.7792 | 0.0082 | 3.6793 | 0.0024 | -0.6595 | 0.0148 | 4.6555 | 0.0033
69583177479805824 | 52.836332 | 25.959207 | 144.031 | 19.2850 | 0.1695 | 0.0039 | 3.4507 | 0.0017 | -2.3999 | 0.0274 | 4.7990 | 0.0254
68051390279853824 | 53.002084 | 23.774632 | 137.339 | 13.3061 | 0.6515 | 0.0073 | 3.6234 | 0.0032 | -0.9824 | 0.0175 | 4.6816 | 0.0045

111Table 2 is published in its entirety in the machine-readable format. The
first six entries, and omitting some of the columns with uncertainties on the
distance estimates, and part of the Gaia DR3 photometry, are shown here for
guidance regarding its form and content.

At average distances of $\approx$45 pc (van Leeuwen 2009; Röser et al. 2011;
Gaia Collaboration et al. 2018) and $\approx$135 pc (e.g., Percival et al.
2005; Heyl et al. 2022), respectively, the Hyades and the Pleiades are among
the most nearby open clusters. The vicinity ensures high quality Gaia
astrometric, photometric, and spectroscopic measurements over a wide range of
spectral types and stellar masses. Age estimates for the Hyades derived from
isochrone fitting are in the range 500 to 800 Myr (Brandner et al. 2023b). For
the Pleiades the age estimates derived from isochrone fitting centre around
100 to 125 Myr (Table 1), which is in good agreement with age estimates based
on the Lithium depletion boundary (e.g., Basri et al. 1996; Stauffer et al.
1998; Dahm 2015; Galindo-Guil et al. 2022).

The stellar samples of candidate members of the open clusters are based on
Gaia Collaboration et al. (2021b) for the Hyades, and on Heyl et al. (2022)
for the Pleiades. In Brandner et al. (2023b) we processed the Hyades sample by
identifying and removing stars with dubious photometry in at least one of the
Gaia DR3 G, BP, or RP bands. Next we used the Gaia Renormalized Unit Weight
(RUWE) parameter to identify likely astrometric binaries. In a final step,
likely photometric (i.e. blended) binaries were identified based on their
distance from the median stellar sequence in colour-magnitude space, taking
into account the uncertainties in absolute brightness and colour. This
resulted in a sample of $\approx$600 bona fide single stars in the Hyades open
cluster. In Brandner et al. (2023a) we derived the astrophysical parameters
for these stars by comparison with the best fitting family of isochrones from
the PAdova and TRieste Stellar Evolution Code (PARSEC, Bressan et al. 2012;
Chen et al. 2014).

Here we apply the same methodology to the Pleiades. The sample of
$\approx$1300 stars within 10 pc of the nominal centre of the Pleiades open
cluster as defined by Heyl et al. (2022) includes 1191 sources with
$-0.2\leq$BP-RP$<$3.9 mag and $-2.0\leq$G${}_{\rm abs}\leq 13.5$ mag. After
rejection of 112 stars with dubious Gaia photometry, and the identification of
171 stars as likely astrometric or photometric binaries, we retain a sample of
908 bona fide single stars. We employ the PARSEC CMD 3.7 web
interface222http://stev.oapd.inaf.it/cmd to obtain version 1.2S grids of non-
rotating isochrones in the Gaia EDR3 photometric system for a range of ages
and [M/H], and for a visual extinction of A${}_{\rm V}=0.12$ mag (Meynet et
al. 1993; van Leeuwen 2009, and assuming A${}_{\rm V}=3.1*$E(B-V))333We
attribute the higher A${}_{\rm V}=0.20$ mag reported by Heyl et al. (2022) due
to their use of strictly solar metallicity [M/H] = 0 isochrones.. Absolute
Gabs magnitudes are based on the photo-geometric distance estimates provided
by Bailer-Jones et al. (2021). $\chi^{2}$ minimisation is used to find the
best fitting isochrones to the 235 brightest (G${}_{\rm abs}\leq 8.0$ mag)
bona fide single upper-main-sequence stars. This yields lg(age [yr]) =
$8.10\pm 0.07$ and [M/H] = $0.07\pm 0.03$ for the Pleiades open cluster.

The colour-magnitude diagram (CMD) of $\approx$900 bona fide single stars in
the Pleiades open cluster, and the family of best fitting PARSEC isochrones is
shown in Fig. 1. Similar to the Hyades, PARSEC isochrones provide an excellent
fit to the observed sequence for G${}_{\rm abs}\lessapprox 9.0$ mag. For
9.0$<$G${}_{\rm abs}<$11.0 mag, the isochrones predict fainter G magnitudes,
respectively bluer BP-RP colours, by up to 0.15 mag. For still fainter, mid-
to late M-dwarfs, the isochrones tend to predict too bright G magnitudes,
respectively too red BP-RP colours (Fig. 2).

For each star, we calculate its distances to the best fitting interpolated
isochrone in colour-magnitude space to assign Teff, $\log$ g, $\log$ L, and
(present-day) mass. The uncertainties in astrophysical parameters take the
uncertainties in the Gaia measurements, and in the [M/H] and age estimates of
the isochrone into account. The derived properties of the bona fide single
stars in the Pleiades are listed in Table 2. The corresponding stellar
properties for the Hyades open cluster are presented in Table 1 in Brandner et
al. (2023a).

##  3 Comparison with Gaia Apsis

###  3.1 Cluster ensemble properties: Age and Metallicity

Figure 3: Gaia DR3 Apsis age estimates as a function of the BP-RP colour for
single stars in the Hyades (left) and Pleiades (right). The blue horizontal
dashed-dotted lines mark the range of ages assigned to the open clusters by
isochrone fitting. For both clusters Apsis age estimates show a strong
correlation with the colour of the star for BP-RP $\leq 0.8$ mag.

Figure 4: [M/H] of single stars in the Hyades (left) and Pleiades (right),
which have both DR3 Apsis GSP-Phot and re-calibrated GSP-Spec estimates.
Marked as a black dashed line is the diagonal along which both [M/H] estimates
agree. The blue open squares mark the range of [M/H] estimates derived from
the best fitting PARSEC isochrones. Not shown are stars with Apsis estimates
of [M/H]$<-0.7$.

A defining characteristics of stellar populations in open clusters is their
small spread in age and abundances (Parmentier et al. 2014; Bovy 2016;
Kopytova et al. 2016; Magrini et al. 2017; Jeffries 2017; Krumholz & McKee
2020). In order to test how well Apsis recovers these essential properties, we
compare the estimates for individual stars with the ensemble average for each
cluster as derived by isochrone fitting. Apsis age estimates originate in the
Final Luminosity Age and Mass Estimator (FLAME) module (Bailer-Jones et al.
2013; Creevey et al. 2022). FLAME in turn relies on input both from the
Generalized Stellar Parametrizer – Photometry (GSP-Phot) and the Generalized
Stellar Parametrizer – Spectroscopy (GSP-Spec).

In Fig. 3 we plot FLAME mean age estimates along with the lower 16th and 84th
quantile for 109 stars in the Hyades and 186 stars in the Pleiades against
their BP-RP colour. FLAME lower and upper age estimates are capped between 0.2
and 13.5 Gyr444We note that this contradicts the lower limit of 0.1 Gyr quoted
in the Gaia DR3 documentation release version 1.2 (7. Feb. 2023), section
11.3.6. For both clusters, the FLAME age estimates reveal a strong correlation
with colour for BP-RP $\leq 0.8$ mag. For BP-RP $>1.0$ mag the majority of age
estimates tend to scatter around ages of $\approx$11 Gyr. Weighted-mean
uncertainties yield $<{\rm age_{FLAME}}>_{\rm wm}=1.73\pm 0.04$ Gyr for the
Hyades and $<{\rm age_{FLAME}}>_{\rm wm}=0.26\pm 0.01$ Gyr for the Pleiades.

The Apsis GSP-Phot and GSP-Spec modules independently provide [M/H] estimates.
In order to test for their consistency, we identified all stars, which have
[M/H] estimates from both modules. In Fig. 4 we plot the [M/H] estimates along
with diagonal dashed lines marking a 1:1 correspondence. GSP-Phot [M/H] are
the original published Apsis values. A recalibration using Code Version 1.6 of
the gdr3apcal Python package555gdr3apcal codebase: https://github.com/mpi-
astronomy/gdr3apcal according to Andrae et al. (2022) results in revised [M/H]
estimates for Pleiades stars with a wider intrinsic spread, and which are on
average 1.0 dex lower than the original Apsis estimates. As a low average
[M/H]$\approx-1.1$ appears non-physical for member stars of the Pleiades
(Funayama et al. 2009; Soderblom et al. 2009), we reject the re-calibration of
the GSP-Phot [M/H]. GSP-Spec [M/H] estimates are re-calibrated according to
Recio-Blanco et al. (2022), and using the polynomial coefficients applicable
to open clusters. This results in an average correction of $\approx$0.09 dex
towards higher [M/H] values.

###  3.2 Individual stellar properties: Teff, log g, mass

Figure 5: Effective temperature estimates for single stars in the Hyades (top)
and Pleiades (bottom) based on GSP-Phot (left) and GSP-Spec (right) as a
function of BP-RP colour. Overplotted in blue is the best fitting PARSEC
isochrone.

Figure 6: Surface gravity estimates for single stars in the Hyades (top) and
Pleiades (bottom) based on GSP-Phot (left) and re-calibrated GSP-Spec (right).
Not shown are stars with Apsis log g [cm/s2] estimates below 3.8. Overplotted
in blue is the best fitting PARSEC isochrone.

In Fig. 5, top, we plot the Apsis deduced effective temperatures for stars in
the Hyades as a function of their BP-RP colour. On the left we show the GSP-
Phot derived Teff, which reveals a monotonic and relatively smooth decline
with increasing colour up to BP-RP$\lessapprox$1.35 mag. In the range
1.35$<$BP-RP$<$3.0 mag, the GSP-Phot temperature estimates bifurcate.
Overplotted in blue is the Teff vs. colour relation according to the best
fitting PARSEC isochrone. Compared to the isochrone, GSP-Phot Teff estimates
are systematically lower by 150 to 250 K for BP-RP$\lessapprox$1.35 mag, while
the bifurcated GSP-Phot Teff estimates for BP-RP$>$1.35 mag seemingly bracket
the temperature according to the isochrone. For BP-RP$>$3.0 mag the GSP-Phot
Teff estimates are systematically higher by $\approx$50 K than implied by the
isochrone. On the right, we show GSP-Spec Teff, which also reveals a
bifurcation for the colour range 1.35$<$BP-RP$<$1.65 mag. For reference, the
corresponding stellar masses according to the isochrone are indicated on the
top abscissa.

Shown in the bottom of Fig. 5 are effective temperature vs. colour for the
Pleiades. Noticeable is the apparent gap in GSP-Phot derived Teff in the range
8500 to 9500 K with the majority of the stars with BP-RP$<$0.3 mag being
assigned T${}_{\rm eff}>$9500 K. For redder colours 0.3$<$BP-RP$<$1.3 mag GSP-
Phot Teff falls by a few 100 K systematically below – and for BP-RP$>$1.7 mag
above – the colour-temperature relation of the best fitting PARSEC isochrone.
We note that a few stars in the BP-RP colour range from 1.0 to 2.0 mag have
GSP-Phot temperature which fall significantly above the isochrone relation.
The GSP-Spec module could deduce effective temperatures only for 115 of the
bona fide single stars in the Pleiades. While the majority follows closely the
Teff–colour relation of the best fitting isochrone, a few stars stand out with
having significantly lower GSP-Spec temperatures.

In Fig. 6, top, we plot the Apsis deduced surface gravity for stars in the
Hyades as a function of their BP-RP colour. On the left we show the GSP-Phot
derived log g. Similar to the Teff estimates, log g follows a smooth,
monotonic relation with colour for BP-RP$\lessapprox$1.35 mag. In the colour
range 1.35$<$BP-RP$<$3.0 mag, GSP-Phot log g estimates bifurcate. For BP-
RP$<$3.0 mag GSP-Phot log g [cm/s2] estimates are typically 0.1 to 0.2 lower
than suggested by the isochrone, while they fall above the isochrone value by
$\approx$0.03 for redder colours. On the right we show the re-calibrated GSP-
Spec log g estimates according to Recio-Blanco et al. (2022). On average the
recalibration lowers log g [cm/s2] by $\approx$0.3 compared to the original
DR3 Apsis estimates (e.g. from 5.0 to 4.7). The log g [cm/s2] values stay
almost constant at 4.7 for 0.7$<$BP-RP$<$2.2 mag, with the exception of a dip
in the range 1.35$<$BP-RP$<$1.70 mag.

Shown in the bottom of Fig. 6 are the corresponding plots for the Pleiades.
GSP-Phot log g [cm/s2] estimates for BP-RP$<$1.7 mag fall below the isochrone
by up to 0.2. In the colour range 1.9$<$BP-RP$<$3.0 mag GSP-Phot log g
estimates bifurcate. For red colours (BP-RP$>$3.0 mag) the GSP-Phot log g
estimates are systematically higher than suggested by the best fitting
isochrone. Akin to the single stars in the Hyades, the re-calibrated GSP-Spec
log g [cm/s2] estimates for Pleiades members tend to cluster around 4.7.

##  4 Discussion of benchmarking results

###  4.1 Age and FLAME

The Apsis FLAME module is tuned using supervised learning (Rybizki et al.
2020) based on templates, which include a model of the Milky Way. Stars get
assigned a population identification (popID) derived from the age bins defined
by Robin et al. (2003) for the different Galactic populations, i.e. thin disk
(popID 0 to 6), thick disk (popID 7), halo (popID 8), plus additional popIDs
for bulge stars, stars belonging to the Magellanic Clouds, or to open
clusters. The apparent age-colour relation in Fig. 3 suggests that the Apsis
age inference is dominated by the model of the underlying Galactic population.
FLAME seems to assign stars a popID, and hence an ”age” primarily based on
their effective temperature.

FLAME uses a solar [M/H] prior (Creevey et al. 2022), which should not be a
major limitation towards its applicability to the near-solar [M/H] of the
Hyades ([M/H] = 0.18$\pm$0.03) and Pleiades ([M/H] = 0.07$\pm$0.03) open
clusters.

Of note is a potential mismatch between different stellar evolutionary models
employed by GSP-Phot, which provides inputs to FLAME, and FLAME itself.
According to Creevey et al. (2022) GSP-Phot employs PARSEC 1.2S Colibri S37
models (Tang et al. 2014; Chen et al. 2015; Pastorelli et al. 2020), i.e the
same set of models we use for the determination of the stellar parameters of
the single stars in the Hyades (Brandner et al. 2023a) and Pleiades. FLAME
itself, though, employs Bag of Stellar Tracks and Isochrones (BaSTI, Hidalgo
et al. 2018). Figure 12 in Hidalgo et al. (2018) highlights the close match of
BaSTI and MIST isochrones (Choi et al. 2016), and the mismatch between BaSTI
and PARSEC isochrones for young late-type stars of approximately solar
abundance.

As discussed in Brandner et al. (2023b), MIST models tend to systematically
predict fainter absolute magnitudes Gabs, respectively predict bluer BP-RP
colours for bona fide single stars in the Hyades with T${}_{\rm eff}\leq 4700$
K than is actually observed by Gaia DR3. PARSEC isochrones provide a
considerably better fit to the observed single star sequence (Brandner et al.
2023a).

A further limitation of Gaia DR3 FLAME is the apparent lower age limit of 200
Myr. This precludes a proper age dating of Pleiades members. In summary,
individual FLAME age estimates seem of very limited informative power.

###  4.2 Metallicity and GSP-Phot/GSP-Spec

In addition to the low resolution spectroscopic information in the BP and RP
bands, GSP-Phot [M/H] estimates also incorporate Gaia parallax measurements,
and use both theoretical model spectra and isochrones (Andrae et al. 2022).
Fouesneau et al. (2022) caution against the use of GSP-Phot [M/H] estimates
“without further investigation”, and state that [M/H] estimates appear on
average too low by 0.2 dex when compared to literature data for sources with
[M/H]$>-1$ dex. This is in agreement with our findings. Both for the Hyades
and Pleiades samples we find that the weighted mean [M/H] estimates of the
bona fide single stars with GSP-Phot [M/H] estimates fall below the estimates
derived from families of best fitting isochrones by 0.15 to 0.2 dex (Fig. 4).

GSP-Phot [M/H] estimates of individual stars, though, exhibit a much larger
scatter than supported by the 16% and 84% quantiles of the Markov-Chain Monte-
Carlo results reported in Apsis, suggesting that the uncertainties in the
parameter estimates for individual stars are underestimated by a factor of 10
to 20. A re-calculation of [M/H] using the Python package gdr3apcal (Code
Version 1.6) makes the discrepancy worse. Thus we agree with the
recommendation by Fouesneau et al. (2022) to not use GSP-Phot [M/H] estimates
of individual stars.

GSP-Spec [M/H] estimates are based on Gaia RVS and grids of model spectra
spanning the full range of Galactic stellar populations. The inherent
estimation bias of this module has been quantified by Recio-Blanco et al.
(2022), who also provide a prescription for correction. The correction moves
the ensemble [M/H] estimates of Apsis to within $\approx$0.1 dex to the [M/H]
values suggested from the isochrone fitting to the Hyades and Pleiades. For
individual stars, in particular in the Hyades, GSP-Spec [M/H] (minimal)
quantile uncertainties appear to be underestimated, indication that the actual
noise floor of the high S/N estimates is considerably higher.

While overall the GSP-Phot and GSP-Spec estimates show a weak correlation,
individual estimates tend to disagree by substantially more than suggested by
the assigned uncertainties. In general, neither GSP-Phot nor GSP-Spec are able
to reproduce the [M/H] estimates derived from the isochrone fitting, which are
marked by blue open squares. The apparent close match of the weighted mean
$<{\rm[M/H]_{GSP-Phot}}>_{\rm wm}=0.1720\pm 0.001$ for the Hyades is
accidental, as it is the result of a selection of stars with both GSP-Phot and
GSP-Spec estimates. A sample of 213 stars in the Hyades with GSP-Phot [M/H]
estimates, and including 96 stars without GSP-Spec estimates, yields
$<{\rm[M/H]_{GSP-Phot}}>_{\rm wm}=0.034\pm 0.001$. Overall, both GSP-Phot and
GSP-Spec tend to systematically underestimate [M/H] of the stars in the Hyades
and Pleiades.

###  4.3 Effective temperature and GSP-Phot/GSP-Spec

As the spectral characteristics of stars are primarily determined by their
effective temperature, it is no surprise that this is also the most robust of
the DR3 Apsis parameter estimates. Still GSP-Phot Teff estimates seem to
systematically fall by 150 K to 250 K below the Teff \- colour relation
suggested by the best fitting isochrones for BP-RP$<$1.35 mag. Surprising are
an apparent bifurcation of Teff estimates for stars in the Hyades with 1.35
mag$<$BP-RP$<$3.0 mag, and the apparent dearth of stars with 8000
K$<$T${}_{\rm eff}<$9500 K in the Pleiades. Apart from the systematic offsets,
the GSP-Phot Teff quantile uncertainties appear plausible (Fig. 5, left).

Overall GSP-Spec estimates are in better agreement with the Teff-colour
relation according to the PARSEC isochrone than GSP-Phot estimates, though
estimates start to diverge for BP-RP$>$1.35 mag (Fig. 5, right).

###  4.4 Surface gravity and GSP-Phot/GSP-Spec

GSP-Phot log g estimates overall follow the log g - colour relation suggested
by the best fitting isochrones. Of notice is a bifurcation in the log g
estimates in the colour range 2.0 mag$<$BP-RP$<$3.0 mag, and a systematic
underestimate of log g for BP-RP$<$1.5 mag (Fig. 6, left).

GSP-Spec log g [cm/s2] estimates, which have been re-calibrated according to
the prescription provided by Recio-Blanco et al. (2022) are almost constant at
4.7 over a wide colour range, with a strange dip to $\approx$4.4 at BP-
RP$\approx$1.4 mag for stars in the Hyades. We suspect that the lack of a
dimension for rotational line broadening in the GSP-Spec model grid might
force the parameter estimates in particular for fast rotating stars in the
Pleiades to higher apparent log g values (Fig. 6, right).

###  4.5 Limitations of our method

While the isochrone fits, and the deduction of stellar astrophysical
parameters from the isochrones provide robust estimates of stellar properties,
there are also some inherent limitations to our approach:

i) the $\chi^{2}$ isochrone fitting is based on the assumption of uniform
extinction to all stars in either of the clusters.666We note that due to their
vicinity, the line-of-sight extinction is very small towards the Hyades (AV=3
mmag), and modest for the Pleiades (AV=120 mmag). Apsis GSP-Phot, on the other
hand, fits individual line-of-sight extinction estimates for each star. While
this potentially provides more accurate fits, it also opens an additional
dimension for fitting errors due to a strong correlation between the
uncertainties in Teff and AV estimates (Andrae et al. 2022).

ii) our bona fide single star sequences are still contaminated by sources with
hidden binarity. These are typically close binary systems with mass ratios
and/or brightness ratios $<$1:2. We expect the contamination to be more severe
for the Pleiades sample due to its $\approx$3 times larger distance, which
reduces the sensitivity to binary detection and rejection based on the RUWE
value. Still, the bias induced by the fainter and redder secondary on the
parameter estimation in such system should be modest, and cannot explain the
systematics encountered in the Apsis estimates.

iii) we employ non-rotating isochrones. Stellar rotation affects in particular
the upper main sequence and main sequence turn-off region in the Hyades and
Pleiades, which are the regions also most sensitive to the isochronal ages.
For the majority of the stars in our sample, the effect on the determined
astrophysical parameters should be minimal, though.

iv) the PARSEC models follow a fixed Y(Z) abundance relation. This, however,
affects both our estimates and Apsis. We note that a deviation of the [Y/Z]
abundance ratio from the relation would in particular change the energy
production rate of the low mass stars (Baraffe & Chabrier 2018)

##  5 Recommendations and outlook

Table 3: Gaia DR3 Apsis quality assessment for main-sequence stars with 11500
K$\geq$T${}_{\rm eff}\geq$3000 K and [M/H]$\approx$0\.

Quantity | typical offset | recommendation2
---|---|---
FLAME age | varies with spectral type, up to $>$10 Gyr | do not use
GSP-Phot [M/H] | ensemble [M/H] 0.15 to 0.20 dex too low | do not use for individual stars
GSP-Phot Teff |  $\approx$200 K too low for BP-RP$<$1.35 mag, bifurcates for BP-RP$>$1.35 mag | ok with $\pm$250 K uncertainty
GSP-Phot log g | 0.1 to 0.2 cm/s2 too low for BP-RP$<$3.0 mag | ok with $\pm$0.2 cm/s2 uncertainty
GSP-Spec [M/H]1 | ensemble [M/H] $\approx$0.1 dex too low | ok with $\pm$0.2 dex uncertainty
GSP-Spec Teff | bifurcates for 1.35$<$BP-RP$<$1.65 mag | ok
GSP-Spec log g1 | non-physical variation with spectral type | do not use

7771Recalibrated according to Recio-Blanco et al. (2022). 2The uncertainties
quoted in the column recommendation should be understood as systematic
uncertainties on top of the quantile estimates provided by Apsis. “ok”
indicates that the estimates are largely in agreement with estimates according
to PARSEC stellar models.

While Gaia offers unprecedented quality and precision in its astrometric,
photometric, and spectroscopic measurements, we identify tension between Gaia
DR3 Apsis estimates and astrophysical properties deduced by us for
$\approx$1500 bona fide single stars in the Hyades and Pleiades open clusters.

Apsis provides very good to fair estimates on the effective temperature of
stars, and potentially on the ensemble [M/H] of homogeneous stellar groups.
Both GSP-Phot and GSP-Spec, though, seem to underestimate the ensemble [M/H]
of two nearby open clusters by 0.1 to 0.2 dex. As recommended by Andrae et al.
(2022), one should not use the GSP-Phot [M/H] for individual stars. GSP-Phot
log g in general is sufficiently precise to deduce luminosity classes, while
GSP-spec log g [cm/s2] for stars in the open clusters converge around 4.7 over
a larger range of spectral types. Our assessment of Gaia DR3 Apsis and the
resulting recommendations are summarized in Table 3.

We also have identified a few potential issues in the underlying assumptions
and training sets used for Gaia DR3 Apsis:

i) GSP-Phot employs PARSEC isochrones. FLAME relies on BaSTI models, whose
isochrones closely resemble MESA isochrones, but disagree with PARSEC
isochrones for solar metallicity stars with masses $\lesssim$0.7 M⊙. This
model mismatch might bias FLAME estimates

ii) The GSP-Spec spectral model grid lacks a dimension for rotational
broadening, hence fast rotating stars might get assigned too high log g
estimates

iii) FLAME age estimates reproduce the training model. The ages start at 100
Myr according to the online documentation, while the youngest (16th quantile)
ages assigned to Pleiades members are 200 Myr. FLAME lacks the ability to
classify the youngest (age $\lesssim$150 Myr) Galactic population (popID = 0
according to Rybizki et al. (2020)).

Addressing these issues might help to improve the quality of Apsis estimates.
Overall, DR3 Apsis seems to assign too small uncertainties to its estimates.
Thus Gaia DR3 Apsis cannot replace more dedicated efforts to determine stellar
properties (like Anders et al. 2022; Berger et al. 2023), or, e.g., provide
the level of accuracy required for the characterisation of exoplanet host
stars.

We expect that the sample of 1500 bona fide single stars in the Hyades
(Brandner et al. 2023a) and Pleiades open clusters with homogeneous
astrophysical parameters estimates derived from PARSEC evolutionary models
might also be highly useful for future benchmarks of theoretical and
observational studies.

###### Acknowledgements.

This work has made use of data from the European Space Agency (ESA) mission
Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing
and Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement. We acknowledge the use of
TOPCAT (Taylor 2005, 2011).

## References

  * Anders et al. (2022) Anders, F., Khalatyan, A., Queiroz, A. B. A., et al. 2022, A&A, 658, A91
  * Andrae et al. (2022) Andrae, R., Fouesneau, M., Sordo, R., et al. 2022, arXiv e-prints, arXiv:2206.06138
  * Bailer-Jones et al. (2013) Bailer-Jones, C. A. L., Andrae, R., Arcay, B., et al. 2013, A&A, 559, A74
  * Bailer-Jones et al. (2021) Bailer-Jones, C. A. L., Rybizki, J., Fouesneau, M., Demleitner, M., & Andrae, R. 2021, AJ, 161, 147
  * Baraffe & Chabrier (2018) Baraffe, I. & Chabrier, G. 2018, A&A, 619, A177
  * Basri et al. (1996) Basri, G., Marcy, G. W., & Graham, J. R. 1996, ApJ, 458, 600
  * Berger et al. (2023) Berger, T. A., Schlieder, J. E., & Huber, D. 2023, arXiv e-prints, arXiv:2301.11338
  * Bouy et al. (2013) Bouy, H., Bertin, E., Moraux, E., et al. 2013, A&A, 554, A101
  * Bouy et al. (2015) Bouy, H., Bertin, E., Sarro, L. M., et al. 2015, A&A, 577, A148
  * Bovy (2016) Bovy, J. 2016, ApJ, 817, 49
  * Brandner et al. (2023a) Brandner, W., Calissendorff, P., & Kopytova, T. 2023a, AJ, 165, 108
  * Brandner et al. (2023b) Brandner, W., Calissendorff, P., & Kopytova, T. 2023b, MNRAS, 518, 662
  * Brandt & Huang (2015) Brandt, T. D. & Huang, C. X. 2015, ApJ, 807, 58
  * Bressan et al. (2012) Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127
  * Chen et al. (2015) Chen, Y., Bressan, A., Girardi, L., et al. 2015, MNRAS, 452, 1068
  * Chen et al. (2014) Chen, Y., Girardi, L., Bressan, A., et al. 2014, MNRAS, 444, 2525
  * Choi et al. (2016) Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102
  * Creevey et al. (2022) Creevey, O. L., Sordo, R., Pailler, F., et al. 2022, arXiv e-prints, arXiv:2206.05864
  * Cutri et al. (2003) Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, VizieR Online Data Catalog, II/246
  * Dahm (2015) Dahm, S. E. 2015, ApJ, 813, 108
  * Delchambre et al. (2022) Delchambre, L., Bailer-Jones, C. A. L., Bellas-Velidis, I., et al. 2022, arXiv e-prints, arXiv:2206.06710
  * Fouesneau et al. (2022) Fouesneau, M., Frémat, Y., Andrae, R., et al. 2022, arXiv e-prints, arXiv:2206.05992
  * Funayama et al. (2009) Funayama, H., Itoh, Y., Oasa, Y., et al. 2009, PASJ, 61, 931
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Babusiaux, C., van Leeuwen, F., et al. 2018, A&A, 616, A10
  * Gaia Collaboration et al. (2021a) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021a, A&A, 649, A1
  * Gaia Collaboration et al. (2021b) Gaia Collaboration, Smart, R. L., Sarro, L. M., et al. 2021b, A&A, 649, A6
  * Galindo-Guil et al. (2022) Galindo-Guil, F. J., Barrado, D., Bouy, H., et al. 2022, A&A, 664, A70
  * Gossage et al. (2018) Gossage, S., Conroy, C., Dotter, A., et al. 2018, ApJ, 863, 67
  * Heyl et al. (2022) Heyl, J., Caiazzo, I., & Richer, H. B. 2022, ApJ, 926, 132
  * Hidalgo et al. (2018) Hidalgo, S. L., Pietrinferni, A., Cassisi, S., et al. 2018, ApJ, 856, 125
  * Høg et al. (2000) Høg, E., Fabricius, C., Makarov, V. V., et al. 2000, A&A, 355, L27
  * Jeffries (2017) Jeffries, R. D. 2017, Mem. Soc. Astron. Italiana, 88, 637
  * Kopytova et al. (2016) Kopytova, T. G., Brandner, W., Tognelli, E., et al. 2016, A&A, 585, A7
  * Krumholz & McKee (2020) Krumholz, M. R. & McKee, C. F. 2020, MNRAS, 494, 624
  * Magrini et al. (2022) Magrini, L., Danielski, C., Bossini, D., et al. 2022, A&A, 663, A161
  * Magrini et al. (2017) Magrini, L., Randich, S., Kordopatis, G., et al. 2017, A&A, 603, A2
  * Mermilliod (1992) Mermilliod, J. C. 1992, Bulletin d’Information du Centre de Donnees Stellaires, 40, 115
  * Meynet et al. (1993) Meynet, G., Mermilliod, J. C., & Maeder, A. 1993, A&AS, 98, 477
  * Parmentier et al. (2014) Parmentier, G., Pfalzner, S., & Grebel, E. K. 2014, ApJ, 791, 132
  * Pastorelli et al. (2020) Pastorelli, G., Marigo, P., Girardi, L., et al. 2020, MNRAS, 498, 3283
  * Percival et al. (2005) Percival, S. M., Salaris, M., & Groenewegen, M. A. T. 2005, A&A, 429, 887
  * Recio-Blanco et al. (2022) Recio-Blanco, A., de Laverny, P., Palicio, P. A., et al. 2022, arXiv e-prints, arXiv:2206.05541
  * Robin et al. (2003) Robin, A. C., Reylé, C., Derrière, S., & Picaud, S. 2003, A&A, 409, 523
  * Röser et al. (2011) Röser, S., Schilbach, E., Piskunov, A. E., Kharchenko, N. V., & Scholz, R. D. 2011, A&A, 531, A92
  * Rybizki et al. (2020) Rybizki, J., Demleitner, M., Bailer-Jones, C., et al. 2020, PASP, 132, 074501
  * Soderblom et al. (2009) Soderblom, D. R., Laskar, T., Valenti, J. A., Stauffer, J. R., & Rebull, L. M. 2009, AJ, 138, 1292
  * Stauffer et al. (1998) Stauffer, J. R., Schultz, G., & Kirkpatrick, J. D. 1998, ApJ, 499, L199
  * Tang et al. (2014) Tang, J., Bressan, A., Rosenfield, P., et al. 2014, MNRAS, 445, 4287
  * Taylor (2011) Taylor, M. 2011, TOPCAT: Tool for OPerations on Catalogues And Tables, Astrophysics Source Code Library, record ascl:1101.010
  * Taylor (2005) Taylor, M. B. 2005, in Astronomical Society of the Pacific Conference Series, Vol. 347, Astronomical Data Analysis Software and Systems XIV, ed. P. Shopbell, M. Britton, & R. Ebert, 29
  * van Leeuwen (2009) van Leeuwen, F. 2009, A&A, 497, 209