# Quarkonium production in A-A and p-A collisions

E. Scomparin  Istituto Nazionale di Fisica Nucleare – Sezione di Torino, Via
Giuria 1, 10125 Torino, Italy

###### Abstract

Thirty years ago, the suppression of quarkonium production in heavy-ion
collisions was first proposed as an unambiguous signature for the formation of
a Quark-Gluon Plasma. Recent results from the LHC run 2 have led to an
unprecedented level of precision on this observable and, together with new
data from RHIC, are providing an accurate picture of the influence of the
medium created in nuclear collisions on the various charmonium (J/$\psi$,
$\psi$(2S)) and bottomonium ($\Upsilon(1S)$, $\Upsilon(2S)$, $\Upsilon(3S)$)
states, studied via their decay into lepton pairs. In this contribution, I
will review the new results presented at Quark Matter 2017, emphasizing their
relation with previous experimental observations and comparing them, where
possible, with theoretical calculations.

###### keywords:

Quarkonium , Quark-Gluon Plasma

##  1 Introduction

In August 1987, at the 6th Quark Matter Conference in Nordkirchen (Germany),
the NA38 Collaboration presented the first evidence for a suppression of the
J/$\psi$ in O-U collisions at 200 GeV/nucleon [1]. By showing that the ratio
between the number of J/$\psi$ decaying to muon pairs and the corresponding
number of continuum events in the same mass region was decreasing by $\sim$35%
between peripheral and central events, charmonium suppression became the
“smoking gun” for the observation of deconfinement, a status retained to a
considerable extent also in more recent times. However, it became clear rather
soon that the interpretation of the experimental observations was far from
being trivial, due to competing suppression sources, which included cold
nuclear matter effects and dissociation processes in a hot hadronic gas.

After thirty years, and several rounds of experiments, starting with the fixed
target program at the CERN SPS and then moving to colliders (RHIC and LHC), a
consensus has been reached on the qualitative interpretation of the charmonium
suppression results in terms of a superposition of color screening and
recombination mechanisms. The first conceptually derives from the Matsui and
Satz original prediction [2], while the recombination effect [3, 4], which
tends to counterbalance color screening, becomes evident at the LHC energy
where the multiplicity of $c\overline{c}$ pairs is very large ($>100$ for
central Pb-Pb collisions). The most evident manifestation of the strong
contribution of recombination processes is visible in the comparison of the
J/$\psi$ nuclear modification factors ($R_{\rm AA}$) between RHIC and LHC
energy, with systematically larger $R_{\rm AA}$ values at the LHC, and in
particular for low transverse momentum ($p_{\rm T}$) values [5]. For
bottomonium, even at the LHC the recombination effects are expected to be
marginal, and a hyerarchy of suppression was observed in run 1 LHC results,
with smaller $R_{\rm AA}$ values for the more weakly bound $\Upsilon$ states
[6].

Although the aforementioned mechanisms can qualitatively explain the
observations, a satisfactory quantitative description is far from being
trivial. In particular, aspects such as the spectral modifications of
quarkonium states in hot matter and their formation time compared to that of
QGP play an important role and are presently much debated in the theory
community. On the experimental side, it is clear that the availability of more
accurate data, in particular for what concerns excited quarkonium states
and/or differential studies of $R_{\rm AA}$ can help to constrain the model
calculations. Connected to that, a precise determination of the open charm
cross section, to calibrate the size of the recombination process, is also
needed.

In the next Sections I will review the recent experimental achievements on
quarkonium production in nuclear collisions, with an emphasis on results
presented at the Quark Matter 2017 Conference and discussing a personal choice
of the main highlights.

##  2 Charmonium production: A-A collisions

In Pb-Pb collisions, inclusive J/$\psi$ production results were presented by
the ALICE Collaboration, covering the forward ($2.5<y<4$) and central
($|y|<0.9$) rapidity regions [7, 8]. Since the kinematic coverage extends to
$p_{\rm T}=0$, the centrality dependence of $R_{\rm AA}$ is dominated by
low-$p_{\rm T}$ J/$\psi$, where recombination effects are expected to be
stronger. In Fig. 1 (left) the preliminary results for the inclusive J/$\psi$
$R_{\rm AA}$ at central rapidity in Pb-Pb collisions at $\sqrt{s_{\rm
NN}}=5.02$ TeV [8] are compared with the recently published forward-$y$
results [9], as a function of centrality. The latter, thanks to the very fine
centrality binning, show a clear saturation of the J/$\psi$ suppression for
values of the number of participant nucleons $N_{\rm part}>100$. The
midrapidity result exhibits the same feature with a hint for a weaker
suppression for very central events. When comparing, as in Fig. 1 (right), the
forward-$y$ results with the corresponding ones from the LHC run 1
($\sqrt{s_{\rm NN}}=2.76$ TeV) [5], agreement is found within uncertainties,
an observation likely indicating that any additional suppression due to the
larger energy density at higher collision energy is counterbalanced by a
stronger recombination effect. The stronger suppression visible in PHENIX data
[10] clearly implies that recombination plays a less important role at the
lower collision energy.

Fig. 1: Left: centrality dependence of the inclusive J/$\psi$ $R_{\rm AA}$ for
Pb-Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV, measured by ALICE, in two
$y$ regions [8, 9]. Right: inclusive J/$\psi$ $R_{\rm AA}$ as a function of
centrality for various $\sqrt{s_{\rm NN}}$ values, from ALICE and PHENIX [5,
9, 10].

Thanks to the large integrated luminosity for the $\sqrt{s_{\rm NN}}=5.02$ TeV
data sample, it is now possible to study the centrality dependence of the
$R_{\rm AA}$ differentially as a function of the J/$\psi$ kinematic variables
[7]. In Fig. 2 (left), a new result from ALICE on the nuclear modification
factor is presented, for four $p_{\rm T}$ bins, as a function of $\langle
N_{\rm part}\rangle$. While for peripheral collisions no significant
dependence on $p_{\rm T}$ is visible, for central collisions a strong
reduction of the suppression can be seen at low $p_{\rm T}$, making the
centrality dependence of $R_{\rm AA}$ in $0.3<p_{\rm T}<2$ GeV/$c$ almost
flat. The comparisons with a theoretical model [11] shows that the main
features of the data can be reproduced, but some tension persists, e.g., at
intermediate $p_{\rm T}$ and centrality. A systematic comparison of such
multi-differential $R_{\rm AA}$ measurements with theoretical calculations
represents a promising way to pin down quantitatively the features of the
various mechanisms at play.

Fig. 2: Left: centrality dependence of the inclusive J/$\psi$ $R_{\rm AA}$ for
Pb-Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV, measured by ALICE for
$2.5<y<4$ and in four $p_{\rm T}$ bins [7], compared with the result of a
theory calculation [11]. Right: the $p_{\rm T}$-dependence of inclusive
J/$\psi$ $v_{2}$ for 20-40% Pb-Pb collisions [7], measured by ALICE in two
rapidity ranges, and compared with the D-meson $v_{2}$ for 30-50% centrality
[12].

As a further step forward, the first observation of a non-zero elliptic flow
($v_{2}$) for the J/$\psi$ in Pb-Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV
was presented by ALICE at the Conference [7]. In Fig. 2 (right), the $p_{\rm
T}$ dependence of $v_{2}$ is shown for the 20-40% centrality interval, in the
ranges $|y|<0.9$ and $2.5<y<4$. A significance up to 7.4$\sigma$ for a non-
zero $v_{2}$ was estimated for the forward-$y$ data. This result may indicate
that a significant fraction of the observed J/$\psi$ comes from charm quarks
which thermalized in the QGP. The comparison, also shown in Fig. 2 (right),
with new results on the D-meson elliptic flow in the 30-50% centrality range
[12], opens up the possibility of investigating the relation between light and
heavy quark flow in the medium.

The study of the weakly bound $\psi(2S)$ represents a potentially important
tool, as one would expect much stronger dissociation effects in the QGP, while
the evaluation of recombination is less straightforward, as it might become
important at a later stage in the collision history. Recent results from CMS
[13, 14], shown in Fig. 3 (left), indicate for Pb-Pb collisions at
$\sqrt{s_{\rm NN}}=5.02$ TeV, a stronger suppression for the $\psi(2S)$ with
respect to the J/$\psi$, quantified through the measurement of the double
ratio of $\psi(2S)$ to J/$\psi$ cross sections between Pb-Pb and pp. This
result, which refers to the region $1.6<|y|<2.4$, $3<p_{\rm T}<30$ GeV/$c$,
shows a different trend with respect to the corresponding one obtained at
$\sqrt{s_{\rm NN}}=2.76$ TeV [15], where a strong enhancement of the double
ratio was observed, in particular for central events. Preliminary ALICE
results [7], covering the region down to zero $p_{\rm T}$, were shown at the
Conference (Fig. 3 (right)). In particular, for $p_{\rm T}<3$ GeV/$c$ and
$2.5<y<4$, values of the double ratio smaller than 1 were measured. The only
proposed mechanism for the enhancement shown by CMS at $\sqrt{s_{\rm
NN}}=2.76$ TeV involves a $\psi(2S)$ regeneration mechanism occurring late in
the collision history, when radial flow is already built-up and therefore the
charmonia have acquired a significant transverse momentum [16]. Also in this
case, it will be important to provide systematic comparisons between data and
calculations to validate this picture.

Fig. 3: Left: centrality dependence of the double ratio of $\psi(2S)$ and
J/$\psi$ cross sections between Pb-Pb and pp, measured by CMS at $\sqrt{s_{\rm
NN}}=2.76$ and 5.02 TeV for $3<p_{\rm T}<30$ GeV/$c$ [14]. Right: same
quantity, measured by ALICE for $p_{\rm T}<3$ GeV/$c$ [7].

##  3 Charmonium production: p-A collisions

In the context of charmonium production, the study of p-A collisions is
relevant for the interpretation of A-A results, since cold nuclear matter
effects, which can be studied via such collisions and are not related to the
production of a QGP, are known to be sizeable. In particular, parton shadowing
in the nucleus, coherent parton energy loss and saturation phenomena related
to the formation of a Color Glass Condensate (CGC) are also interesting from a
more general point of view, to investigate aspects of QCD that are up to now
not completely settled. Studies of p-Pb collisions during the LHC run 1
($\sqrt{s_{\rm NN}}=5.02$ TeV) have shown that J/$\psi$ nuclear modification
factors as low as $R_{\rm pPb}\sim 0.6$ can be reached at forward rapidity
($y_{\rm cms}\sim 3.5$) and low $p_{\rm T}$, a value similar to the one
observed in Pb-Pb [17]. In other words, this implies that by correcting the
measured $R_{\rm AA}$ for cold nuclear matter effects, one would end up with
$R_{\rm AA}\sim 1$ at low $p_{\rm T}$, indicating that suppression and
recombination effects for low-$p_{\rm T}$ J/$\psi$ in Pb-Pb are of a similar
size.

With the very recent availability of high statistics p-Pb data from the LHC
run 2 ($\sqrt{s_{\rm NN}}=8.16$ TeV), the possibility of producing high-
quality multi-differential results has opened up. For the moment, first
results on the J/$\psi$ $R_{\rm pPb}$ were presented by ALICE [7, 18] as a
function of rapidity and $p_{\rm T}$ (see Fig. 4). They are compatible with
previous results at $\sqrt{s_{\rm NN}}=5.02$ TeV (not shown [17]) and are in
good agreement with theory models which include various cold nuclear matter
effects. In addition, the $p_{\rm T}$ coverage was now extended from 8 up to
20 GeV/$c$. These results indicate that the effect of J/$\psi$ break-up in
nuclear matter (not included in the models) is weak. Such effect was known to
be non-negligible at lower energy, where the $c\overline{c}$ pair spends more
time in nuclear matter. This fact has received another confirmation by
preliminary STAR data [19] in p-Au collisions shown at the Conference (Fig. 5
(left)), which indicate that only including a J/$\psi$ break-up contribution
$\sigma_{\rm abs}=4.2$ mb in the theoretical models a good agreement can be
found.

Fig. 4: Rapidity (left) and $p_{\rm T}$ (right) dependence of the inclusive
J/$\psi$ $R_{\rm pPb}$ for p-Pb collisions at $\sqrt{s_{\rm NN}}=8.16$ TeV,
shown by ALICE [7, 18], and compared to theory calculations (see [18] for the
corresponding references).

Concerning the $\psi(2S)$, its weaker binding energy may lead to a different
behaviour with respect to J/$\psi$. In principle, one would expect at the LHC
energies no significant effect from break-up in nuclear matter also for the
$\psi(2S)$, because of the short crossing times in the nucleus for the
$c\overline{c}$ pair. However, results at both RHIC and LHC energies [20, 21,
22], shown in Fig. 5 (right), indicate that the double ratio between
$\psi(2S)$ and J/$\psi$ in p-Pb and pp is significantly smaller than 1 and
decreases for more central collisions. The only explanation for this
observation is a suppression effect due to the interaction with the system
produced in the collision, which can be a dense hadronic system or even a
deconfined medium, if the energy density reached in the collision is large
enough. The results from the high statistics p-Pb data from the LHC run 2 will
likely help in constraining the models describing this suppression process.

Fig. 5: Left: $p_{\rm T}$ dependence of the J/$\psi$ $R_{\rm pAu}$, shown by
STAR [19] and compared with the results of theory calculations. Right: double
ratio of $\psi(2S)$ and J/$\psi$ cross sections between p-A (or 3He–Au) and
pp, as a function of the comoving particle density, estimated starting from
the collision centrality by using AMPT simulations [20].

##  4 Bottomonium production: A-A collisions

In recent years, the study of bottomonia has received considerable attention.
The relatively small multiplicity of $b$ quarks, even at the LHC energies,
makes the recombination effects rather weak, allowing an easier isolation of
phenomena related to color screening. In addition, the $\Upsilon(1S)$, with a
binding energy larger than 1 GeV/$c$, represents a good probe of the very hot
QGP formed at the highest available collision energies. Furthermore, the
simultaneous study of the 1S, 2S and 3S states, which have increasingly weaker
binding energy ($\sim 0.2$ GeV for the 3S) is expected to provide a good
sensitivity to the characteristics of the medium. Exploratory studies
performed at RHIC [23] have given hints for a significant suppressions of the
$\Upsilon$ states, while results from the LHC run 1, from ALICE [24] and in
particular from CMS [6], have established a clear hyerarchy in the suppression
in Pb-Pb collisions, with $R_{\rm AA}(\Upsilon(1S))\sim 0.45$, the more weakly
bound 2S state having $R_{\rm AA}\sim$ 0.1, and the 3S state no longer
visible.

At this Conference, a certain number of new results from run 2 were presented
for the first time. In Fig. 6 (left) the centrality dependence of the nuclear
modification factor for the various $\Upsilon$ states, measured by CMS [25,
26] at $\sqrt{s_{\rm NN}}=5.02$ TeV, is shown, while in Fig. 6 (right) a
comparison of the $\Upsilon(1S)$ $R_{\rm AA}$ at $\sqrt{s_{\rm NN}}=2.76$ and
5.02 TeV is reported [25, 6]. The very strong suppression of the 2S and 3S
states confirms the observation carried out by CMS on run 1 data. A slightly
stronger suppression at the higher energy is visible for the $\Upsilon(1S)$,
but it is not significant with respect to the measurement uncertainties. The
observed $R_{\rm AA}\sim 0.3$ for the 1S state in central events points to a
strong effect also on this tightly bound resonance. However, it is known that
due to the feed-down to the 1S state from higher mass (and more weakly bound)
states, it is not straightforward to conclude on the suppression of the
directly produced $\Upsilon(1S)$ (see discussion in the next Section).

Fig. 6: Left: centrality dependence of the $R_{\rm AA}$ for the various
$\Upsilon$ resonances in Pb-Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV,
shown by CMS [25, 26]. Right: comparison of the $\Upsilon(1S)$ $R_{\rm AA}$,
as a function of centrality, at $\sqrt{s_{\rm NN}}=5.02$ and 2.76 TeV [25,
26].

Considering together bottomonium results from CMS and ALICE, a good coverage
in the region $0<y_{\rm cms}<4$ is obtained. In Fig. 7 the rapidity dependence
of the $\Upsilon(1S)$ $R_{\rm AA}$ is shown at the two energies that have been
studied at the LHC [25, 27]. While at $\sqrt{s_{\rm NN}}=5.02$ TeV an overall
flat distribution can be seen, hints for a stronger suppression at large
$y_{\rm cms}$ are visible at $\sqrt{s_{\rm NN}}=2.76$ TeV. The latter
behaviour is typical in case of the presence of recombination effects, where
the higher multiplicity of heavy quarks at central rapidity would induce a
stronger regeneration of bottomonia. However, the non-negligible size of
systematic uncertainties induces to consider this possibility with caution,
also because the recombination effects should in principle be stronger at the
higher collision energy. The comparison of the results with a theoretical
model based on a hydrodynamical description of the plasma [28] shows the
opposite effect, with a weaker suppression at large $y$, which clearly leads
to a tension with the experimental results at $\sqrt{s_{\rm NN}}=2.76$ TeV.

Fig. 7: The rapidity dependence of the $\Upsilon(1S)$ $R_{\rm AA}$ at
$\sqrt{s_{\rm NN}}=2.76$ (left) and 5.02 TeV (right) [25, 26]. Compilation of
data from ALICE and CMS, compared with the result of a theoretical model based
on an anisotropic hydrodynamic description of the QGP [27, 25, 28].

In addition to the rich sample of results from the LHC experiments, new
bottomonium results from Au-Au collisions at $\sqrt{s_{\rm NN}}=0.2$ TeV were
presented by STAR at the Conference, separately for the $\Upsilon(1S)$ and for
the sum of the 2S and 3S states [29]. For the $\Upsilon(2S+3S)$, the $R_{\rm
AA}$, compared to CMS results for the two resonances [6], shows an indication
for a weaker suppression at RHIC energies up to $N_{\rm part}\sim 200$, while
for central collisions the suppression effects become quite similar.
Qualitatively, this observation could be connected with the increase of the
energy density with centrality, followed by a saturation of the suppression
effect from central collisions at RHIC up to LHC energy. For the
$\Upsilon(1S)$ (Fig. 8 (left)), the centrality dependence of $R_{\rm AA}$ is
very similar for STAR and CMS results, an observation which could be explained
if the suppression were mainly due to the feed-down effect from $\Upsilon(2S)$
and $\Upsilon(3S)$, with little or no effect on the direct $\Upsilon(1S)$ at
both RHIC and LHC. However, to confirm this scenario, an accurate knowledge of
cold nuclear matter effects and of feed-down contributions from higher mass
resonances, together with quantitative model comparisons, are necessary.

Fig. 8:  Left: preliminary results on the centrality dependence of the
$\Upsilon(1S)$ $R_{\rm AA}$, measured by STAR in Au-Au collisions at
$\sqrt{s_{\rm NN}}=0.2$ TeV [29] and compared with the corresponding CMS
results (stars) [6]. Right: rapidity dependence of $R_{\rm pPb}$ for
$\Upsilon(1S)$ in p-Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV, with
results from ATLAS, LHCb and ALICE [30, 31, 32].

##  5 Bottomonium production: p-A collisions

Similarly to the charmonium case, also for bottomonia cold nuclear matter
effects can in principle be sizeable, and need to be estimated via proton-
nucleus collisions. For the moment, only the LHC run 1 results are available
and their quality, as can be seen in Fig. 8 (right) is still lower than for
charmonia. In particular, ATLAS mid-rapidity results in p-Pb collisions at
$\sqrt{s_{\rm NN}}=5.02$ TeV [30] favour $R_{\rm pPb}$ values close to 1,
while at forward (p-going) and backward (Pb-going) rapidity, LHCb [31] and
ALICE [32] results, although compatible, span a rather large range of possible
values. The trend from small (or no) cold nuclear matter effects at negative
$y$ towards a suppression at positive $y$, clearly observed for J/$\psi$,
seems visible also for the $\Upsilon(1S)$. Comparisons with theory models (not
shown) are still not conclusive, due to the large uncertainties of the
experimental results. From a qualitative point of view, one may attempt to
make use of p-Pb data to evaluate the presence of a suppression signal for
directly produced $\Upsilon(1S)$ in Pb-Pb. Restricting these considerations to
the region around mid-rapidity, the ATLAS p-Pb data allow, at 1$\sigma$ level,
$R_{\rm pPb}\sim 0.8$, which, assuming a dominance of shadowing effects, would
imply for Pb-Pb collisions, $R_{\rm PbPb}^{\rm CNM}\sim 0.8^{2}=0.64$. On the
other hand, CMS measures, in Pb-Pb collisions at the same energy, for
inclusive $\Upsilon(1S)$ production, $R_{\rm PbPb}^{\rm incl}=0.36$ [26].
Recent estimates of feed-down effects, based mostly on LHCb data (see e.g.,
[33]), indicate a $\sim 30$% effect, which would lead to $R_{\rm PbPb}^{\rm
direct}\sim 0.36/0.7=0.51$, a value smaller than CNM extrapolations. Even if
this evaluation is admittedly qualitative, the conclusion that can be drawn is
that we are likely observing a suppression of direct $\Upsilon(1S)$, which,
according to lattice calculations, would imply a QGP temperature at least
$~{}\sim 2T_{\rm c}$ at the LHC energies [34].

Finally, also for $\Upsilon(2S)$ and $\Upsilon(3S)$ no precise evaluations of
$R_{\rm pPb}$ are available for the moment. A relative measurement (double
ratio between the excited and the 1S state between p-Pb and pp) was published
by CMS [35]. It clearly indicates that in p-Pb collisions the 2S and 3S states
are more suppressed than the $\Upsilon(1S)$. As for the $\psi(2S)$ vs J/$\psi$
case, one does not expect a different $R_{\rm pPb}$ between the various states
due to either shadowing or interactions with nuclear matter, and the only
mechanism that may explain this observation is the break-up of the 2S and 3S
states via interactions with the medium (either partonic or hadronic) produced
in the p-Pb collision.

##  6 Conclusions

After many years, the study of quarkonium production continues to be a source
of new and relevant results for our understanding of the QGP created in heavy-
ion collisions. Several highlights were presented at Quark Matter 2017: for
charmonium production the first unambiguous observation of a strong elliptic
flow for the J/$\psi$ at low $p_{\rm T}$ has been established, and precise new
data on J/$\psi$ suppression and regeneration have become available, thanks to
the first analyses of the LHC run 2 data in Pb-Pb collisions. Also the
$\psi(2S)$ results are becoming more accurate, and even if their $\sqrt{s_{\rm
NN}}$-dependence is still not completetly understood, they offer the
possibility of a deeper insight on charmonium properties in the medium.
Concerning bottomonia, detailed information on $\Upsilon(1S)$ and
$\Upsilon(2S)$ $R_{\rm AA}$ is now available at both $\sqrt{s_{\rm NN}}=2.76$
and 5.02 TeV, and an unambiguous evidence for a hierarchy of suppression was
reached. The description of the $y$-dependence of $R_{\rm AA}$ is still not
fully satisfactory, and might need some further refinement in the models.
Finally, a first set of precise results from RHIC has become available, and
will surely represent an important ingredient towards a unified description of
the bottomonium phenomenology from low to high energy.

## References

  * [1] M. C. Abreu, et al., Z. Phys. C38 (1988) 117.
  * [2] T. Matsui, H. Satz, Phys. Lett. B178 (1986) 416.
  * [3] P. Braun-Munzinger, J. Stachel, Phys. Lett. B490 (2000) 196–202.
  * [4] R. L. Thews, M. Schroedter, J. Rafelski, Phys. Rev. C63 (2001) 054905.  arXiv:hep-ph/0007323.
  * [5] J. Adam, et al., JHEP 05 (2016) 179.  arXiv:1506.08804.
  * [6] V. Khachatryan, et al., Phys. Lett. B.  arXiv:1611.01510.
  * [7] M. Tarhini, Measurement of charmonium production at forward rapidity in Pb-Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV with ALICE, QM2017 proceedings.
  * [8] T. J. Bustamante, ALICE measurement of the J/$\psi$ nuclear modification factor at mid-rapidity in Pb-Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV, QM2017 proceedings.
  * [9] J. Adam, et al., Phys. Lett. B766 (2017) 212–224.  arXiv:1606.08197.
  * [10] A. Adare, et al., Phys. Rev. C84 (2011) 054912.  arXiv:1103.6269.
  * [11] X. Zhao, R. Rapp, Nucl. Phys. A859 (2011) 114–125.  arXiv:1102.2194.
  * [12] A. Barbano, D-meson nuclear modification factor and elliptic flow measurements in Pb-Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV with ALICE at the LHC, QM2017 proceedings.
  * [13] J. M. Blanco, $\psi(2S)$ and J/$\psi$ modification in p-Pb and Pb-Pb collisions at 5.02 TeV with CMS, QM2017 proceedings.
  * [14] A. M. Sirunyan, et al., Submitted to: Phys. Rev. Lett.arXiv:1611.01438.
  * [15] V. Khachatryan, et al., Phys. Rev. Lett. 113 (26) (2014) 262301.  arXiv:1410.1804.
  * [16] X. Du, R. Rapp, J. Phys. Conf. Ser. 779 (1) (2017) 012042.  arXiv:1609.04868.
  * [17] J. Adam, et al., JHEP 06 (2015) 055.  arXiv:1503.07179.
  * [18] Preliminary Physics Summary: Inclusive J/$\psi$ production at forward rapidity in p-Pb collisions at $\sqrt{s}_{\rm{NN}}=8.16$ TeV, Tech. Rep. ALICE-PUBLIC-2017-001, CERN, Geneva (Feb 2017).
  * [19] T. Todoroki, Measurements of charmonium production in pp, p-Au and Au-Au collisions at $\sqrt{s_{\rm NN}}=200$ GeV with the STAR experiment, QM2017 proceedings.
  * [20] S. Lim, PHENIX measurements of open and hidden heavy-flavor in pp, p-Al and p/d/3He-Au collisions at $\sqrt{s_{\rm NN}}=200$ GeV across a wide range of rapidity, QM2017 proceedings.
  * [21] A. Adare, et al., Phys. Rev. C95 (3) (2017) 034904.  arXiv:1609.06550.
  * [22] B. Abelev, et al., JHEP 12 (2014) 073.  arXiv:1405.3796.
  * [23] L. Adamczyk, et al., Phys. Lett. B735 (2014) 127–137, [Erratum: Phys. Lett.B743,537(2015)].  arXiv:1312.3675.
  * [24] B. Abelev, et al., Phys. Lett. B738 (2014) 361–372.  arXiv:1405.4493.
  * [25] C. Flores, Bottomonia results from the LHC Run 1 and 2 with CMS, QM2017 proceedings.
  * [26] Measurement of Nuclear Modification Factors of $\Upsilon\textrm{(nS)}$ Mesons in PbPb Collisions at $\sqrt{s_{{}_{\rm{NN}}}}=5.02~{}\mathrm{TeV}$, Tech. Rep. CMS-PAS-HIN-16-023, CERN, Geneva (2017).
  * [27] I. Das, Bottomonium measurements at forward $y$ in Pb-Pb and p-Pb collisions with ALICE at the LHC, QM2017 proceedings.
  * [28] B. Krouppa, M. Strickland, Universe 2 (3) (2016) 16.  arXiv:1605.03561.
  * [29] Z. Ye, $\Upsilon$ measurements in p-p, p-Au and Au-Au at $\sqrt{s_{\rm NN}}=200$ GeV with the STAR experiment, QM2017 proceedings.
  * [30] Measurement of $\Upsilon(\mathrm{nS})$ production with $p$+Pb collisions at $\sqrt{s_{\rm NN}}=5.02~{}\mathrm{TeV}$ and $pp$ collisions at $\sqrt{s}=2.76~{}\mathrm{TeV}$, Tech. Rep. ATLAS-CONF-2015-050, CERN, Geneva (Sep 2015).
  * [31] R. Aaij, et al., JHEP 1407 (2014) 094.  arXiv:1405.5152.
  * [32] B. B. Abelev, et al., Phys. Lett. B740 (2015) 105–117.  arXiv:1410.2234.
  * [33] R. Aaij, et al., JHEP 11 (2012) 031.  arXiv:1209.0282.
  * [34] Y. Burnier, O. Kaczmarek, A. Rothkopf, JHEP 12 (2015) 101.  arXiv:1509.07366.
  * [35] S. Chatrchyan, et al., JHEP 04 (2014) 103.  arXiv:1312.6300.


# Formation of ultracold ion pairs through long-range Rydberg molecules

Michael Peper  Laboratory of Physical Chemistry, ETH Zürich, 8093 Zürich,
Switzerland  Johannes Deiglmayr  Felix Bloch Institute, University of Leipzig,
04103 Leipzig, Germany<EMAIL_ADDRESS>
###### Abstract

We propose a new approach to excite ion-pair states of ultracold dimers. The
central idea is a two-step process where first long-range Rydberg molecules
are formed by photoassociation, which are then driven by stimulated emission
towards the ion-pair state, a process bearing features of a photo-induced
harpooning reaction. We assess the feasibility of this approach through a
detailed experimental and theoretical study on a specific system, $p$-wave-
scattering dominated long-range Rydberg molecules in caesium, and discuss
potential applications for the study of strongly correlated plasmas consisting
of oppositely charged particles of equal or similar mass.

  * September 2019

##  1 Introduction

In recent years, a solid understanding of the properties of long-range Rydberg
molecules has been established and the developed theoretical tools explain
experimental observations reliably, often even quantitatively [1, 2]. Long-
range Rydberg molecules are bound states of an atom in a Rydberg state and a
ground-state atom within the orbit of the Rydberg electron. The binding
results from the Fermi-contact interaction [3] of the (almost) free Rydberg
electron scattering off the neutral ground-state atom. This interaction was
first observed as a pressure-dependent shift of the transition frequencies in
the spectroscopy of thermal vapours [4, 5]. In Fermi’s original pseudo-
potential derivation, later extended by Omont [6] to include higher partial
waves, the interaction potential between the Rydberg and the ground-state atom
is given (in atomic units) by the expression

$V=2\pi a_{s}\left|\Psi(R)\right|^{2}+6\pi
a_{p}\left|\nabla\Psi(R)\right|^{2}\dots\;,$ (1)

where $\left|\Psi(R)\right|^{2}$ is the probability density of the Rydberg
electron at the position $R$ of the ground-state atom, $\nabla\Psi(R)$ is the
gradient of the wavefunction of the Rydberg electron, and the $s$-wave
scattering length $a_{s}$ and the $p$-wave scattering volume $a_{p}$ are
parameters directly derived from the electron-neutral interaction potential.
The proportionality of the interaction strength to the probability density of
the Rydberg electron leads to oscillatory interaction potentials supporting
vibrational bound states, whose internuclear equilibrium separations are
comparable to the size of the classical Rydberg orbit [7].

We propose to exploit the properties of long-range Rydberg molecules for
exploring a new route to highly-excited ion-pair states. Spectroscopy of ion-
pair states and extrapolation of their Rydberg series to the ionization energy
can yield precise values for electron affinities [8, 9]. Dissociation of
weakly bound, ultracold ion pairs could provide a unique source of ultracold
anions, e.g., for sympathetic cooling of antihydrogen [10]. Ultracold clouds
of dissociated ion pairs would also open up the ways to study neutral plasmas
in an unexplored regime: in contrast to electron-cation plasmas studied
previously [11, 12], an anion-cation plasma, consisting of oppositely charged
components with identical or similar masses, would feature dynamics on
experimentally accessible time scales of microseconds [13].

The interaction potential between an anion and a cation at large distance is
dominated by their Coulomb interaction just like in the case of a Rydberg atom
where the negatively charged electron is bound to a positively charged ion
core. The binding energy of vibrational levels in the ion-pair potential is
given by the Rydberg formula

$E=-hc\frac{R_{\mu}}{n_{\mathrm{HR}}^{2}}\;,$ (2)

where $R_{\mu}=R_{\infty}\frac{\mu}{m_{\mathrm{el}}}$ is the mass-corrected
Rydberg constant, $n_{\mathrm{HR}}$ is an effective vibrational quantum
number, $R_{\infty}$ is Rydberg’s constant in wavenumbers, $\mu$ is the
reduced mass of the two bound particles, $m_{\mathrm{el}}$ is the rest mass of
the electron, and $h$ and $c$ are the Planck constant and the speed of light,
respectively. For a bound cation-anion system, the value of $R_{\mu}$ is much
larger than for the Rydberg states of atoms or molecules and the vibrational
levels in the ion-pair-state potential have thus been named heavy-Rydberg
states [14]. Because $R_{\mu}$ is very large, $n_{\mathrm{HR}}$ of a heavy-
Rydberg state is much larger than the principal quantum number $n$ of an
electronically-excited Rydberg state with the same binding energy.

At short internuclear distances, the electronic configuration of the ion-pair
state mixes with excited electronic configurations of the neutral-parent
dimer, possibly leading to strong perturbations of the electronic structure.
Previous excitation schemes of molecular heavy-Rydberg states relied on such
short-range mixings to excite the ion-pair states from the ground state of the
parent molecule in a supersonic molecular beam (see, e.g., [14, 15]) or in
ultracold gases [16, 17]. However, these schemes have not yet yielded dense
samples of heavy-Rydberg molecules, mostly because the exploited perturbations
are weak and their accurate prediction is a significant challenge for current
quantum-chemical methods.

Figure 1:  a) Schematic route for the excitation of heavy-Rydberg ion-pair
states: a $p$-wave dominated scattering state of a Rydberg electron and a
ground-state atom is deexcited to the electronic ground state of Cs- by
stimulated emission driven via a mid-IR laser pulse; b) electronic structure
of Cs-: the dipole-allowed transition from the initial 1P1 scattering state (a
shape-resonance) to the 1S0 electronic ground state is indicated by an arrow;
c) Coulomb potential of the Cs- \- Cs+ system and close-lying atomic
asymptotes of Cs2 [18], the red line marks the equilibrium separation $R\sim
1500\,a_{0}$ of a long-range Rydberg molecule in 352P3/2.

Here we propose to excite the ion-pair states of a dimer via molecular long-
range Rydberg states. In the past, we have studied extensively the long-range
Rydberg states of Cs2 [19, 20] and thus use them as model system here. In a
semiclassical picture, the binding of long-range Rydberg molecules is based on
the scattering of the almost-free Rydberg electron off a ground-state atom
within its orbit. By stimulating the radiative transition from the electron-
atom scattering state to the ground state of Cs-, as depicted in Fig. 1 a),
the long-range molecular Rydberg state is transferred into a bound state of
Cs- and Cs+. Note that in the case of a free electron, this transition would
correspond to the reverse of the photo-detachment transition observed at a
wavelength of $\sim 2.6\,\mu$m in Cs- (Fig. 1 b)) [21].

The selection rules for electric dipole transitions allow only for transitions
between the 1P1 $p$-wave component of the scattering state and the 1S0
electronic ground state of Cs- (see Fig. 1 b)). Thus, the initial long-range
Rydberg molecule has to be chosen such that its binding has a significant
contribution from the 1P1 scattering channel. The existence of bound states
with this characteristic, called “Butterfly molecules” because of the shape of
their spatial wavefunction, has been demonstrated both theoretically [22] and
experimentally [23]. Transitions from a molecular Butterfly state to an ion-
pair state may thus provide a novel route towards the excitation of heavy-
Rydberg states [23].

The Franck-Condon principle will strongly favour transitions to Cs-Cs+ ion-
pair states with an outer turning point close to the equilibrium distance of
the molecular long-range Rydberg state (typically around 1000 to 2000
$a_{0}$). The equilibrium distance of the ion-pair state is then about half of
this distance, as known from the general properties of Rydberg states [24].
The relevant part of the Coulomb-interaction potential is shown in Fig. 1 c)
together with the close-lying atomic dissociation asymptotes of Cs2.
Inspection of Fig. 1 c) reveals a region of internuclear distances $R$ from
$\sim 500\,a_{0}$ to $>2500\,a_{0}$ where the potential-energy curve of the
ion-pair state does not cross any dissociation threshold and should remain
unperturbed.

In this article we present a detailed numerical study of this system, using
Cs2 as an example, emphasizing the importance of singlet-triplet mixing [19]
in the molecular long-range Rydberg state. We will separately discuss the
following steps: i) photoassociation of a pair of colliding atoms into
metastable long-range Rydberg molecules with an electronic character that has
significant contributions from 1P1 scattering and ii) driving stimulated
emission from the long-range Rydberg to the ion-pair state. This photoinduced
charge-transfer process bears features of a harpoon reaction, where the
Rydberg atom captures the ground-state atom by donating its electron to the
neutral and a more strongly bound ion pair is formed. Note that the direct
photoassociation into an electronic state with ion-pair character is only
possible at much shorter internuclear distances than considered here, where
covalent and ionic molecular states are strongly mixed [16, 17].

##  2 Photoassociation of butterfly-like long-range Rydberg molecules

In 2002, Hamilton et al. predicted that shape resonances in $p$-wave electron-
alkali-atom scattering cause molecular long-range Rydberg states with large
binding energies and electronic wavefunctions resembling (in a specific
graphical representation) the shape of a butterfly [22]. Shape resonances
result from metastable scattering complexes bound behind a centrifugal barrier
which can strongly modify the scattering interactions of colliding particles.
On resonance, the scattering electron is temporarily bound in an excited state
of the anion which provides an opportunity to deexcite the collision complex
to the anionic ground state. Experimental signatures of butterfly molecules
have been observed in rubidium close to the 25 2P3/2 asymptote, where an
analysis of the rotationally-resolved photoassociation spectra revealed that
the excited molecular Rydberg states had internuclear separations between 100
and 300 $a_{0}$ [23]. More recently, molecular bound states with strong
contributions to their binding from $p$-wave scattering have been observed
close to $n$ 2S1/2 [25] and $n$ 2Dj [26] Rydberg states of rubidium. The
$p$-wave character of these “butterfly-like” states, bound just below low-$l$
Rydberg states, is smaller than the one of the true butterfly molecules
discussed above. However, the equilibrium distances of these molecules are on
the order of 1000 $a_{0}$, which makes such states a more promising
intermediate state for transferring an ultracold atomic gas into a gas of
weakly bound ion pairs.

Here we investigate $p$-wave-scattering dominated trilobite-like molecular
states in caesium close to $n$ 2P3/2 states. The electron-caesium scattering
features broad $p$-wave shape resonances close to zero collision energy [21],
and thus strong contributions from $p$-wave scattering to the binding of long-
range Rydberg molecules are to be expected. The experimental setup has been
described previously [19]: we release ultracold samples of Cs atoms from a
compressed magneto-optical trap (MOT) with typical densities of 1012 atoms/cm3
and temperatures of 40 $\mu$K, as determined by absorption imaging. The atoms
are then prepared in a specific hyperfine state ($F=3$ or 4) of the 6 2S1/2
ground state by optical pumping. After turning off all optical and magnetic
trapping fields, single-photon photoassociation into long-range Rydberg
molecules is performed. The UV radiation driving the photoassociation
transitions is obtained by frequency-doubling the output of a continuous-wave
(cw) single-mode ring dye laser at $\sim$639 nm and forming pulses of 3 $\mu$s
length using an acousto-optic modulator. The laser frequency is calibrated
with the help of a wavemeter and a frequency comb [27]. The formation of
Rydberg molecules is detected by applying pulsed electric potentials to a pair
of electrodes surrounding the photoexcitation region, creating an electric
field with sufficient magnitude to field ionize the Rydberg molecules. Field-
ionization of long-range Cs2 Rydberg molecules produces an electron, a Cs+
ion, and a neutral Cs atom, because the molecule instantaneously dissociates
upon removal of the Rydberg electron, which mediates the binding interaction.
We detect the resulting Cs+ ions by accelerating them towards a micro-channel
plate (MCP) detector and sample the detector current using a fast
oscilloscope. The extraction field also accelerates Cs${}_{2}^{+}$ ions, which
result from vibronic autoionization of long-range Rydberg molecules [28],
towards the detector. Cs+ and Cs${}_{2}^{+}$ ions are easily separated by
their different arrival times at the detector. The respective ion-yields,
plotted in Figure 2 a), are obtained by integrating the oscilloscope traces in
narrow integration windows around the mean arrival times of Cs+ and
Cs${}_{2}^{+}$ ions.

Figure 2:  Left panel: experimental photoassociation spectra in the vicinity
of the 62S1/2($F=4$)-$n^{2}$P3/2 asymptote, comparing the yield of Cs+ ions
resulting from field-ionization of Rydberg atoms and molecules (solid line) to
the yield of Cs${}_{2}^{+}$ ions resulting from associative ionization of
Rydberg molecules (dashed line). Right panel: calculated potential energy
curves with ${}^{1,3}\Sigma^{+}$ symmetry (solid lines) and ${}^{3}\Sigma^{+}$
symmetry (dashed lines). Asterisks mark experimental resonances (left panel)
arising from photoassociation of molecules localized in the outer-most wells
(right panel). Resonances attributed to photoassociation of $p$-wave dominated
molecules are marked in the left panel by coloured lines, where the lower
vertical line indicates the experimental line position and the upper vertical
line the calculated binding energy of the assigned molecular state. The
calculated binding energies are also indicated in the right panel by
horizontal lines.

In Fig. 2 (left panel), a photoassociation spectrum taken in the vicinity of
the 35 2P3/2 asymptote is presented. The right-hand panel of this figure
depicts the molecular structure close to this asymptote, calculated on the
basis of the Hamiltonian of Eiles and Greene [29]. This Hamiltonian includes
$s$\- and $p$-wave scattering as well as all relevant interactions of
electronic and nuclear spins. The electron-Cs scattering parameters are
adopted from Ref. [30], and the atomic basis set includes the manifolds
$n=31-33$ and all low-$\ell$ states within this range. The inclusion of the
nuclear spin of the ground-state atom leads to the appearance of two classes
of molecular states: ${}^{3}\Sigma^{+}$ states with a binding resulting only
from scattering in the triplet states of Cs-, and ${}^{1,3}\Sigma^{+}$ states
with contributions from both triplet and singlet scattering [31, 19]. These
labels are only approximate because, e.g., the spin-orbit interaction mixes
$\Sigma$ and $\Pi$ character. The remaining good quantum number in the non-
rotating molecule is $\Omega$, the projection of the total angular momentum
(excluding rotation) on the internuclear axis.

Two strong resonances marked by asterisks are assigned to the lowest
vibrational state in the outermost well of the ${}^{3}\Sigma^{+}$ and
${}^{1,3}\Sigma^{+}$ states, which are almost degenerate in $\Omega$. The
resonance positions agree very well with our previous analysis of these
molecular resonances over a large range of $n$ values considering only
$s$-wave scattering [19]. The calculation reveals that the degeneracy of the
potential-energy curves in $\Omega$ is lifted at shorter internuclear
separations, where contributions from $p$-wave scattering become important.
The depth of the inner wells of the ${}^{1,3}\Sigma^{+}$ states centered
around 1480 $a_{0}$ becomes strongly $\Omega$-dependent. These inner wells
support only a single bound state each, whose zero-point energy is determined
from a harmonic approximation of the potential. The bound states in the
$\Omega$-split potentials are expected to yield a regular series of
resonances. Indeed, a regular series of molecular resonances is observable in
the experimental photoassociation spectra close to 35 2P3/2 (Fig. 2, left
panel). The calculated series of binding energies agrees qualitatively with
the observed structure of photoassociation resonances. The tentative
assignment shown in Fig. 2 reduces the overall deviation, but a more detailed
theoretical study is necessary to explain the remaining deviations, as well as
the mechanism for the observed autoionization limiting the molecular lifetimes
to few microseconds. Very similar electronic structures have been predicted
below the 322P3/2 [30] and the 332P1/2 [29] asymptote in caesium, which we
have confirmed in our calculations. The choice of $n$ determines the range of
binding energies and equilibrium distances of the addressed ion-pair states,
which scale as $1/n^{2}$ and $n^{2}$, respectively. Simple estimates based on
the normalisation of the Rydberg wavefunction, the $n^{2}$ dependence of the
classical Rydberg orbit, and the pair-distance distribution in a homogeneous
gas,k show that photoassociation rates should increase linearly with $n$ and
the rates for stimulated emission into the ion-pair state, discussed in the
next section, should scale as $n^{-3}$. We chose to continue our exploration
using molecules bound close to the 35 2P3/2 asymptote as an example from which
rough estimates for other values of $n$ can be obtained.

##  3 Transfer to an ion-pair state

Encouraged by the qualitative agreement between the theoretical model and the
experimental observations presented above, we use the calculated molecular
structure to evaluate the feasibility of transferring the long-range Rydberg
molecules into an ion-pair state. An analysis of the calculated electronic
wavefunctions around $R\sim 1480\,a_{0}$ shows that these states have strong
contributions both from singlet- and triplet-scattering channels, and that the
interaction is indeed dominated by $p$-wave scattering. An investigation of
the spatial component of the electronic wavefunctions at the molecular
equilibrium separation, depicted for $\Omega=5/2$ in Fig. 3, reveals that
these bound states are localized close to the position of maximal gradient of
the wavefunction of the Rydberg electron, i.e., at a node of its probability
density, as expected for a $p$-wave dominated molecular state, see Eq. (1).
While this property facilitates the optical dipole transitions from the
molecular long-range Rydberg state to ion-pair states, it also results in
small transition dipole moments.

Figure 3:  Cut through the electronic-density distribution of the
${}^{1,3}\Sigma^{+}$ state with $\Omega=5/2$ (solid blue line in the right
panel of Fig. 2) at $R=1480\,a_{0}$ in a plane containing the internuclear
axis (chosen along the $z$-axis). A red dot marks the position of the
perturbing ground-state atom. The inset shows the region around the perturber
at an enlarged scale.

In order to estimate the rates for stimulated emission, we focus on these
$(\Omega)\,^{1,3}\Sigma^{+}\,(v=0)$ states (vibrational wavefunction for
$\Omega=5/2$ shown in Fig. 4). The $R$-dependent transition dipole moment
between the long-range Rydberg state and the ion-pair state is obtained by
numerical evaluation of the matrix element of the dipole operator between the
wavefunction of the Rydberg electron (see, e.g., Fig. 3) and the
$(6s)$-single-electron wavefunction of Cs- obtained from a Hartree-Fock-based
calculation [32] in cylindrical coordinates centered at the position of the
ground-state atom under conservation of nuclear and electronic spins. Note
that in the employed electric dipole approximation, only the singlet component
of the long-range Rydberg molecules wavefunction couples to the 1S0 ground
state of Cs-. This neglects the mixing of 3S1 and 1P1 scattering channels by
the spin-orbit interaction in Cs- ([33]) and the calculation is thus likely to
underestimate the transition rate. Only parallel transitions ($\Delta
m_{l}=0$) yield finite transition-dipole moments, perpendicular transition-
dipole moments ($|\Delta m_{l}|=1$) are strongly suppressed. This is not
surprising, because in the ground state, $m_{l}$ is zero and the $|m_{l}|=1$
component of the Rydberg-electron’s wavefunction has a node on the
internuclear axis which was chosen as the $z$-axis. The dependence of the
transition dipole moment on $\Omega$ is found to be weak, with $\mu_{z}$
decreasing by about a factor of two from $\Omega=1/2$ to $\Omega=7/2$. We thus
chose the intermediate $\Omega=5/2$ as exemplary state to continue our
analysis.

Vibrational wavefunctions of initial and final molecular states as well as the
expectation value of the transition-dipole moment are derived by numerical
integration using Numerov’s algorithm. The predicted spontaneous decay rate
(Einstein $A$ coefficient $\ll{1\,\textrm{s}^{-1}}$) is much smaller than the
rate for radiative decay of the long-range Rydberg molecules, which can be
well approximated by the decay rate for an isolated Rydberg atom (e.g., for Cs
352P3/2, Ref. [34] calculated a total decay rate, including spontaneous decay
and transitions induced by black-body radiation, of 28 kHz). It is therefore
necessary to enhance the transition rate, e.g., by stimulated emission in an
external light field. The required wavelength lies in the mid-infrared region
of the optical spectrum (around 2.6 $\mu$m) where intense, broadly tunable cw
laser sources are rare. We thus evaluate the transition rates for pulsed laser
sources, such as a pulse-amplified cw laser system combined with difference-
frequency mixing [35], and calculate the Rabi frequency for the driven
transition. For realistic laser parameters given in the caption of Fig. 4, the
product of Rabi frequency and pulse length reaches unity, indicating that the
transition can be saturated.

Within the Franck-Condon window, final heavy-Rydberg states between
$n_{\mathrm{HR}}\sim 9\,350$ and $n_{\mathrm{HR}}\sim 9\,600$ can be
addressed, corresponding to an energy span of about 240 GHz. The expected
vibrational spacing of the addressed heavy-Rydberg states is around 1 Ghz,
which is well resolvable by a mid-IR laser pulse of a few nanoseconds length
(transform-limited bandwidth typically below 150 MHz). Transitions will be
observable by a loss in the yield of field-ionized Rydberg molecules and the
appearance of Cs- anions on resonance. The addressed heavy-Rydberg states of
Cs with $n_{\mathrm{HR}}\sim 10\,000$ have a significantly increased
sensitivity to electric fields: classically, the cation-anion system at an
internuclear separation of $R\sim 1\,500a_{0}$ constitutes an electric dipole
with a dipole moment $\mu=1\,500ea_{0}$, where $e$ is the elementary charge.
Assuming random orientations of this dipole in an external electric field of
10 mV/cm, the resulting transition broadening amounts to $\sim$40 MHz. The
Inglis-Teller limit, i.e., the electric-field strength where neighbouring
Rydberg manifolds start to overlap, is reached at a field strength of
$\sim$250 mV/cm [14]. Given typical residual fields in current experiments of
$\leq 10$ mV/cm [36] and the envisioned laser linewidth of $\sim$150 MHz, the
electric-field induced broadening should not hinder the spectroscopic study of
heavy-Rydberg states.

Figure 4:  Left panel: (red line, left axis) calculated $R$-dependent
transition dipole moment ($\Delta m_{l}=0$) between the electronic state
${}^{1,3}\Sigma^{+}$ ($\Omega=5/2$), whose spatial wavefunction is presented
in Fig. 3, and the 1S0 ground state of Cs-; (blue line, right axis)
vibrational wavefunction of the long-range Rydberg molecule in $v=0$; (gray
line, right axis) vibrational wavefunction of a Cs2 heavy-Rydberg state in
$n_{\mathrm{HR}}=9\,500$. Right panel: Calculated product of Rabi rate and
pulse width ($\tau_{\mathrm{laser}}=5$ ns) as a function of the vibrational
level of the final heavy-Rydberg state $n_{\mathrm{HR}}$ for a pulse energy of
10 $\mu$J and a beam waist radius of 300 $\mu$m.

Direct photoionization of the Rydberg molecule by the mid-IR radiation
competes with stimulated emission. The photoionization rates for a long-range
Rydberg molecule are very similar to the rates for isolated Rydberg atoms,
which have been calculated by Ovsiannikov et al. [37]. Based on their
calculations, the ionization probability for the laser parameters given in the
caption of Fig. 4 (a Gaussian pulse with temporal full width at half maximum
of 5 ns and a peak intensity of about $2\cdot 10^{6}$ W/cm2) is about 3%,
significantly smaller than the probability for stimulated decay (which is
about unity for the given parameters).

The lifetime of heavy-Rydberg states is in general limited by couplings to
valence states of the neutral parent species at short internuclear distances,
which can cause dissociation or autoionization of the ion-pair state as
observed in other heavy-Rydberg systems [14, 15]. Estimating lifetimes of the
excited ion-pair states goes beyond the scope of this paper. However, the
lifetimes should increase with increasing rotational angular momentum of the
ion pair.

##  4 Discussion and outlook

In this article we have presented a novel route to excite heavy-Rydberg states
with high values of $n_{\mathrm{HR}}$ via photoassociation of long-range
Rydberg molecules and stimulated emission into the ion-pair state. By
performing detailed calculations for a specific system (molecular states close
to the 6 2S1/2-35 2P3/2 asymptote in caesium), we find this approach to be
feasible for realistic experimental parameters. We note that there is nothing
specific about the chosen system, except for the criteria discussed above, and
that we are confident that the same approach can be implemented in a variety
of systems, employing photoassociation close to, e.g., $n$ 2S1/2, $n$ 2Pj or
$n$ 2Dj asymptotes in homonuclear (Cs2, Rb2, K2, …) or heteronuclear alkali
dimers [38], as long as the formed long-range Rydberg states have significant
contributions from $p$-wave scattering. The scheme can be similarly applied to
more complex chemical species, such as, e.g., dimers composed of the alkaline
earth metals with stable anions (Ca, Sr, Ba). For example, the formation of
long-range Sr2 Rydberg molecules has already been demonstrated [39]. Because
the ground state of alkaline earth anions is a 2Pj state, efficient transfer
from $s$-wave-scattering dominated long-range Rydberg molecules to ion-pair
states can be expected in these systems. The lower electron affinity of
alkaline-earth atoms will, however, require to stimulate the transitions by
radiation in the terahertz region.

The approach presented in this work opens up the ways to study structure and
dynamics of heavy-Rydberg states in a previously unexplored region of
internuclear separations. As discussed in the introduction, a dense gas of ion
pairs may transform into a plasma with equal-mass charges (for homonuclear
systems) or charges with variable mass ratios (for heteronuclear systems).
This transformation may occur spontaneously, as routinely observed in dense
gases of Rydberg atoms [40, 41, 42, 43], or be stimulated by DC-electric-
field-induced dissociation [15] or multi-photon dissociations by low-frequency
radio-frequency (RF) fields [44]. A theoretical study by Robicheaux et. al
[13] predicts the equilibration dynamics of plasmas with equal-mass charges to
be on experimentally accessible timescales in the range of microseconds. The
non-equilibrium dynamics of such a strongly coupled system are expected to
exhibit strong spatial correlations between particles and show physical
phenomena relevant to the description of giant gas planets like Jupiter or
fusion devices [11]. The neutral plasma might even be confined in the RF field
of a hybrid atom-ion trap [45], because the RF parameters can be easily chosen
such that the masses of all components fall within the trap’s stability window
[46].

We thank Frédéric Merkt for his continuous and generous support. We thank
Matthew Eiles, Chris Greene, Francis Robicheaux, Christian Fey, Frederic
Hummel, Herwig Ott, and Hossein Sadeghpour for fruitful discussions. This work
is supported financially by the ETH Research Grant ETH-22 15-1 and the NCCR
QSIT of the Swiss National Science Foundation.

## References

  * [1] J. P. Shaffer, S. T. Rittenhouse, and H. R. Sadeghpour.  Ultracold Rydberg molecules.  Nature Communications, 9(1):1965, 2018.
  * [2] M. T. Eiles.  Trilobites, butterflies, and other exotic specimens of long-range Rydberg molecules.  Journal of Physics B: Atomic, Molecular and Optical Physics, 52(11):113001, 2019.
  * [3] E. Fermi.  Sopra lo Spostamento per Pressione delle Righe Elevate delle Serie Spettrali.  Il Nuovo Cimento, 11(3):157–166, 1934.
  * [4] C. Füchtbauer, G. Joos, and O. Dinkelacker.  Über Intensität, Verbreiterung und Druckverschiebung von Spektrallinien, insbesondere der Absorptionslinie 2537 des Quecksilbers.  Annalen der Physik, 376(9-12):204–227, 1923.
  * [5] E. Amaldi and E. Segrè.  Effetto della Pressione Sui Termini Elevati Degli Alcalini.  Il Nuovo Cimento, 11(3):145–156, 1934.
  * [6] A. Omont.  On the theory of collisions of atoms in Rydberg states with neutral particles.  Journal de Physique, 38(11):1343–1359, 1977.
  * [7] C. H. Greene, A. S. Dickinson, and H. R. Sadeghpour.  Creation of Polar and Nonpolar Ultra-Long-Range Rydberg Molecules.  Physical Review Letters, 85(12):2458–2461, 2000.
  * [8] G. Herzberg and Ch. Jungen.  Rydberg series and ionization potential of the H2 molecule.  Journal of Molecular Spectroscopy, 41(3):425–486, 1972.
  * [9] M. Beyer and F. Merkt.  Communication: Heavy-Rydberg states of HD and the electron affinity of the deuterium atom.  The Journal of Chemical Physics, 149(3):031102, 2018.
  * [10] G. Cerchiari, A. Kellerbauer, M. S. Safronova, U. I. Safronova, and P. Yzombard.  Ultracold Anions for High-Precision Antihydrogen Experiments.  Physical Review Letters, 120(13):133205, 2018.
  * [11] T. C. Killian, T. Pattard, T. Pohl, and J. M. Rost.  Ultracold neutral plasmas.  Physics Reports, 449(4–5):77–130, 2007.
  * [12] M. Lyon and S. L. Rolston.  Ultracold neutral plasmas.  Reports on Progress in Physics, 80(1):017001, 2016.
  * [13] F. Robicheaux, B. J. Bender, and M. A. Phillips.  Simulations of an ultracold, neutral plasma with equal mass for every charge.  Journal of Physics B: Atomic, Molecular and Optical Physics, 47(24):245701, 2014.
  * [14] E. Reinhold and W. Ubachs.  Heavy Rydberg states.  Molecular Physics, 103(10):1329–1352, 2005.
  * [15] S. Mollet and F. Merkt.  Dissociation dynamics of ion-pair states of Cl2 at principal quantum numbers beyond 1500.  Physical Review A, 82(3):032510, 2010.
  * [16] A. Kirrander, S. Rittenhouse, M. Ascoli, E. E. Eyler, P. L. Gould, and H. R. Sadeghpour.  Approach to form long-range ion-pair molecules in an ultracold Rb gas.  Physical Review A, 87(3):031402, 2013.
  * [17] S. Markson and H. R. Sadeghpour.  A model for charge transfer in ultracold Rydberg ground-state atomic collisions.  Journal of Physics B: Atomic, Molecular and Optical Physics, 49(11):114006, 2016.
  * [18] A. Kramida, Yu. Ralchenko, J. Reader, and and NIST ASD Team.  NIST Atomic Spectra Database (ver. 5.7), [Online]. Available: https://physics.nist.gov/asd [2017, April 9]. National Institute of Standards and Technology, Gaithersburg, MD., 2019.
  * [19] H. Saßmannshausen, F. Merkt, and J. Deiglmayr.  Experimental Characterization of Singlet Scattering Channels in Long-Range Rydberg Molecules.  Phys. Rev. Lett., 114(13):133201, 2015.
  * [20] H. Saßmannshausen, J. Deiglmayr, and F. Merkt.  Long-range Rydberg molecules, Rydberg macrodimers and Rydberg aggregates in an ultracold Cs gas.  Eur. Phys. J. Spec. Top., 225(15-16):2891–2918, 2016.
  * [21] M. Scheer, J. Thøgersen, R. C. Bilodeau, C. A. Brodie, H. K. Haugen, H. H. Andersen, P. Kristensen, and T. Andersen.  Experimental Evidence that the 6s6p 3PJ states of Cs- are shape resonances.  Physical Review Letters, 80(4):684–687, 1998.
  * [22] E. L. Hamilton, C. H. Greene, and H. R. Sadeghpour.  Shape-resonance-induced long-range molecular Rydberg states.  Journal of Physics B: Atomic, Molecular and Optical Physics, 35(10):L199–L206, 2002.
  * [23] T. Niederprüm, O. Thomas, T. Eichert, C. Lippe, J. Pérez-Ríos, C. H. Greene, and H. Ott.  Observation of pendular butterfly Rydberg molecules.  Nature Communications, 7:12820, 2016.
  * [24] T. F. Gallagher.  Rydberg Atoms.  Cambridge University Press, Cambridge, 1994.
  * [25] F. Engel, T. Dieterle, F. Hummel, C. Fey, P. Schmelcher, R. Löw, T. Pfau, and F. Meinert.  Precision Spectroscopy of Negative-Ion Resonances in Ultralong-Range Rydberg Molecules.  Physical Review Letters, 123(7):073003, 2019.
  * [26] J. L. MacLennan, Y. Chen, and G. Raithel.  Deeply bound ($24{{D}}_{{J}}+5{{S}}_{1/2}$) ${}^{87}\mathrm{Rb}$ and ${}^{85}\mathrm{Rb}$ molecules for eight spin couplings.  Physical Review A, 99(3):033407, 2019.
  * [27] J. Deiglmayr, H. Herburger, H. Saßmannshausen, P. Jansen, H. Schmutz, and F. Merkt.  Precision measurement of the ionization energy of Cs I.  Phys. Rev. A, 93(1):013424, 2016.
  * [28] T. Niederprüm, O. Thomas, T. Manthey, T. M. Weber, and H. Ott.  Giant Cross Section for Molecular Ion Formation in Ultracold Rydberg Gases.  Physical Review Letters, 115(1):013003, 2015.
  * [29] M. T. Eiles and C. H. Greene.  Hamiltonian for the inclusion of spin effects in long-range Rydberg molecules.  Physical Review A, 95(4):042515, 2017.
  * [30] S. Markson, S. T. Rittenhouse, R. Schmidt, J. P. Shaffer, and H. R. Sadeghpour.  Theory of Ultralong-Range Rydberg Molecule Formation Incorporating Spin-Dependent Relativistic Effects: Cs(6s)–Cs(np) as Case Study.  ChemPhysChem, 17(22):3683–3691, 2016.
  * [31] D. A. Anderson, S. A. Miller, and G. Raithel.  Angular-momentum couplings in long-range Rb2 Rydberg molecules.  Physical Review A, 90(6):062518, 2014.
  * [32] R. D. Cowan.  The Theory of Atomic Structure and Spectra.  Number 3 in Los Alamos Series in Basic and Applied Sciences. University of California Press, Berkeley - California, 1981.
  * [33] A. A. Khuskivadze, M. I. Chibisov, and I. I. Fabrikant.  Adiabatic energy levels and electric dipole moments of Rydberg states of Rb2 and Cs2 dimers.  Physical Review A, 66(4):042709, 2002.
  * [34] V. D. Ovsiannikov, I. L. Glukhov, and E. A. Nekipelov.  Rates of blackbody radiation-induced transitions from Rydberg states of alkali atoms.  Journal of Physics B: Atomic, Molecular and Optical Physics, 44(19):195010, 2011.
  * [35] U. Jacovella, J. A. Agner, H. Schmutz, J. Deiglmayr, and F. Merkt.  Infrared spectroscopy of molecular ions in selected rotational and spin-orbit states.  The Journal of Chemical Physics, 145(1):014301, 2016.
  * [36] M. Peper, F. Helmrich, J. Butscher, J. A. Agner, H. Schmutz, F. Merkt, and J. Deiglmayr.  Precision measurement of the ionization energy and quantum defects of 39K I.  Physical Review A, 100(1):012501, 2019.
  * [37] V. D. Ovsiannikov, I. L. Glukhov, and E. A. Nekipelov.  Radiative lifetime and photoionization cross-section for Rydberg states in alkali-metal atoms.  Optics and Spectroscopy, 111(1):25, 2011.
  * [38] M. T. Eiles.  Formation of long-range Rydberg molecules in two-component ultracold gases.  Physical Review A, 98(4):042706, 2018.
  * [39] B. J. DeSalvo, J. A. Aman, F. B. Dunning, T. C. Killian, H. R. Sadeghpour, S. Yoshida, and J. Burgdörfer.  Ultra-long-range Rydberg molecules in a divalent atomic system.  Physical Review A, 92(3):031403, 2015.
  * [40] P. J. Tanner, J. Han, E. S. Shuman, and T. F. Gallagher.  Many-Body Ionization in a Frozen Rydberg Gas.  Physical Review Letters, 100(4):043002, 2008.
  * [41] M. Robert-de-Saint-Vincent, C. S. Hofmann, H. Schempp, G. Günter, S. Whitlock, and M. Weidemüller.  Spontaneous Avalanche Ionization of a Strongly Blockaded Rydberg Gas.  Physical Review Letters, 110(4):045004, 2013.
  * [42] H. Saßmannshausen, F. Merkt, and J. Deiglmayr.  Pulsed excitation of Rydberg-atom-pair states in an ultracold Cs gas.  Phys. Rev. A, 92(3):032505, 2015.
  * [43] G. T. Forest, Y. Li, E. D. Ward, A. L. Goodsell, and D. A. Tate.  Expansion of an ultracold Rydberg plasma.  Physical Review A, 97(4):043401, 2018.
  * [44] A. Arakelyan, J. Nunkaew, and T. F. Gallagher.  Ionization of Na Rydberg atoms by a 79-GHz microwave field.  Physical Review A, 94(5):053416, 2016.
  * [45] J. Deiglmayr, A. Göritz, T. Best, M. Weidemüller, and R. Wester.  Reactive collisions of trapped anions with ultracold atoms.  Physical Review A, 86(4): 043438, 2012.
  * [46] R. Wester.  Radiofrequency multipole traps: Tools for spectroscopy and dynamics of cold molecular ions.  Journal of Physics B: Atomic, Molecular and Optical Physics, 42(15):154001, 2009.


# Adjustment with Three Continuous Variables

Brian Knaeble

(Date: June 2017)

###### Abstract.

Spurious association between $X$ and $Y$ may be due to a confounding variable
$W$. Statisticians may adjust for $W$ using a variety of techniques. This
paper presents the results of simulations conducted to assess the performance
of those techniques under various, elementary, data-generating processes. The
results indicate that no technique is best overall and that specific
techniques should be selected based on the particulars of the data-generating
process. Here we show how causal graphs can guide the selection or design of
techniques for statistical adjustment. R programs are provided for researchers
interested in generalization.

##  1\. Introduction

Statistical adjustment has meant different things to different people over the
years. Demming (1964) wrote a book titled Statistical Adjustment of Data,
focusing on the use of known facts to lesson errors in measurement.
Epidemiologists routinely make adjustments (Kasim, 2012), and cancer
researchers report incidence and mortality rates after adjusting for variables
such as age and ethnicity (Howlader et al., 2016). Whether or not to adjust is
a question relevant even to those studying workforce discrimination (Schieder
and Gould, 2016).

The slope coefficient from simple (bivariate) linear regression may be
adjusted to account for a third lurking variable (Lu, 2009). Multiple
regression is one of many ways to make an adjustment. Stratification over the
third variable is another possibility, and there are many varieties of
matching procedures, should the original regressor be dichotomous. When all
three variables are continuous adjustment can be accomplished using partial
and semi-partial correlation.

In some instances conditioning on a third variable can increase bias (Pearl,
2009b). Given an accurate causal graph (also known as a Directed Acyclic Graph
or DAG) the back-door criterion can be used to determine an admissible set of
covariates on which to condition for unbiased estimation of a causal effect
(Pearl, 2009a). Some argue for conditioning on as many pre treatment variables
as possible (Rubin, 2009), but this is not without controversy. See the
introduction of Ding and Miratrix (2015) for an overview.

Here we have conducted numerous simulations, restricting our attention to the
simplest case of only three variables. We assume linearity and continuity, and
we study 33 different data generating processes to assess the performance of
multiple regression compared with simple regression and methods relating to
partial correlation. For detailed description of assumptions and studied
adjustment techniques see Section 2. For graphical displays of results see
Section 3. Commentary on the results is provided in Section 4. For further
discussion and conclusions see Section 5.

##  2\. Methods

There are 25 non-cyclic graphs of three variables upon which to generate data.
The classic confounding case is shown here:
$X$$Y$$W$ Each of the three arrows can be independently modified so as to
vanish or reverse direction. For example, $W$ could be an intermediate
variable as shown here:
$X$$Y$$W$ Most would agree given the latter causal graph that conditioning on
$W$ is inappropriate, because a change in $X$ can not produce a change in $Y$
when $W$ is held fixed.

We write $X\rightarrow Y$ to represent a direct causal effect of $X$ on $Y$.
We use a dash to represent the absence of a direct causal effect, e.g. $X-Y$.
We write $X\leftarrow Y$ when $Y$ causes $X$. In this latter case the true
causal effect of $X$ on $Y$ is zero. Analogous notation and interpretation
holds for relationships with $W$. $X$ and $Y$ are defined or redefined if
necessary so that any causal effect between them is positive not negative.
When there are two causal effects involving $W$, we default to them both being
positive, and we allow for the effect between $W$ and $X$ to be positive or
negative. We signal a negative effect with a lower case n, e.g.
$X\leftarrow^{n}W$. We use the adjective “twisted” when the direction of the
effect between $W$ and $X$ differs from the direction of the effect between
$W$ and $Y$. We do not allow cycles.

For consistency, data is generated assuming normality and unit variance, and
associated with each arrow is a true slope coefficient of magnitude
$1/\sqrt{3}$. To maintain unit variance for a response variable, defined as
any variable at the tip of an arrow, the error term has standard deviation
$\sqrt{2}/\sqrt{3}$ when there is one cause, and $1/\sqrt{3}$ for each when
there are two causes. The sample size can be set arbitrarily, and here we have
chosen to set the sample size at $n=30$.

The following R commands,

    x=rnorm(30,0,1)
    y=(1/sqrt(3))*x+rnorm(30,0,sqrt(2)/sqrt(3))
    w=(1/sqrt(3))*x+(1/sqrt(3))*y+rnorm(30,0,1/sqrt(3))

generate data for this graph:
$X$$Y$$W$ In general, following our assumptions, each graph has its own data
generating process.

We assume that researchers, blind to the data generating process, will attempt
to estimate the causal effect of $X$ on $Y$ by 1. regressing $Y$ on $X$
(Simple Regression), 2. regressing $Y$ on $X$ and $W$ (Multiple Regression),
3. regressing $Y$ on the residuals from a regression of $X$ on $W$ (Residual
$X$), 4. regressing the residuals from a regression of $Y$ on $W$, onto $X$
(Residual $Y$), 5. regressing the residuals from a regression of $Y$ on $W$,
onto the residuals from a regression of $X$ on $W$ (Residual $X$ and $Y$), or
6. regressing $Y$ onto the fitted values of a regression of $X$ onto $W$
(Fitted $X$). We show in the appendix that Multiple Regression, Residual $X$,
and Residual $X$ and $Y$ are equivalent, each having conditioned on $W$. The
results thus display estimates from only Simple Regression, Multiple
Regression, Residual $Y$, and Fitted $X$.

##  3\. Results

For each data generating process, one thousand datasets were produced, and the
four adjustment techniques applied to each. The true, total causal effect is
defined from the data generating process, using the chain rule to multiply
slope coefficients of a unidirectional causal chain, and adding effects across
paths when appropriate. The true, total causal effect is represented
graphically with a vertical black line. Figures are classified by how $W$
relates to both $X$ and $Y$ within the data generating process, while within-
figure-plots are classified by how $X$ relates to $Y$ within the data
generating process.

Figure 1. Confounding: $X\leftarrow W\rightarrow Y$. Figure 2. Twisted
Confounding: $X\leftarrow^{n}W\rightarrow Y$. Figure 3. $W$ as a collider:
$X\rightarrow W\leftarrow Y$. Figure 4. $W$ as a twisted collider:
$X\rightarrow^{n}W\leftarrow Y$. Figure 5. Instrumental W: $X\leftarrow W-Y$.
Figure 6. Weak Confounding: $X-W\rightarrow Y$. Figure 7. $W$ Post Treatment:
$X\rightarrow W-Y$. Figure 8. $W$ Post Response: $X-W\leftarrow Y$. Figure 9.
$W$ Irrelevant: $X-W-Y$. Figure 10. Forward Chain: $X\rightarrow W\rightarrow
Y$; cycle excluded, replaced by twisted forward chain. Figure 11. Backward
Chain: $X\leftarrow W\leftarrow Y$; cycle excluded, replaced by twisted
backward chain.

##  4\. Comments

With an unobserved confounder $U$ as shown here
$X$$Y$$W$$U$ we would expect Fitted $X$ to outperform Simple Regression
(Morgan and Winship, 2007, Chapter 7). Fitted $X$ matches the causal structure
of Figure 5, where it is the only consistently unbiased technique. Fitted $X$
performs poorly otherwise.

Residual $Y$ is surprisingly robust against a mis specified direction for the
causal effect between $X$ and $Y$. Take Figure 1 as an example. There we see
for each plot that the bias of Residual $Y$ is bounded by approximately $0.3$.
This is a better bound than any of the alternatives. When estimating the
causal effect of $X$ on $Y$ a spurious result is possible when $Y$ causes $X$
in reality. The risk for such a spurious result can be reduced in some
situations by using Residual $Y$.

Another observed benefit of Residual $Y$ is its generally small variance. With
a large sample and a known causal graph, theory says (Pearl, 2009b;
VanderWeele and Shpitser, 2011) when $W\rightarrow X$ then Multiple Regression
should be used in place of Simple Regression and when $X\rightarrow W$ then
Simple Regression should be used instead of Multiple Regression, and this is
consistent with our results. Our results show the utility of additional
adjustment techniques when samples are smaller or when there is uncertainty
regarding the structure of the causal graph.

##  5\. Discussion

We have displayed results for sample size $n=30$ using box plots so as to show
bias and sampling error. Those interested in large sample asymptotic theory
can focus on the medians of the box plots, while those interested in small
sample application can focus on the width or span of the box plots. We have
assumed normality and linearity throughout. In some cases real data may meet
these assumptions, and in other cases the data may be transformed so as to
meet these assumptions, but often data will not satisfy these assumptions.
Also, the data generating process may involve more than three variables, some
or all possibly categorical. The R code used for our simulations should be
adapted and generalized before application to a specific analysis. Our simple
results are meant to guide intuition and spur new thought.

It is hard to believe that any proposed causal graph perfectly matches
reality. Given a set of candidate graphs, we can conduct simulations to assess
the performance of various adjustment techniques across data generating
processes associated with the candidate graphs, before selecting an adjustment
technique in accordance with costs associated with our specific problem. If
costs rise dramatically only for very large errors then the conservative
choice of Residual Y adjustment may be preferable to Multiple Regression
adjustment, especially when reverse causality is possible. On the other hand,
when there is a high degree of confidence in a causal graph then the
simulations may reveal a single best adjustment techniques that matches
theoretical recommendations for very precise estimation.

Pearl (2009a, p. 123) states that “adjustment amounts to partitioning the
population into groups that are homogeneous with respect to (the third
variable)”, while Ding and Miratrix (2015) titled their paper “To Adjust or
Not to Adjust?”. Much concern (VanderWeele and Shpitser, 2011) has been given
to the problem of determining which covariates to condition on. Rubin (2009)
argues against leaving observed covariates unbalanced. It has been suggested
that ignoring covariate data is contrary to Bayesianism (Gellman, 2009).
Armistead (2014) shows how to use such covariate information. Here we have
simply presented the results of basic simulations in support of the notion
that we should ask not only whether we should adjust or not, but also how we
should best adjust.

##  6\. Appendix

######  Lemma 6.1.

The estimates for “Residual $X$” and “Residual $X$ and $Y$” are both identical
to the estimate for “Multiple Regression”.

###### Proof.

Let $\bf{x}$,$\bf{w}$, and $\bf{y}$ be observed vectors of data associated
with $n$ observations of $(X,W,Y)$. We assume mean zero vectors of data
without loss of generality so $\bf{x}$,$\bf{w}$, and $\bf{y}$ are within a
linear, $(n-1)$-dimensional, subspace of $\mathbb{R}^{n}$ orthogonal to a
vector of ones. All that follows takes place within this subspace. We assume
also that $\mathbf{w}\neq k_{1}\mathbf{x}$ for any $k_{1}$.

Let $\mathbf{p}_{1}$ be the projection of $\mathbf{y}$ onto the span of
$\mathbf{x}$ and $\mathbf{w}$. Let $\mathbf{p}_{2}$ be the projection of
$\mathbf{x}$ onto $\mathbf{w}$, and note that $\mathbf{p}_{2}=k_{2}\mathbf{w}$
for some $k_{2}$. Let $\mathbf{p}_{3}$ be the projection of $\mathbf{y}$ onto
$\mathbf{w}$, and note that $\mathbf{p}_{3}=k_{3}\mathbf{w}$ for some $k_{3}$.

Let $\alpha_{1}$ be the Multiple Regression slope estimate associated with
$X$. By least-squares we have
$\alpha_{1}\mathbf{x}+c_{1}\mathbf{w}=\mathbf{p}_{1}$ for some $c_{1}$. Let
$\mathbf{p}_{4}$ be the projection of $\mathbf{y}$ onto the span of
$\mathbf{x}-\mathbf{p}_{2}$ and $\mathbf{w}$. Note that
$\mathbf{p}_{4}=\mathbf{p}_{1}$.

Let $\alpha_{2}$ be the Residual $X$ slope estimate. We have
$\alpha_{2}(\mathbf{x}-\mathbf{p}_{2}))+c_{2}\mathbf{w}=\mathbf{p}_{4}$, which
can be solved with $\alpha_{2}=\alpha_{1}$ and $c_{2}=c_{1}+\alpha_{1}k_{2}$.

Let $\alpha_{3}$ be the Residual $X$ and $Y$ slope estimate. Note that
$\mathbf{p}_{1}-\mathbf{p}_{3}$ is the projection of
$\mathbf{y}-\mathbf{p}_{3}$ onto the span of $\mathbf{x}-\mathbf{p}_{2}$ and
$\mathbf{w}$. We have
$\alpha_{3}(\mathbf{x}-\mathbf{p}_{2}))+c_{3}\mathbf{w}=\mathbf{p}_{1}-\mathbf{p}_{3}$,
which (since $\mathbf{p}_{3}=k_{3}\mathbf{w}$) can be solved with
$\alpha_{3}=\alpha_{2}$ and $c_{3}=c_{2}+k_{3}$.

∎

    ### R Program for simulations
    k=1000 ### number of samples
    v1=numeric(k)
    v2=numeric(k)
    v3=numeric(k)
    v4=numeric(k)
    for (i in 1:k) {
    ### insert particular data generating process below
    x=rnorm(30,0,1)
    w=(1/sqrt(3))*x+rnorm(30,0,sqrt(2)/sqrt(3))
    y=(1/sqrt(3))*w+(1/sqrt(3))*x+rnorm(30,0,1/sqrt(3))
    ### the above example is for a forward chain
    ### the true value is (1/sqrt(3))*(1/sqrt(3))+(1/sqrt(3))
    f=lm(x~w)$fitted.values
    g=lm(y~w)$residuals
    v1[i]=summary(lm(y~x))$coefficients[2,1]
    v2[i]=summary(lm(y~x+w))$coefficients[2,1]
    v3[i]=summary(lm(g~x))$coefficients[2,1]
    v4[i]=summary(lm(y~f))$coefficients[2,1]
    }
    ### R program for graphing box plots
    par(las=1)
    par(mar=c(5, 9, 4, 2) + 0.1)
    boxplot(v1,v2,v3,v4,
    xlab=c("Estimate for causal effect of X on Y"),
    horizontal=T,
    names=c(v1="Simple Regression",v2="Multiple Regression",
    v3="Residual Y",v4="Fitted X"),
    boxwex=.3)
    abline(v=(1/sqrt(3))*(1/sqrt(3))+(1/sqrt(3)))

## References

  * Armistead (2014) Armistead, T. 2014. Resurrecting the Third Variable: A Critique of Pearl’s Causal Analysis of Simpson’s Paradox. The American Statistician. Volume 68, Issue 1.
  * Carroll and Huff (1980) Carroll, K. and Huff M. 1980. Nutrition and food science, present knowledge and utilization Vol. 3. Plenum Press.
  * Demming (1964) Deming, E. 1964. Statistical Adjustment of Data. Dover Publications (November 17, 2011).
  * Ding and Miratrix (2015) Ding, P. and Miratrix, L. 2015. To Adjust or Not to Adjust? Sensitivity Analysis of M-Bias and Butterfly-Bias. Journal of Causal Inference.
  * Gellman (2009) Gellman, A. 2009. Resolving disputes between J. Pearl and D. Rubin on causal inference. andrewgelman.com. Posted by Andrew on 5 July 2009, 12:32pm.
  * Howlader et al. (2016) Howlader N, Noone AM, Krapcho M, Miller D, Bishop K, Kosary CL, Yu M, Ruhl J, Tatalovich Z, Mariotto A, Lewis DR, Chen HS, Feuer EJ, Cronin KA (eds). 2016. SEER Cancer Statistics Review, 1975-2014. National Cancer Institute, Bethesda, MD, based on November 2016 SEER data submission, posted to the SEER web site, April 2017.
  * Kasim (2012) Kasim, K. 2012. Basic Concepts of Modern Epidemiology. Lambert Academic Publishing.
  * Krit et al. (1982) Kritchevsky, D., Tepper, S. A., Czarnecki, S. K. and Klurfeld, D. M. 1982. Atherogenicity of animal and vegetable protein. Influence of the lysine to arginine ratio. Atherosclerosis. 41, 429-431.
  * Kritchevsky et al. (1983) Kritchevsky, D., Tepper, S. A., Czarnecki, S. K., Klurfeld, D. M., and Story, J. A. 1983. Current Topics in Nutrition and Disease, Volume 8: Animal and Vegetable Proteins in Lipid Metabolism and Atherosclerosis, 85-100.
  * Lu (2009) Lu, Christine Y. 2009. Observational studies: a review of study designs, challenges and strategies to reduce confounding. International journal of clinical practice, volume 63, number 5, p 691-697.
  * Morgan and Winship (2007) Morgan, S.L. and Winship, C. 2007. Counterfactuals and Causal Inference, Methods and Principles for Social Research. Cambridge.
  * Pearl (2009a) Pearl, Judea. 2009. Causality: models, reasoning and inference. Cambridge University Press.
  * Pearl (2009b) Pearl, Judea. 2009. Causal inference in statistics: An overview. Statistical Surveys, volume 3, p 96-146.
  * Rubin (2009) Rubin, D. 2009. Should observational studies be designed to allow lack of balance in covariate distributions across treatment groups? Statistics in Medicine, 28, 1420-1423.
  * Schieder and Gould (2016) Schieder, Jessica and Gould, Elise. 2016. Women’s work and the gender pay gap. Economic Policy Institute report.
  * Terpstra, Hermus, and West (1983) Terpstra, A. H. M., Hermus, R. J. J. and West, C. E. 1983. Animal and Vegetable Proteins in Lipid Metabolism and Athersclerosis. Alan R. Liss, Inc.
  * VanderWeele and Shpitser (2011) VanderWeele, T.J. and Shpitser, I. 2011. A new criterion for confounder selection. Biometrics, 67, 1406-1413.
  * Wasserman (2004) Wasserman, L. (2004). All of Statistics, A Concise Course in Statistical Inference. Springer.