# Detection of Carbon Monoxide in the Atmosphere of WASP-39b Applying Standard
Cross-Correlation Techniques to JWST NIRSpec G395H Data

Emma Esparza-Borges Instituto de Astrofísica de Canarias, E-38200 La Laguna,
Tenerife, Spain  Departamento de Astrofísica, Universidad de La Laguna,
E-38206 La Laguna, Tenerife, Spain  Mercedes López-Morales Center for
Astrophysics ${\rm\mid}$ Harvard & Smithsonian, 60 Garden St, Cambridge, MA
02138, USA  Jéa I. Adams Redai Center for Astrophysics ${\rm\mid}$ Harvard &
Smithsonian, 60 Garden St, Cambridge, MA 02138, USA  Enric Pallé Instituto de
Astrofísica de Canarias, E-38200 La Laguna, Tenerife, Spain  Departamento de
Astrofísica, Universidad de La Laguna, E-38206 La Laguna, Tenerife, Spain
James Kirk Department of Physics, Imperial College London, Prince Consort
Road, London, SW7 2AZ, UK  Núria Casasayas-Barris Instituto de Astrofísica de
Canarias, E-38200 La Laguna, Tenerife, Spain  Natasha E. Batalha NASA Ames
Research Center, Moffett Field, CA 94035, USA  Benjamin V. Rackham 51 Pegasi b
Fellow  Department of Earth, Atmospheric and Planetary Sciences, Massachusetts
Institute of Technology, 77 Massachusetts Avenue, Cambridge, MA 02139, USA
Kavli Institute for Astrophysics and Space Research, Massachusetts Institute
of Technology, Cambridge, MA 02139, USA  Jacob L. Bean Department of Astronomy
& Astrophysics, University of Chicago, Chicago, IL 60637, USA  S. L. Casewell
Centre for Exoplanet Research, School of Physics and Astronomy, University of
Leicester, University Road, Leicester, LE1 7RH, UK  Leen Decin Institute of
Astronomy, KU Leuven, Celestijnenlaan 200D, 3001 Heverlee, Belgium  Leonardo
A. Dos Santos Space Telescope Science Institute, 3700 San Martin Drive,
Baltimore, MD 21218, USA  Antonio García Muñoz Université Paris-Saclay,
Université Paris Cité, CEA, CNRS, AIM, 91191, Gif-sur-Yvette, France  Joseph
Harrington Planetary Sciences Group, Department of Physics and Florida Space
Institute, University of Central Florida, Orlando, FL, USA  Kevin Heng Ludwig
Maximilian University of Munich, Faculty of Physics, Germany  Renyu Hu Jet
Propulsion Laboratory, California Institute of Technology, Pasadena, CA 91109,
USA  Division of Geological and Planetary Sciences, California Institute of
Technology, Pasadena, CA 91125, USA  Luigi Mancini Department of Physics,
University of Rome “Tor Vergata”, Rome, Italy  Max Planck Institute for
Astronomy, Heidelberg, Germany  INAF – Turin Astrophysical Observatory, Pino
Torinese, Italy
Karan Molaverdikhani Universitäts-Sternwarte, Fakultät für Physik, Ludwig-
Maximilians-Universität München, Scheinerstr. 1, 81679 München, Germany
Exzellenzcluster ‘Origins’, Boltzmannstr. 2, 85748 Garching, Germany
Giuseppe Morello Department of Space, Earth and Environment, Chalmers
University of Technology, SE-412 96 Gothenburg, Sweden;  Instituto de
Astrofísica de Canarias, E-38200 La Laguna, Tenerife, Spain  Nikolay K.
Nikolov Space Telescope Science Institute, 3700 San Martin Drive, Baltimore,
MD 21218, USA  Matthew C. Nixon Department of Astronomy, University of
Maryland, College Park, MD, USA  Seth Redfield Astronomy Department and Van
Vleck Observatory, Wesleyan University, Middletown, CT 06459, USA  Kevin B.
Stevenson Johns Hopkins APL, Laurel, MD, USA  Hannah R. Wakeford University of
Bristol, HH Wills Physics Laboratory, Tyndall Avenue, Bristol, UK  Munazza K.
Alam Carnegie Earth & Planets Laboratory, 5241 Broad Branch Road NW,
Washington, DC 20015, USA  Björn Benneke Department of Physics and Institute
for Research on Exoplanets, Universite de Montréal, Montreal, QC, Canada
Jasmina Blecic Department of Physics, New York University Abu Dhabi, Abu
Dhabi, UAE  Center for Astro, Particle and Planetary Physics (CAP3), New York
University Abu Dhabi, Abu Dhabi, UAE  Nicolas Crouzet Leiden Observatory,
Leiden University, P.O. Box 9513, 2300 RA Leiden, The Netherlands  Tansu
Daylan Department of Astrophysical Sciences, Princeton University, 4 Ivy Lane,
Princeton, NJ 08544  LSSTC Catalyst Fellow  Julie Inglis Division of
Geological and Planetary Sciences, California Institute of Technology,
Pasadena, CA, 91125  Laura Kreidberg Max Planck Institute for Astronomy,
Heidelberg, Germany  Dominique J.M. Petit dit de la Roche Département
d’Astronomie, Université de Genève, Sauverny, Switzerland  Jake D. Turner NHFP
Sagan Fellow  Department of Astronomy and Carl Sagan Institute, Cornell
University, Ithaca, New York 14853, USA

###### Abstract

Carbon monoxide was recently reported in the atmosphere of the hot Jupiter
WASP-39b using the NIRSpec PRISM transit observation of this planet, collected
as part of the JWST Transiting Exoplanet Community Early Release Science (JTEC
ERS) Program. This detection, however, could not be confidently confirmed in
the initial analysis of the higher resolution observations with NIRSpec G395H
disperser. Here we confirm the detection of CO in the atmosphere of WASP-39b
using the NIRSpec G395H data and cross-correlation techniques. We do this by
searching for the CO signal in the unbinned transmission spectrum of the
planet between 4.6 and 5.0 $\mu$m, where the contribution of CO is expected to
be higher than that of other anticipated molecules in the planet’s atmosphere.
Our search results in a detection of CO with a cross-correlation function
(CCF) significance of $6.6\sigma$ when using a template with only
${\rm{}^{12}C^{16}O}$ lines. The CCF significance of the CO signal increases
to $7.5\sigma$ when including in the template lines from additional CO
isotopologues, with the largest contribution being from ${\rm{}^{13}C^{16}O}$.
Our results highlight how cross-correlation techniques can be a powerful tool
for unveiling the chemical composition of exoplanetary atmospheres from
medium-resolution transmission spectra, including the detection of
isotopologues.

Exoplanet atmospheres (487) — Exoplanet atmospheric composition (2021) — Hot
Jupiters (753) — Astronomical methods (1043)

††facilities: JWST(NIRSpec)††software:  batman (Kreidberg, 2015), emcee
(Foreman-Mackey et al., 2013b), ExoTIC-LD (Grant & Wakeford, 2022), george
(Ambikasaran et al., 2015), petitRADTRANS (Mollière et al., 2019), PyAstronomy
(Czesla et al., 2019, https://github.com/sczesla/PyAstronomy), spectres
(Carnall, 2017), scipy (Virtanen et al., 2020), Tiberius (Kirk et al., 2017,
2018, 2019, 2021)

##  1 Introduction

The detection of carbon monoxide ($\rm CO$) and water ($\rm H_{2}O$) in gas
giant exoplanet atmospheres has long been regarded as key to understanding how
these objects form. From atmospheric $\rm CO$ and $\rm H_{2}O$ measurements it
is possible to estimate the planets carbon-to-oxygen (C/O) ratios. The
observed C/O ratios can then help trace the location in the stellar
protoplanetary disk where the planets formed, and how they accreted the solids
and gas that compose their cores and their atmospheres (see e.g., Öberg et
al., 2011; Mordasini et al., 2016; Madhusudhan et al., 2017; Cridland et al.,
2019). Isotopologue ratios have also been suggested as good tracers of planet
formation conditions and atmospheric evolution (e.g., Clayton & Nittler, 2004;
Mollière & Snellen, 2019).

$\rm H_{2}O$ has been readily detectable in planetary atmospheres for years,
given the Hubble Space Telescope ($\rm HST$) capacity to observe one of its
main bands centered around 1.4 $\mu$m (e.g., Deming et al., 2013; Crouzet et
al., 2014; Wakeford et al., 2018; Benneke et al., 2019). $\rm CO$, on the
other hand, has been harder to detect because of the lack of suitable space-
based instrumentation targeting the spectral ranges of interest. Nevertheless,
a few detections have been possible from the ground for exoplanets around very
bright stars, using high-resolution observations (R $\sim$ $\rm
30,000-100,000$) centered on the $\rm CO$ 2.3 $\mu$m band (e.g., Snellen et
al., 2010; Rodler et al., 2012; Brogi et al., 2012; de Kok et al., 2013; Brogi
et al., 2014; van Sluijs et al., 2022), and for young exoplanets via direct
imaging (e.g., Konopacky et al., 2013; Barman et al., 2015; Petit dit de la
Roche et al., 2018). Recent studies have also reported
${\rm{}^{12}C^{16}O/^{13}C^{16}O}$ ratios in the atmospheres of a young
accreting super-jupiter (Zhang et al., 2021b) and a brown dwarf (Zhang et al.,
2021a), using high-resolution observations centered on the CO 2.3$\mu$m band.

With the successful launch of $\rm JWST$, it is now possible to observe the CO
band at 2.3 $\mu$m and its fundamental band centered at 4.7 $\mu$m with
several of $\rm JWST$’s onboard instruments. This, combined with $\rm JWST$’s
capacity to observe fainter targets than previous missions could, opens the
possibility to search for $\rm CO$ in a large number of exoplanets, enabling
future composition population studies.

Early $\rm JWST$ observations have already produced detections of $\rm CO$ in
exoplanetary atmospheres. Specifically, Miles et al. (2022) detected the 2.3
$\mu$m and 4.7 $\mu$m CO bands in the young, directly imaged gas giant planet
VHS 1256-1257 b using the NIRSpec IFU setup, and Rustamkulov et al. (2023)
reported a detection of CO at 4.7 $\mu$m in the transiting gas giant exoplanet
WASP-39b using NIRSpec PRISM 111Low-resolution (R $\sim$ $20-300$) double-pass
prism disperser. (R $\sim$ $20-300$). This last CO detection could not be
confirmed by the initial analysis of the observations of WASP-39b with NIRSpec
G395H 222High-resolution (R $\sim$ $2700$) grating disperser. (Alderson et
al., 2023), although a more recent analysis by Grant et al. (2023) has
detected CO in these data using a new technique that searches for the CO
signal in pre-defined sets of wavelength sub-bands.

WASP-39b (Faedi et al., 2011) is a 0.28 $\rm M_{Jup}$, 1.27 $\rm R_{Jup}$
exoplanet with an equilibrium temperature of about 1200 K, orbiting a G7 V
star. Four transits of WASP-39b were observed as part of the JWST Director’s
Discretionary Time Transiting Exoplanet Community Early Release Science
Program (ERS-1366 PI: N. Batalha; Stevenson et al., 2016; Bean et al., 2018)
with the NIRSpec PRISM and G395H, NIRCam WFSS (Wide Field Slitless
Spectroscopy), and NIRISS SOSS (Single Object Slitless Spectroscopy)
instrument modes, which combined detected a number of chemical species,
including Na ($19\sigma$), H2O ($33\sigma$), CO2 ($28\sigma$), CO ($7\sigma$)
(Rustamkulov et al., 2023; Alderson et al., 2023; Ahrer et al., 2022;
Feinstein et al., 2022) and the unexpected photochemical byproduct SO2
($4.8\sigma$) (Tsai et al., 2023; Alderson et al., 2023). The detection of CO
and CO2, combined with the absence of CH4 in the transmission spectrum of
WASP-39b, reveals that the planet has a super-solar metallicity.

Here we set out to confirm the detections of $\rm CO$ in Rustamkulov et al.
(2023) and Grant et al. (2023) and search for CO isotopologues using the JTEC
ERS Program NIRSpec G395H dataset and cross-correlation techniques. Cross-
correlation techniques have been ubiquitously used in astronomy for decades
(e.g., Tonry & Davis, 1979), but they have only been used for exoplanet
atmospheres studies through transmission spectroscopy at high ($R\sim
100,000$) resolutions (e.g., Snellen et al., 2010) or direct imaging at medium
($R\sim 1,000$) resolutions (e.g., Konopacky et al., 2013). The superb
stability and data quality of JWST observations allows us to apply this
technique to medium resolution transmission spectra.

Figure 1: Transmission spectrum of WASP-39b obtained from the JTEC ERS Program
NIRSpec G395H NRS2 observation using the Tiberius pipeline. In the lower part
of the plot and in the right y-axis, we show the relative opacity
contributions of $\rm H_{2}O$, $\rm CO$ and $\rm CO_{2}$, which are the main
features expected in that wavelength range. The vertical, purple-shaded area
highlights the part of the spectrum where the contribution of $\rm CO$ is
expected to be predominant over $\rm H_{2}O$ and $\rm CO_{2}$. That is the
portion of the spectrum used in the CCF analysis.

##  2 Observations and data reduction

JWST observed a full transit of WASP-39b with the NIRSpec G395H mode (Jakobsen
et al., 2022; Birkmann et al., 2022) as part of the JTEC ERS Program. The
target was monitored for 8 hours and 36 min between 30 July 2022 21:45 and 31
July 2022 06:21 UTC, using the Bright Object Time Series mode (BOTS), with the
1.6$\arcsec$ $\times$ 1.6$\arcsec$ slit, SUB2048 subarray, and the NRSRAPID
readout pattern. The observation resulted in a total of 465 integrations, with
70 groups per integration, with the spectra of the target dispersed across the
NRS1 and NRS2 detectors. A particular feature of this dataset is a tilt event
that produced a mid-transit jump in flux (see Alderson et al., 2023, for
details).

We reduced the NRS1 and NRS2 spectra separately using the Tiberius pipeline
(Kirk et al., 2017, 2018, 2019, 2021), but here we focus on the NRS2 detector,
which covers the CO fundamental band. Tiberius works from the group-level 2D
rateints images produced by the STScI Science Calibration Pipeline (Rigby et
al., 2022, v1.6.2), from which we remove the 1/f noise contributions (Jakobsen
et al., 2022; Birkmann et al., 2022) per column using the median of all off-
target pixels in each column. Images were corrected for bad and saturated
pixels manually, using custom-made masks that flagged all pixels deviating
more than 5$\sigma$ from running medians along pixel rows and replacing their
count values by interpolating the counts from nearby pixels. We traced the
spectra fitting Gaussian distributions at each column, and then extracted each
spectrum using 4-pixel wide apertures centered on the trace. Residual
background flux not captured by the 1/f correction was removed by fitting
linear polynomials along each column of the detector after masking a 3-pixel-
wide region centered on the trace. To remove residual cosmic rays in each
stellar spectrum, we compared each spectrum to the median of all the spectra,
iteratively removing 8$\sigma$ or larger outliers and replacing them by linear
interpolation with neighboring pixels. Finally, we corrected for any potential
wavelength shift in each spectrum by cross-correlating them with the first
spectrum and putting all the spectra in a common wavelength grid by linearly
resampling them. We used wavelength maps from the JWST Calibration Pipeline to
produce the wavelength solution, as in Alderson et al. (2023).

We produced a white light curve (WLC) from the extracted spectra by summing
the flux per spectrum over the full wavelength range of the detector
(3.823–5.177 $\mu$m), and fitted it with a transit + Gaussian Processes (GP)
systematics model using batman (Kreidberg, 2015), george (Ambikasaran et al.,
2015), and a Monte Carlo Markov Chain (MCMC) implemented using emcee (Foreman-
Mackey et al., 2013a). The model fitted for the planet-to-star radii ratio
(${\rm R_{p}/R_{s}}$), orbital inclination (${\rm i}$), epoch of transit
(${\rm T_{0}}$), semimajor axis in units of the stellar radius (${\rm
a/R_{s}}$), the ${\rm u_{1}}$ quadratic law limb-darkening coefficient, and an
exponential squared GP kernel with the x- and y-pixel shifts, the FWHM of the
spectra, and sky background as parameters. We used a Gaussian prior centered
on 0.015 to fit for ${\rm a/R_{s}}$, and uniform priors for all other fitted
parameters. For this fit, we kept the orbital period fixed to the value
reported by Mancini et al. (2018) and the eccentricity fixed to zero. We fixed
the ${\rm u_{2}}$ limb darkening coefficient to the value computed by ExoTIC-
LD for the quadratic limb darkening law.

For the spectroscopic light curves (SLC), we generated 2044 light curves at
1-pixel resolution within the NRS2 wavelenth range and fitted them using a
similar procedure to the one described above, but holding ${\rm i}$, ${\rm
T_{0}}$ and ${\rm a/R_{s}}$ fixed to the values obtained from the WLC fit, and
applied a systematics correction from the WLC fit to aid in fitting the mirror
tilt event. We then fitted each spectroscopic light curve using polynomials to
describe residual systematics, not GPs. We chose a combination of four linear
polynomials describing the change in x-position and y-position of the star on
the detector, FWHM of the stellar trace and the measured background flux.

Table A.1 summarizes the adopted system parameters, the priors for the
parameters fitted in the white and spectroscopic light curves, and the
posteriors for the WLC fit. Figure 1 shows the obtained transmission spectrum
of WASP-39b between 3.823–5.177 $\mu$m. The data in the figure are provided in
Table 1.

We note that in our analysis we use the spectrum extracted per pixel column,
not per resolution element as defined by the Nyquist sampling limit333JWST
User Documentation (JDox). Baltimore, MD: Space Telescope Science Institute;
https://jwst-docs.stsci.edu, since both samplings produce consistent results,
but the per pixel column sampling is more sensitive to the signal of CO and
its isotopologues as demonstrated in Appendix B.

Table 1: Wavelength range per point, transit depth, TD, and transit depth uncertainty, $\Delta$TD, of the WASP-39b transmission spectrum derived using Tiberius. The full table is available electronically in the online version of this article. $\lambda_{1}$ ($\mu$m) |  $\lambda_{2}$ ($\mu$m) | TD (ppm) |  $\Delta$TD (ppm)
---|---|---|---
3.822427 | 3.823099 | 0.02095 | 0.00044
3.823099 | 3.823770 | 0.02100 | 0.00045
3.823770 | 3.824442 | 0.02143 | 0.00040
… | … | … | …

##  3 CO Template Generation

The fundamental band of CO covers air wavelengths between 4.28 and 6.50
$\mu$m. The portion of the band in the rest frame of WASP-39b covered by the
NIRSpec NRS2 detector is highlighted by the purple line in Fig. 1. For our
search, however, we limit the spectral range to 4.6 – 5.0 $\mu$m to minimize
contamination from ${\rm CO_{2}}$ below 4.6 $\mu$m and ${\rm H_{2}O}$ above
5.0 $\mu$m, where the signals from those molecules dominate.

WASP-39 has a barycentric systemic velocity of $\rm-58.4421$ ${\rm
km~{}s^{-1}}$ (Mancini et al., 2018), and the JWST barycentric velocity at the
time of the observations was -28.888 ${\rm km~{}s^{-1}}$, as recorded in the
image headers under the keyword velosys. Therefore, any atmospheric ${\rm CO}$
signal from WASP-39b is expected to appear shifted by about -87.330 ${\rm
km~{}s^{-1}}$.

We generated pure CO atmospheric transmission spectrum templates of WASP-39b
with petitRADTRANS (Mollière et al., 2019), using the lbl (line-by-line)
radiative transfer mode with a resolving power of
$\lambda/\Delta\lambda=10^{6}$, and adopting both the main $\rm{}^{12}C^{16}O$
isotopologue and the all CO isotopes444 i.e.,
12C16O,13C16O,12C18O,12C17O,13C18O and 13C17O following abundances of
$98.65~{}\%$, $1.11~{}\%$, $0.20~{}\%$, $3.7\times 10^{-2}~{}\%$, $2.2\times
10^{-3}~{}\%$ and $4.1\times 10^{-4}~{}\%$, respectively. high-resolution line
lists ”CO_main_iso” and ”CO_all_iso” from HITEMP (Rothman et al., 2010). For
the planet, we adopted a radius of 1.27 $\rm R_{Jup}$, a surface gravity of
407.38 cm $\rm s^{-2}$ following Faedi et al. (2011), an isothermal
temperature profile of 1166 K (Mancini et al., 2018), and atmospheric bottom
and top pressures of $10^{2}$ and $10^{-12}$ bars, respectively.

The top panels of Figure 2 show the resultant R=$\rm 10^{6}$ ${\rm CO}$ model
spectra for the main CO isotopologue (left) and for all CO isotopologues
(right). To match the resolution of the observations, we resampled the
template ${\rm CO}$ spectrum produced by petitRADTRANS to the sampling of the
observations using spectres555https://spectres.readthedocs.io/ (Carnall,
2017), which produced the spectra shown in the middle panel of Fig. 2. The
bottom panel of Fig. 2 shows the resampled templates, normalized using a 4th-
order polynomial fit. The portion of these templates between 4.6 and 5.0
$\rm\mu$m, highlighted in purple, is the one used in the cross-correlation
search for CO.

As a final step before starting the cross-correlation search, we normalized
the observed transmission spectrum of WASP-39b between 4.6 and 5.0 $\mu$m by
subtracting from the spectrum its mean value in that wavelength range.

Figure 2: Atmospheric models for the main CO isotopologue (a) and all CO
isotopologues (b). Top: High-resolution ($10^{6}$) petitRADTRANS-computed
transmission spectrum model of an atmosphere composed solely of CO. Middle: CO
spectrum model degraded to match the resolution of the observed transmission
spectrum. The region used for the cross-correlation search is highlighted in
purple, with resolving powers of ${\rm R=3470}$ and ${\rm 3803}$ between 4.6
and 5.0 $\rm\mu$m, respectively. The black dashed line shows the 4-order
polynomial fit that is used for the normalization. Bottom: CO template used in
the cross-correlation analysis, calculated by subtracting the 4-order
polynomial fit to the ${\rm R=3470-3803}$ CO model.

##  4 Cross-Correlation Search Results

We carried out cross-correlations between the normalized observed spectrum and
the normalized CO main-isotope and CO all-isotopes templates using the
crosscorr function in PyAstronomy (Czesla et al.,
2019)666https://github.com/sczesla/PyAstronomy, shifting the templates with
respect to the observed spectrum over a range of systemic velocities between
$\pm$ 1000 ${\rm km~{}s^{-1}}$, in 6 ${\rm km~{}s^{-1}}$ steps, which is about
15 times smaller than the pixel-to-pixel radial-velocity sampling of the
spectrum. We also tested the cross-correlation using a smaller, 1 ${\rm
km~{}s^{-1}}$ step, obtaining analogous results. We calculated uncertainties
for each correlation value (CV) by propagating the uncertainties in the
observed transmission spectrum of WASP-39b as

$\rm\Delta CV(RV)=\sum_{\lambda}\Delta S(\lambda)T(RV,\lambda)~{},$ (1)

where $\rm\Delta CV(RV)$ is the uncertainty of each point of the cross-
correlation function (CCF), $\rm\Delta S(\lambda)$ is the uncertainty of each
point of the transmission spectrum, and $\rm T(RV,\lambda)$ is the template
used in the cross-correlations at each RV shift.

To put the resulting CVs into more tangible units, we converted the CVs to
signal-to-noise ratio (SNR) values as follows: We defined out-of-peak
continuum regions where the CVs are expected to be dominated by noise. These
regions are between [$-452$, $-262$] ${\rm km~{}s^{-1}}$ and [$+98$, $+288$]
${\rm km~{}s^{-1}}$, defined as 190 ${\rm km~{}s^{-1}}$ wide regions on each
side of the central peak of the CCF, starting at velocities $\pm$ 3$\sigma$
away from that central peak, assuming that the distribution of velocities in
the peak is Gaussian, see Fig. 3. The conversion of the CVs to SNR values is
then calculated as

$\rm SNR(per~{}point)=\frac{|CV(per~{}point)-\mu_{cont}|}{\sigma_{cont}}~{},$
(2)

where $\mu_{cont}$ is the mean of the CVs in the continuum regions defined
above, and $\sigma_{cont}$ is the standard deviation of the CVs in those
regions.

The resultant CCFs in SNR units are shown in Figure 3, where the continuum
regions are indicated by the horizontal black lines in the top panels. The
CCFs for the CO main-isotopologue and the CO all-isotopologues templates both
show SNR $>$ 8 peaks centered near $-80~{}{\rm km~{}s^{-1}}$. This velocity
shift is consistent with the barycentric velocity of the system at the time of
the observations as discussed in Sect. 3.

Figure 3: Cross-correlation function (CCF) of the WASP39b transmission
spectrum with the CO main isotopologue (left panels) and the CO all
isotopologues templates (right panels) in the 4.6–5.0 $\rm\mu$m range using a
radial velocity step of 6 km/s (blue dotted line). The black horizontal bars
indicate the continuum regions used to estimate the noise of the CCF, and
calculate the significance of the central peak. In panels (a) and (c) we
compare the CCFs resulting from the self cross-correlation of the template
(light blue area) and the spectrum (dark blue area), which were scaled
arbitrarily. Panels (b) and (d) show the central peak of the CCF compared to
the histogram of the measurements of the CCF peak radial velocity through
random sampling. The right y-axis shows the percentage of samples. The red-
continuum line shows the mean value of the RV measurements distribution. The
black-dashed line shows the systemic velocity relative to the JWST at the time
of the observation. The light-blue lines show a set of Gaussian fits among the
sampled CCFs.

To measure the significance and the actual central velocity of the main peaks
in each CCF we fit a Gaussian to the SNR-converted CCF using the least squares
routine scipy.optimize.curvefit, and fit for the SNR value of the central
peak, ${\rm S_{peak}}$, its central velocity, ${\rm RV_{peak}}$, and the FWHM
of the distribution of velocities, ${\rm\sigma_{RVs}}$.

As a last step, we estimate uncertainties on the values of ${\rm S_{peak}}$
and ${\rm RV_{peak}}$, calculated as described above, by generating 30,000
Monte Carlo random samples of the observed transmission spectrum of WASP-39b
in Figure 1, where each point is simulated using a uniform probability
distribution of values within the range of errors observed in the spectrum. We
then calculate cross-correlations for each simulated spectrum, with each of
the two CO templates and measure the significance of the peak and central
velocity of each peak as described above. The bottom panels of Figure 3 show
the central CCF peaks for each CO template with a sample of the CCFs obtained
from individual MCMC realizations in light blue, and the resultant histograms
of the central radial velocities of the peak in red. The solid blue lines in
the figures correspond to the CCFs obtained from the observed spectrum, and
vertical dashed lines correspond to the velocity of the WASP-39 system at the
time of the observations. The retrieved distributions of the three parameters,
${\rm S_{peak}}$, ${\rm RV_{peak}}$, and ${\rm\sigma_{RV}}$ for each CO
template are shown in Figure 4.

Our analysis reveals CCF peaks centered at ${\rm-79.4\pm 9.8~{}km~{}s^{-1}}$
with a significance of ${\rm 6.6~{}_{-1.50}^{+2.02}~{}\sigma}$ when using the
CO main-isotopologue template, and at
${\rm-75.7~{}^{+9.2}_{-9.4}~{}km~{}s^{-1}}$ with a significance of ${\rm
7.5~{}_{-1.66}^{+2.30}~{}\sigma}$ when using the CO all-isotopologues
template. Both velocities are consistent within about $1\sigma$ with the
computed velocity of the system at the time of the observations of
${\rm-87.33~{}km~{}s^{-1}}$. In Appendix C, we show an additional analysis of
the observed transmission spectrum using a template composed of the two most
abundant CO isotopologues: $\rm{}^{12}C^{16}O$ and $\rm{}^{13}C^{16}O$. The
CCF results are similar to the ones obtained using the CO all-isotopologues
template. The significance of the detection, being $\rm 7.33^{+2.26}_{-1.04}$
sigma, increases compared to the CO main-isotopologue result but is slightly
less than the all-isotopologues significance, indicating that the all-
isotopologues template produces a better match.

###  4.1 Secondary Cross-Correlation Function Peaks

In addition to the central CCF peaks in agreement with the systemic velocity
of WASP-39, the top panels in Figure 3 show secondary CCF peaks at around -600
${\rm km~{}s^{-1}}$ and +600 ${\rm km~{}s^{-1}}$, which also appear
significant. To elucidate the origin of those peaks we cross-correlated the
observed transmission spectrum of WASP-39b with itself, and each CO template
with itself at the same resolution as the observations, over the same range of
velocities, ${\rm\pm~{}1000~{}km~{}s^{-1}}$, and with the same ${\rm
6~{}km~{}s^{-1}}$ steps. The top panels in Figure 3 show the result of the
template-to-template CCFs in cyan, and of the data-to-data CCF in light blue.
The template-to-template and data-to-data CCFs have been shifted by the ${\rm
RV_{peak}}$ values measured in the previous section so we can compare them
directly to the data-to-template CCFs.

The template-to-template CCFs show the expected peak at ${\rm
0~{}km~{}s^{-1}}$, and two symmetric peaks at a mean distance of
${\rm\pm~{}624~{}km~{}s^{-1}}$ from the central peak, similar to the ones
observed in the data-to-template CCFs. These peaks are expected in the
template-to-template CCFs of CO because of the periodic patterns in the
distribution of lines in the spectrum of this molecule (see Figure 2). The
data-to-data CCF shows the expected peak at ${\rm 0~{}km~{}s^{-1}}$, but no
other significant peaks. The presence of the same secondary peaks in the data-
to-template and template-to-template CCFs reinforces the conclusion that the
signal detected in the transmission spectrum of WASP-39b comes from CO.

Figure 4: Distribution of the best-fit Gaussian parameters (significance of
the central peak, ${\rm S_{peak}}$, the central peak velocity, ${\rm
RV_{peak}}$, and the FWHM of the distribution of velocities,
${\rm\sigma_{RVs}}$) retrieved from the sampled transmission spectra that were
randomly generated following a Gaussian distribution centered in the observed
transmission spectrum and considering the errors as the Gaussian standard
deviation in each point. We performed 30000 iterations. Panels (a) and (b)
show the results using the main CO isotopologue and all CO isotopologues
templates, respectively.

##  5 Conclusions and Discussion

We report the detection of CO in the atmosphere of WASP-39b via a cross-
correlation analysis of the JTEC ERS Program NIRSpec G395H dataset. Our
analysis produces a 6.6$\sigma$ CCF detection of the fundamental band of CO
centered around 4.8 $\mu$m when using a 12C16O template and a 7.5$\sigma$
detection when using a template with all CO isotopologues included in
petitRADTRANS.

Our reported detection significance levels are relative to custom-designed
baseline regions of the CCF, and as such cannot be readily used to estimate CO
abundances. In addition, some stellar CO absortion is expected in WASP-39b’s
host star, as pointed out by Rustamkulov et al. (2023), and this may dilute
the strength of the planetary CO signal, yielding to underestimated detection
significances, an effect previously reported by Deming & Sheppard (2017).

Our result confirms the detections of CO in the atmosphere of WASP-39b
reported by Rustamkulov et al. (2023) and Grant et al. (2023). This detection,
combined with the detection of ${\rm CO_{2}}$ (JWST Transiting Exoplanet
Community Early Release Science Team et al., 2023), and the absence of ${\rm
CH_{4}}$ (Alderson et al., 2023), reinforces the conclusion of super-solar
atmospheric metallicity for this planet (Rustamkulov et al., 2023; Alderson et
al., 2023; Ahrer et al., 2022; Feinstein et al., 2022; Tsai et al., 2023).

In addition, we show that standard CCF analysis techniques can be used to
identify the presence of molecular features in exoplanet atmospheres using
transmission spectra with a resolving power of just a few thousand, as is the
case of JWST observations in high-resolution modes, provided that the spectra
are of high-enough quality. CO has the advantage that, at the native
resolution of JWST, it still presents a band shape with resolved lines. Still,
this technique can be used to search for further molecules in the atmosphere
of WASP-39b and also to identify atmospheric molecular species in other
planets observed with JWST using these modes.

We highlight the presence of secondary peaks in the CCFs described in section
4.1, which we also attribute to CO (see Appendix B). Even though we did not
exploit the idea in this letter, it is worth exploring if secondary peaks in
CCFs can be used to increase the confidence of future CO detections, or other
molecules, or to place stronger non-detection constraints.

Finally, we highlight how the significance of our CO detection increases from
6.6$\sigma$ to 7.5$\sigma$ when adding a number of isotopologues to the 12C16O
template. To evaluate the influence of individual isotopologues, we performed
an additional cross-correlation analysis using a $\rm{}^{12}C^{16}O$ ($\rm
98.7\%$) + $\rm{}^{13}C^{16}O$ ($\rm 1.3\%$) template. In this case the
significance of the CCF detection is 7.3$\sigma$, indicating that, as
expected, 13C16O is the second most common form of CO. Although the
isotopologue detections are admittedly marginal, these results suggest that
with future higher signal-to-noise NIRSpec data it will be possible to search
for individual isotopologues of CO or other molecules by looking at how the
significance of CCF detections increases when adding one isotopologue at a
time.

##  6 Acknowledgments

This work is based on observations made with the NASA/ESA/CSA JWST. The data
were obtained from the Mikulski Archive for Space Telescopes at the Space
Telescope Science Institute, which is operated by the Association of
Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127
for JWST. The specific observations can be accessed via
http://dx.doi.org/10.17909/gdm2-0q65 (catalog doi:10.17909/gdm2-0q65). These
observations are associated with program JWST-ERS-01366. Support for program
JWST-ERS-01366 was provided by NASA through a grant from the Space Telescope
Science Institute. This work was supported by grant JWST-ERS-01366.033-A.
E.E-B and E.P. acknowledge funding from the Spanish Ministry of Economics and
Competitiveness through project PGC2018-098153-B-C31. E.E-B acknowledges
financial support from the European Union and the State Agency of
Investigation of the Spanish Ministry of Science and Innovation (MICINN) under
the grant PRE2020-093107 of the Pre-Doc Program for the Training of Doctors
(FPI-SO) through FSE funds. J.K. acknowledges financial support from Imperial
College London through an Imperial College Research Grant.

## References

  * Ahrer et al. (2022) Ahrer, E.-M., Stevenson, K. B., Mansfield, M., et al. 2022, arXiv e-prints, arXiv:2211.10489.  https://arxiv.org/abs/2211.10489
  * Alderson et al. (2023) Alderson, L., Wakeford, H. R., Alam, M. K., et al. 2023, Nature, 614, 664, doi: 10.1038/s41586-022-05591-3
  * Ambikasaran et al. (2015) Ambikasaran, S., Foreman-Mackey, D., Greengard, L., Hogg, D. W., & O’Neil, M. 2015, IEEE Transactions on Pattern Analysis and Machine Intelligence, 38, 252, doi: 10.1109/TPAMI.2015.2448083
  * Barman et al. (2015) Barman, T. S., Konopacky, Q. M., Macintosh, B., & Marois, C. 2015, ApJ, 804, 61, doi: 10.1088/0004-637X/804/1/61
  * Bean et al. (2018) Bean, J. L., Stevenson, K. B., Batalha, N. M., et al. 2018, PASP, 130, 114402, doi: 10.1088/1538-3873/aadbf3
  * Benneke et al. (2019) Benneke, B., Wong, I., Piaulet, C., et al. 2019, ApJ, 887, L14, doi: 10.3847/2041-8213/ab59dc
  * Birkmann et al. (2022) Birkmann, S. M., Ferruit, P., Giardino, G., et al. 2022, A&A, 661, A83, doi: 10.1051/0004-6361/202142592
  * Brogi et al. (2014) Brogi, M., de Kok, R. J., Birkby, J. L., Schwarz, H., & Snellen, I. A. G. 2014, A&A, 565, A124, doi: 10.1051/0004-6361/201423537
  * Brogi et al. (2012) Brogi, M., Snellen, I. A. G., de Kok, R. J., et al. 2012, Nature, 486, 502, doi: 10.1038/nature11161
  * Carnall (2017) Carnall, A. C. 2017, arXiv e-prints, arXiv:1705.05165.  https://arxiv.org/abs/1705.05165
  * Clayton & Nittler (2004) Clayton, D. D., & Nittler, L. R. 2004, ARA&A, 42, 39, doi: 10.1146/annurev.astro.42.053102.134022
  * Cridland et al. (2019) Cridland, A. J., van Dishoeck, E. F., Alessi, M., & Pudritz, R. E. 2019, A&A, 632, A63, doi: 10.1051/0004-6361/201936105
  * Crouzet et al. (2014) Crouzet, N., McCullough, P. R., Deming, D., & Madhusudhan, N. 2014, ApJ, 795, 166, doi: 10.1088/0004-637X/795/2/166
  * Czesla et al. (2019) Czesla, S., Schröter, S., Schneider, C. P., et al. 2019, PyA: Python astronomy-related packages.  http://ascl.net/1906.010
  * de Kok et al. (2013) de Kok, R. J., Brogi, M., Snellen, I. A. G., et al. 2013, A&A, 554, A82, doi: 10.1051/0004-6361/201321381
  * Deming & Sheppard (2017) Deming, D., & Sheppard, K. 2017, ApJ, 841, L3, doi: 10.3847/2041-8213/aa706c
  * Deming et al. (2013) Deming, D., Wilkins, A., McCullough, P., et al. 2013, ApJ, 774, 95, doi: 10.1088/0004-637X/774/2/95
  * Faedi et al. (2011) Faedi, F., Barros, S. C. C., Anderson, D. R., et al. 2011, A&A, 531, A40, doi: 10.1051/0004-6361/201116671
  * Feinstein et al. (2022) Feinstein, A. D., Radica, M., Welbanks, L., et al. 2022, arXiv e-prints, arXiv:2211.10493.  https://arxiv.org/abs/2211.10493
  * Foreman-Mackey et al. (2013a) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013a, PASP, 125, 306, doi: 10.1086/670067
  * Foreman-Mackey et al. (2013b) —. 2013b, PASP, 125, 306, doi: 10.1086/670067
  * Grant & Wakeford (2022) Grant, D., & Wakeford, H. R. 2022, Exo-TiC/ExoTiC-LD: ExoTiC-LD v3.0.0, v3.0.0, Zenodo, doi: 10.5281/zenodo.7437681
  * Grant et al. (2023) Grant, D., Lothringer, J. D., Wakeford, H. R., et al. 2023, arXiv e-prints, arXiv:2304.11994, doi: 10.48550/arXiv.2304.11994
  * Jakobsen et al. (2022) Jakobsen, P., Ferruit, P., Alves de Oliveira, C., et al. 2022, A&A, 661, A80, doi: 10.1051/0004-6361/202142663
  * JWST Transiting Exoplanet Community Early Release Science Team et al. (2023) JWST Transiting Exoplanet Community Early Release Science Team, Ahrer, E.-M., Alderson, L., et al. 2023, Nature, 614, 649, doi: 10.1038/s41586-022-05269-w
  * Kirk et al. (2019) Kirk, J., López-Morales, M., Wheatley, P. J., et al. 2019, AJ, 158, 144, doi: 10.3847/1538-3881/ab397d
  * Kirk et al. (2017) Kirk, J., Wheatley, P. J., Louden, T., et al. 2017, MNRAS, 468, 3907, doi: 10.1093/mnras/stx752
  * Kirk et al. (2018) —. 2018, MNRAS, 474, 876, doi: 10.1093/mnras/stx2826
  * Kirk et al. (2021) Kirk, J., Rackham, B. V., MacDonald, R. J., et al. 2021, AJ, 162, 34, doi: 10.3847/1538-3881/abfcd2
  * Konopacky et al. (2013) Konopacky, Q. M., Barman, T. S., Macintosh, B. A., & Marois, C. 2013, Science, 339, 1398, doi: 10.1126/science.1232003
  * Kreidberg (2015) Kreidberg, L. 2015, PASP, 127, 1161, doi: 10.1086/683602
  * Madhusudhan et al. (2017) Madhusudhan, N., Bitsch, B., Johansen, A., & Eriksson, L. 2017, MNRAS, 469, 4102, doi: 10.1093/mnras/stx1139
  * Mancini et al. (2018) Mancini, L., Esposito, M., Covino, E., et al. 2018, A&A, 613, A41, doi: 10.1051/0004-6361/201732234
  * Miles et al. (2022) Miles, B. E., Biller, B. A., Patapis, P., et al. 2022, arXiv e-prints, arXiv:2209.00620.  https://arxiv.org/abs/2209.00620
  * Mollière & Snellen (2019) Mollière, P., & Snellen, I. A. G. 2019, A&A, 622, A139, doi: 10.1051/0004-6361/201834169
  * Mollière et al. (2019) Mollière, P., Wardenier, J. P., van Boekel, R., et al. 2019, A&A, 627, A67, doi: 10.1051/0004-6361/201935470
  * Mordasini et al. (2016) Mordasini, C., van Boekel, R., Mollière, P., Henning, T., & Benneke, B. 2016, ApJ, 832, 41, doi: 10.3847/0004-637X/832/1/41
  * Öberg et al. (2011) Öberg, K. I., Murray-Clay, R., & Bergin, E. A. 2011, ApJ, 743, L16, doi: 10.1088/2041-8205/743/1/L16
  * Petit dit de la Roche et al. (2018) Petit dit de la Roche, D. J. M., Hoeijmakers, H. J., & Snellen, I. A. G. 2018, A&A, 616, A146, doi: 10.1051/0004-6361/201833384
  * Rigby et al. (2022) Rigby, J., Perrin, M., McElwain, M., et al. 2022, arXiv e-prints, arXiv:2207.05632.  https://arxiv.org/abs/2207.05632
  * Rodler et al. (2012) Rodler, F., Lopez-Morales, M., & Ribas, I. 2012, ApJ, 753, L25, doi: 10.1088/2041-8205/753/1/L25
  * Rothman et al. (2010) Rothman, L. S., Gordon, I. E., Barber, R. J., et al. 2010, J. Quant. Spec. Radiat. Transf., 111, 2139, doi: 10.1016/j.jqsrt.2010.05.001
  * Rustamkulov et al. (2023) Rustamkulov, Z., Sing, D. K., Mukherjee, S., et al. 2023, Nature, 614, 659, doi: 10.1038/s41586-022-05677-y
  * Snellen et al. (2010) Snellen, I. A. G., de Kok, R. J., de Mooij, E. J. W., & Albrecht, S. 2010, Nature, 465, 1049, doi: 10.1038/nature09111
  * Stevenson et al. (2016) Stevenson, K., Agol, E., Angerhausen, D., et al. 2016, A Preparatory Program to Identify the Single Best Transiting Exoplanet for JWST Early Release Science, HST Proposal. Cycle 24, ID. #14642
  * Tonry & Davis (1979) Tonry, J., & Davis, M. 1979, AJ, 84, 1511, doi: 10.1086/112569
  * Tsai et al. (2023) Tsai, S.-M., Lee, E. K. H., Powell, D., et al. 2023, Nature, 617, 483, doi: 10.1038/s41586-023-05902-2
  * van Sluijs et al. (2022) van Sluijs, L., Birkby, J. L., Lothringer, J., et al. 2022, arXiv e-prints, arXiv:2203.13234.  https://arxiv.org/abs/2203.13234
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Wakeford et al. (2018) Wakeford, H. R., Sing, D. K., Deming, D., et al. 2018, AJ, 155, 29, doi: 10.3847/1538-3881/aa9e4e
  * Zhang et al. (2021a) Zhang, Y., Snellen, I. A. G., & Mollière, P. 2021a, A&A, 656, A76, doi: 10.1051/0004-6361/202141502
  * Zhang et al. (2021b) Zhang, Y., Snellen, I. A. G., Bohn, A. J., et al. 2021b, Nature, 595, 370, doi: 10.1038/s41586-021-03616-x

##  Appendix A White and Spectroscopic light curves adopted parameters and
priors

Table A.1: Adopted parameters for the white light curve (WLC) fitting, including whether the value was fixed in the fitting or if a uniform ($\cal{U}$) or Gaussian ($\cal{G}$) prior was used instead. For uniform priors, the lower and upper bounds are given in brackets. For Gaussian priors, the mean and standard deviation of the priors are given in brackets. We also list the posterior distribution of each parameter after the fitting. The posteriors of the WLC fitting were used as starting values for the spectroscopic light curves (SLC) fitting. Parameter | WLC starting value(1) | WLC Prior type | WLC Posterior / SLC starting value | SLC Prior type
---|---|---|---|---
P [days] | 4.0552941 (34) | Fixed | – | Fixed
e | 0 | Fixed | – | Fixed
$\rm T_{0}$ [BJD] | 2 455 342.96913 (63) |  $\cal{U}$($T_{0}-2$, $T_{0}+2$) | 2459791.61206 (89) | Fixed
i [$\rm\deg$] | $87.32$ |  $\cal{U}$(80,90) | $87.81^{+0.20}_{-0.18}$ | Fixed
a/$R_{s}$ | $11.55$ |  $\cal{G}$(11.55,0.17) | $11.425\pm 0.015$ | Fixed
$\rm R_{p}/R_{s}$ | 0.14052 |  $\cal{U}$(0., 0.3) | 0.14636 ${}^{+0.00068}_{-0.00071}$ |  $\cal{U}$(0.,0.3)
(1) Values adopted from Mancini et al. (2018).

##  Appendix B Spectrum Sampling Tests

To determine whether to use the observed pixel-column sampled or Nyquist limit
sampled spectrum of WASP-39b as input to our CCF analysis in section 4, we
followed the same procedure described in that section, but replacing the
observed spectrum of WASP-39b by synthetic spectra.

We generated synthetic CO spectra of WASP-39b following similar steps to the
ones described in section 3. We generated a CO all isotopologues transmission
spectrum of WASP-39b with petitRADTRANS (Mollière et al., 2019) at resolving
power R=$10^{6}$ following as described in section 3, and then lowered the
resolution of the spectrum to match the 1-pixel resolution of our NIRSpec
G395H NRS2 observations using the spectres module. From that spectrum we
produced four different synthetic spectra: 1-pixel sampled; 2-pixel sampled
starting on the first pixel of the array, hereafter 2-pixel (0), 2-pixel
sampled starting on the second pixel of the array, hereafter 2-pixel (1), and
3-pixel sampled. We normalized each of the synthetic spectra using a fourth
order polynomial, as in section 3, then selected the portion of the spectrum
between 4.6 and 5.0 $\mu$m used in our CO search, and added noise using the
distribution of noise in the observed spectrum between those wavelengths. As a
last step we offsetted each synthetic spectrum in wavelength by -87.33 ${\rm
kms^{-1}}$ to match the expected systemic velocity.

We cross-correlated each synthetic spectrum with the CO main-isotopologue and
CO all-isotopologues templates described in Section 3, following the same
procedure as with the real data described in section 4. The results are
summarized in Figure 5. The figure shows that the CO signal is detected in all
cases and with a CCF shape that is very similar to the one produced by the
real data, including a main peak centered at the systemic velocity, and
secondary peaks at a mean distance of $\pm$ 624 ${\rm kms^{-1}}$ from the
central peak. In addition, the SNRs of the central peaks are consistently
larger in the cases when we use the all isotopologues template, confirming the
increase in signal that we see in the real data when we use a template with
multiple isotopologues.

The simulations show that 1-pixel sampling produces the same CCF results, but
with higher SNR than the other samplings. Therefore we decided to use 1-pixel
sampling in our analysis. We note that the peak SNRs of the CCFs in the
simulation are significantly larger that the values obtained in the real data,
but we attribute that difference to the simulated spectra only containing CO,
which in the real data, there are also contributions from ${\rm H_{2}O}$,
${\rm CO_{2}}$ lines, stellar ${\rm CO}$ lines (Deming & Sheppard, 2017), in
that wavelength region that can dilute the CO signal.

Figure 5: Cross-correlation functions (CCFs) resulting from the cross-
correlation of the synthetic transmission spectra with the CO main
isotopologue (left panels) and the CO all isotopologues templates (right
panels) in the 4.6–5.0 $\rm\mu$m range using a radial velocity step of 6 km/s.
The first row panels show the resulting CCFs for the 1-pixel sampled
transmission spectrum. The second and third rows panels show the resulting
CCFs for the 2-pixel sampled transmission spectra, starting in the first and
second pixel, respectively. The third row panels show the resulting CCFs for
the 3-pixel sampled transmission spectrum. The black horizontal bars indicate
the continuum regions used to calculate the signal-to-noise ratio of the CCF,
and calculate the significance of the central peak. The black-dashed line
shows the systemic velocity relative to the JWST at the time of the
observation. Additionally, the median signal-to-noise ratio value of the main
peak is printed in each panel.

##  Appendix C $\rm{}^{12}C^{16}O$ \+ $\rm{}^{13}C^{16}O$ cross-correlation
results

We performed an additional cross-correlation analysis of the observed
transmission spectrum using a template composed of the two most abundant CO
isotopologues: $\rm{}^{12}C^{16}O$ and $\rm{}^{13}C^{16}O$. We generated the
template following the procedure detailed in Section 3, setting the relative
abundances of the isotopologues to match those of Earth, i.e., $\rm 98.7\%$
$\rm{}^{12}C^{16}O$ and $\rm 1.3\%$ $\rm{}^{13}C^{16}O$. Figure 6 shows the
model spectrum used as template in this case. The CCF results, shown in Figure
7), are similar to the ones obtained using the CO all-isotopologues template.
The significance of the detection, being $\rm 7.33^{+2.26}_{-1.04}$ sigma,
increases compared to the CO main-isotopologues result but is slightly below
the all-isotopologues significance, indicating that the all-isotopologues
template produces a better match.

Figure 6: Atmospheric models for the combination of $\rm{}^{12}C^{16}O$ ($\rm
98.7\%$) and $\rm{}^{13}C^{16}O$ ($\rm 1.3\%$). Top: High-resolution
($10^{6}$) petitRADTRANS-computed transmission spectrum model. Middle: CO
spectrum model degraded to match the sampling of the observed transmission
spectrum. The region used for the cross-correlation search is highlighted in
purple. The black dashed line shows the 4-order polynomial fit that is used
for the normalization. Bottom: Normalized CO template used in the cross-
correlation analysis.  Figure 7: (a) Cross-correlation function (CCF) of the
WASP39b transmission spectrum with the $\rm{}^{12}C^{16}O$ ($\rm 98.7\%$) +
$\rm{}^{13}C^{16}O$ ($\rm 1.3\%$) template. We compare the CCFs resulting from
the self cross-correlation of the template (light blue area) and the spectrum
(dark blue area), which were scaled arbitrarily. (b) Distribution of the best-
fit Gaussian parameters (significance of the central peak, ${\rm S_{peak}}$,
the central peak velocity, ${\rm RV_{peak}}$, and the FWHM of the distribution
of velocities, ${\rm\sigma_{RVs}}$) retrieved from the sampled transmission
spectra that were randomly generated following a Gaussian distribution
centered in the observed transmission spectrum and considering the errors as
the Gaussian standard deviation in each point. We performed 30000 iterations.


# Probing CPT invariance with top quarks at the LHC

A. Belyaev,1,2, L. Cerrito,3,4 E. Lunghi,5 S. Moretti,1,2,6 and N. Sherrill7
1School of Physics and Astronomy, University of Southampton, Southampton, SO17
1BJ, United Kingdom  2Particle Physics Department, Rutherford Appleton
Laboratory, Chilton, Didcot, Oxon OX11 0QX, UK  3INFN Sezione di Roma Tor
Vergata, Rome 00133, Italy  4Dipartimento di Fisica, Universitá di Roma Tor
Vergata, Rome 00133, Italy  5Physics Department, Indiana University,
Bloomington, IN 47405, USA  6Department of Physics and Astronomy, Uppsala
University, Box 516, SE-751 20 Uppsala, Sweden  7Department of Physics and
Astronomy, University of Sussex, Brighton BN1 9QH, UK

###### Abstract

The first model-independent sensitivity to CPT violation in the top-quark
sector is extracted from ATLAS and CMS measurements of the top and antitop
kinematical mass difference. We find that the temporal component of a CPT-
violating background field interacting with the top-quark vector current is
restricted within the interval $[-0.13,0.29]$ GeV at 95% confidence level.

††preprint: APS/123-QED

CPT invariance is a symmetry of the Standard Model (SM) and generic prediction
of local, unitary, and Lorentz invariant quantum field theories in Minkowski
spacetime [1, *Pauli:1955, *Luders:1957bpq]. Among its consequences are the
equalities of particle and antiparticle masses and lifetimes.

Connecting theoretical predictions of CPT invariance violation to
particle–antiparticle CPT tests has been carried out for all particle species
of the SM with exception of the top quark. This work addresses this gap,
establishing the first sensitivity to top-sector CPT violation from top and
antitop mass reconstructions in ATLAS and CMS Collaboration measurements.

Kostelecký and Potting demonstrated that spontaneous CPT violation in string
theory may result in remnant observables (likely suppressed by the energy
scale of this violation) which can be possibly tested in current experiments
[4, *kp95]. The subsequent development of effective field theory (EFT)
descriptions [6, *ck98, *cg97, *cg98] generated intense interest in testing
fundamental spacetime symmetries and have produced hundreds of constraints
across a variety of systems [10]. Within this EFT framework, known as the
Standard Model Extension (SME), are operators involving top-quark fields.

Lorentz- and CPT-violating effects described by the SME can be viewed as
generalized tensor background fields. In view of the discovery of the Higgs
boson [11, 12] described by a homogeneous and isotropic scalar background with
$SU(2)\times U(1)$ quantum numbers, searching for new physics in the form of
additional backgrounds is well motivated.

In this study, we suggest using measurements of the kinematically
reconstructed top and antitop mass difference, $\Delta m^{\rm
kin}_{t\bar{t}}$, an observable which is uniquely sensitive to CPT violation,
to extract the first constraints on various coefficients responsible for it in
the top-quark sector of the SME. Dedicated measurements of the kinematically
reconstructed mass difference have been performed by DØ [13], CDF [14], ATLAS
[15], and CMS [16, 17].111In Ref. [18], earlier DØ and CDF top mass
measurements [19, 20] were used to place a conservative upper limit on the top
and antitop mass difference. All results are consistent with SM expectations
within the experimental uncertainties, and may be translated into constraints
on top quark coefficients for CPT violation.

In our study, we aim to test gauge invariant and renormalizable CPT-violating
SME operators of the form [6, 21]

$\displaystyle\mathcal{L}^{\text{CPT}-}=$ $\displaystyle-(a_{Q})_{\mu
AB}\bar{Q}_{A}\gamma^{\mu}Q_{B}-(a_{U})_{\mu AB}\bar{U}_{A}\gamma^{\mu}U_{B}$
$\displaystyle-(a_{D})_{\mu AB}\bar{D}_{A}\gamma^{\mu}D_{B}.$ (1)

The odd number of operator Lorentz indices implies a change of sign under a
CPT transformation. Therefore, $\mathcal{L}^{\text{CPT}-}$ is odd under CPT,
which is reflected in its superscript. We are interested in probing operators
involving the third generation ($A=B=3$), where

$Q_{3}=\begin{pmatrix}t\\\ b\end{pmatrix}_{L},\quad U_{3}=t_{R},\quad
D_{3}=b_{R}.$ (2)

The relevant coefficients for CPT-violating operators are therefore
$(a_{Q})_{\mu 33},(a_{U})_{\mu 33}$, and $(a_{D})_{\mu 33}$. Note that the
Lagrangian mass parameter $m_{t}$ remains identical for top and antitop quarks
in the presence of CPT violation, in accordance with Greenberg’s theorem [22].

A reduction of the number of independent coefficients in Eq. (Probing CPT
invariance with top quarks at the LHC) is possible under suitable
approximations. In the limit of a zero bottom quark mass $m_{b}$ (in
comparison with $m_{t}$), the phases of $D_{3}$ and $Q_{3}$ fields can be
independently changed. Thus, the position-dependent field redefinition
$D_{3}\rightarrow\exp[-i(a_{D})_{\mu 33}x^{\mu}]D_{3}$ can be used to remove
the last term in Eq. (Probing CPT invariance with top quarks at the LHC). An
analogous redefinition of $Q_{3}$ and $U_{3}$ with the same phase
$\exp[-i(a_{Q})_{\mu 33}x^{\mu}]$ allows to eliminate the term with $Q_{3}$
and shifts the coefficient $(a_{Q})_{\mu 33}$ into the $U_{3}$ term, such that
$(a_{U})_{\mu 33}\rightarrow(a_{U})_{\mu 33}-(a_{Q})_{\mu 33}$. This
transformation preserves the structure of all SM terms and yields an
equivalent Lagrangian density expressed in the mass eigenstate basis:

$\mathcal{L}_{\rm top}^{\text{CPT}-}=b_{\mu}\bar{t}_{R}\gamma^{\mu}t_{R},$ (3)

where $b_{\mu}\equiv[(a_{Q})_{\mu 33}-(a_{U})_{\mu 33}]$ and only $t_{R}$
fields appear. Under these field redefinitions, CPT violation is isolated to
the top-quark sector and quantified through the differential propagation of
$t_{R}$ relative to $t_{L}$. 222One could have chosen $\exp[-i(a_{U})_{\mu
33}x^{\mu}]$ as phase shift for both $Q_{3}$ and $U_{3}$ which would lead to
$\mathcal{L^{\prime}}_{\rm
top}^{\text{CPT}-}=b^{\prime}_{\mu}\bar{t}_{L}\gamma^{\mu}t_{L},$ (4) with
$b^{\prime}_{\mu}\equiv[(a_{U})_{\mu 33}-(a_{Q})_{\mu 33}]=-b_{\mu}$ and only
$t_{L}$ fields involved. For this equivalent formulation the experimental
limits on $b_{\mu}$ will be identical to those for Eq. (3).

Performing the variational procedure including the conventional top kinetic
terms results in a modified Dirac equation

$\left[i\not{\partial}+\tfrac{1}{2}(1-\gamma_{5})\not{b}-m_{t}\right]t=0.$ (5)

Plane-wave solutions imply a quartic equation in $p^{\mu}=(E_{t},\vec{p})$
with four distinct solutions linear in $b_{\mu}$:

$\displaystyle p^{2}=\begin{cases}m_{t}^{2}-p\cdot b\pm[(p\cdot
b)^{2}-m_{t}^{2}b^{2}]^{1/2}\hskip 5.69054pt\rm{(top)}\\\ m_{t}^{2}+p\cdot
b\pm[(p\cdot b)^{2}-m_{t}^{2}b^{2}]^{1/2}\hskip
5.69054pt\rm{(antitop)}\end{cases}$ (6)

neglecting higher-order terms in $b_{\mu}$. The first and second
unconventional terms in each row are associated with the vector and
pseudovector pieces of Eq. (3), respectively, and the $\pm$ signs denote
states of opposite helicities. The difference between particle and
antiparticle solutions is obtained by $b_{\mu}\rightarrow-b_{\mu}$, reflecting
the CPT-odd property of Eq. (3) and the effect of CPT conjugation on the
plane-wave solutions. Note that potential CPT-violating corrections to top and
antitop decay widths are suppressed relative to free-propagation effects (6)
by the square of the weak coupling constant and are neglected.

The presence of $b_{\mu}$ implies both charge- and helicity-dependent energy
eigenvalues. Top $p$ and antitop $\bar{p}$ momenta and kinematical masses
$m_{t}^{\rm{kin}}\equiv\sqrt{p^{2}},m_{\bar{t}}^{\rm{kin}}\equiv\sqrt{\bar{p}^{2}}$
are reconstructed through the charge and four-momentum conservation of final-
state decay products. In the conventional CPT-invariant case,
$m_{t}^{\rm{kin}}=m_{\bar{t}}^{\rm{kin}}=m_{t}$. The kinematical mass
difference

$\Delta m^{\rm
kin}_{t\bar{t}}(p,\lambda_{p},\bar{p},\lambda_{\bar{p}},m_{t},b)\equiv
m_{t}^{\rm{kin}}-m_{\bar{t}}^{\rm{kin}}$ (7)

parametrizes a CPT-violating asymmetry, where $\lambda_{p}$
($\lambda_{\bar{p}}$) are the top (antitop) helicities. In principle,
measurements of $\Delta m^{\rm kin}_{t\bar{t}}$ could be used to extract
$b_{\mu}$. However, this is generically non-trivial because $\Delta m^{\rm
kin}_{t\bar{t}}$ is time dependent. Implicit time dependence enters via event-
by-event reconstructions of the top and antitop four-momenta. Explicit time
dependence enters through $b_{\mu}$ directly since the relevant experiments
are performed in non-inertial Earth-based laboratories. As a result, the
coefficient $b_{\mu}$ is modulated as a function of the laboratory velocity
and rotation rate [23]. It is convenient and standard practice to introduce
the approximately inertial sun-centered frame (SCF) where the coefficients for
CPT violation carry indices $\mu=\\{T,X,Y,Z\\}$ and may be approximated as
constants [24, *Bluhm:2003un, *Kostelecky:2003fs, *Kostelecky:2002hh]. In this
setting, the leading laboratory signatures are given by single harmonics of
the Earth’s sidereal rotation frequency $\omega_{\oplus}\approx
2\pi/(23\;\rm{h}\;56\;{min})$. Despite these complications, we demonstrate
that suitably averaged observables permit analyses involving distinct proper
subsets of $\\{b_{T},b_{X},b_{Y},b_{Z}\\}$.

The ATLAS and CMS Collaborations have reported measurements of $\langle\Delta
m^{\rm kin}_{t\bar{t}}\rangle$ [15, 16, 17], where $\langle\rangle$ indicates
averaging over all events. ATLAS used a sample of $t\bar{t}$ events in the
single charged lepton + jets decay mode, selected from 4.7 fb-1 of $pp$
collisions at $\sqrt{s}=7$ TeV [15]. The data were regularly collected over
several months in 2011 [28]. The kinematical mass difference $\langle\Delta
m^{\rm kin}_{t\bar{t}}\rangle$ was determined from a maximum likelihood fit to
the per-event top and antitop candidate mass difference, reconstructed in the
ATLAS detector frame. The reported result

$\displaystyle\langle\Delta m^{\rm kin}_{t\bar{t}}\rangle^{\rm ATLAS}=0.67\pm
0.61_{\rm stat}\pm 0.41_{\rm syst}\;{\rm GeV}$ (8)

is consistent with zero within uncertainties. The measurement involved
averaging over several sidereal days and sampling of the full $t\bar{t}$ phase
space. This measurement therefore has negligible sensitivity to the top
polarizations and spatial components of $b_{\mu}$. The invariance of the
temporal component $b_{0}$ under the rotation connecting the ATLAS detector
and SCF frames implies $b_{0}=b_{T}$ and Eq. (7) takes the form

$\langle\Delta m^{\rm
kin}_{t\bar{t}}\rangle^{\rm{ATLAS}}\approx-\frac{b_{T}}{m_{t}}\frac{\langle
E_{t}+E_{\bar{t}}\rangle}{2},$ (9)

where $\langle E_{t}+E_{\bar{t}}\rangle$ is the average over the phase space
of the sum of top and antitop energies characteristic to the dataset. To
leading order in $b_{\mu}$ the energies $E_{t},E_{\bar{t}}$ are the
conventional eigenenergies.

The CMS analyses also selected events where one $W$ boson decays hadronically
and the other leptonically. The data were split into $\ell^{+}$ and $\ell^{-}$
samples. Using the ideogram likelihood method [29], $m_{t}^{\rm kin}$ and
$m_{\bar{t}}^{\rm kin}$ were reconstructed from the two samples, respectively,
in the CMS collider frame, from which their difference was obtained. The 2017
analysis used a data sample [30] corresponding to $pp$ collisions at
$\sqrt{s}=8~{}$TeV and an integrated luminosity of $19.6\pm 0.5$ fb-1. The
data were collected over several months in 2012, yielding [31]

$\displaystyle\langle\Delta m^{\rm kin}_{t\bar{t}}\rangle^{\rm CMS}=-0.15\pm
0.19_{\rm stat}\pm 0.09_{\rm syst}\;{\rm GeV}.$ (10)

Since the tops and antitops were selected from different events, CMS is
sensitive to a sum of the averages,
$\langle E_{t}\rangle$ \+ $\langle E_{\bar{t}}\rangle$, and therefore

$\langle\Delta m^{\rm
kin}_{t\bar{t}}\rangle^{\rm{CMS}}\approx-\frac{b_{T}}{m_{t}}\frac{\langle
E_{t}\rangle+\langle E_{\bar{t}}\rangle}{2}.$ (11)

|  $t\bar{t}$ |  $t\bar{t}\to\ell\nu jjb\bar{b}$ tot |  $t\bar{t}\to\ell\nu jjb\bar{b}$ fid
---|---|---|---
$\langle E_{t}+E_{\bar{t}}\rangle_{@7~{}\text{TeV}}$ |  706.3 |  708.9 |  658.4
$\langle E_{t}+E_{\bar{t}}\rangle_{@8~{}\text{TeV}}$ |  738.9 |  742.2 |  674.4
$\langle E_{t}+E_{\bar{t}}\rangle_{@13~{}\text{TeV}}$ |  878.8 |  883.7 |  725.2
$\langle<EMAIL_ADDRESS>|  892.5 |  898.7 |  729.1
Table 1: The values of $\langle E_{t}+E_{\bar{t}}\rangle$ in GeV at various
LHC energies for $pp\to t\bar{t}$ processes: a) $pp\to t\bar{t}$ ($t\bar{t}$);
b) $pp\to t\bar{t}\to\ell\nu jjb\bar{b}$ with no cuts applied (tot); c) $pp\to
t\bar{t}\to\ell\nu jjb\bar{b}$ with fiducial CMS cuts applied (fid). See
further details in the text. The events here are generated with CalcHEP.

In order to extract an upper bound on $b_{T}$, we need to calculate the
average $\langle E_{t}+E_{\bar{t}}\rangle$ for $t\bar{t}$ events in the
fiducial region considered by the two experiments. For the same set of cuts,
however, there will be no difference between the sum of the averages $\langle
E_{t}\rangle$ \+ $\langle E_{\bar{t}}\rangle$ and average of the sum $\langle
E_{t}+E_{\bar{t}}\rangle$. The evaluation of $\langle
E_{t}+E_{\bar{t}}\rangle$ for $t\bar{t}$ events at various center-of-mass
energies is performed with the aid of the Monte Carlo (MC) generator CalcHEP
[32] and cross-checked using Madgraph [33] (interfaced with Pythia8 [34] and
the detector simulator Delphes [35]). The results we find are summarized in
Table 1.
The value of $b_{T}$ including uncertainty reads

$\displaystyle b_{T}=$ $\displaystyle-\frac{2\langle\Delta m^{\rm
kin}_{t\bar{t}}\rangle m_{t}}{\langle E_{t}+E_{\bar{t}}\rangle}$
$\displaystyle\times\left[1\pm\sqrt{\left(\frac{\delta\langle\Delta m^{\rm
kin}_{t\bar{t}}\rangle}{\langle\Delta m^{\rm
kin}_{t\bar{t}}\rangle}\right)^{2}+\left(\frac{\delta\langle
E_{t}+E_{\bar{t}}\rangle}{\langle
E_{t}+E_{\bar{t}}\rangle}\right)^{2}}\right].$ (12)

We estimate the uncertainty ${\delta\langle E_{t}+E_{\bar{t}}\rangle}$ on the
predicted value $\langle E_{t}+E_{\bar{t}}\rangle$ from our MC simulations by
varying the factorization scale within $m_{t}/2$ and $2m_{t}$, as well as
using different PDF sets. We find it to be below 5% for each uncertainty. We
have also evaluated the effect of the top-quark width on the $\langle
E_{t}+E_{\bar{t}}\rangle$ value by comparing the results for $2\to 2$ versus
$2\to 6$ processes at parton level (the second versus the third column in
Table 1) and find an increase of only about $+0.5\%$ on $\langle
E_{t}+E_{\bar{t}}\rangle$ when the width is taken into account, which is
simply related to a kinematical effect. A more detailed analysis of these
uncertainties is not relevant for this study since the uncertainty on $b_{T}$
is completely dominated by the experimental uncertainty on
$\delta\langle\Delta{m_{t\bar{t}}^{\rm kin}}\rangle$, which is of order 100%.
Therefore, even a very conservative value of ${\delta\langle
E_{t}+E_{\bar{t}}\rangle}$ of order 10% will affect the overall uncertainty
for $b_{T}$ determination at the level of only 1%, as one can see from Eq.
(12).

Combining experimental statistical and systematic uncertainties in quadrature
and using Eq. (12), we find the following exclusion limits on $b_{T}$:

$\displaystyle b_{T}\in\begin{cases}[-1.10,0.41]\;\text{GeV}&\text{ATLAS @ 7
TeV}\\\ [-0.13,0.29]\;\text{GeV}&\text{CMS @ 8 TeV}\\\ \end{cases}\;.$ (13)

Outside of these intervals $b_{T}$ is excluded at 95% confidence level.

| 7 TeV | 8 TeV | 13 TeV | 13.6 TeV
---|---|---|---|---
| tot | fid | tot | fid | tot | fid | tot | fid
$\langle|\Delta_{Z}|\rangle$ | 72 | 68 | 76 | 70 | 97 | 79 | 100 | 80
$\langle|C_{X}|\rangle$ | 74 | 69 | 79 | 70 | 103 | 84 | 103 | 81
$\langle|S_{X}|\rangle$ | 418 | 329 | 451 | 361 | 590 | 416 | 603 | 405
Table 2: Phase space averages (both in the whole phase space and in the CMS
fiducial region) for the quantities described in the text in units of GeV. The
events here are generated with MadGraph and showered/hadronized with Pythia8.
Finally, Delphes is used to simulate detector effects.

The spatial components $b_{X,Y,Z}$ can be constrained with the same dataset
but require dedicated analyses. To good approximation, the laboratory-frame
coefficients for both ATLAS and CMS are related to the SCF coefficients via
the rotation

$\begin{pmatrix}b_{1}\\\ b_{2}\\\
b_{3}\end{pmatrix}=\begin{pmatrix}-b_{Z}\sin\chi\cos\psi+\cos(\omega_{\oplus}T_{\oplus})(b_{X}\cos\chi\cos\psi+b_{Y}\sin\psi)+\sin(\omega_{\oplus}T_{\oplus})(b_{Y}\cos\chi\cos\psi-
b_{X}\sin\psi)\\\
b_{Z}\cos\chi+\sin\chi[b_{X}\cos(\omega_{\oplus}T_{\oplus})+b_{Y}\sin(\omega_{\oplus}T_{\oplus})]\\\
-b_{Z}\sin\chi\sin\psi+\cos(\omega_{\oplus}T_{\oplus})(b_{X}\cos\chi\sin\psi-
b_{Y}\cos\psi)+\sin(\omega_{\oplus}T_{\oplus})(b_{Y}\cos\chi\sin\psi+b_{X}\cos\psi)\end{pmatrix},$
(14)

where the colatitude is $\chi=43.7^{\rm o}$ and the beamline orientation
north-of-east is $\psi=-11.3^{\rm o}$. The time $T_{\oplus}$ is identified
with the event time. Note that this choice of laboratory frame is identical
for ATLAS and CMS: the 2-axis points in the upward vertical direction and the
1-axis points towards (away from) the center of the LHC ring for ATLAS (CMS).

The kinematical mass difference for single event in terms of the SCF
coefficients is given by

$\displaystyle\Delta m_{t\bar{t}}^{\rm kin}$
$\displaystyle=-\frac{1}{2m_{t}}\left[b_{T}\Delta_{T}+b_{Z}\Delta_{Z}+\sum_{A=X,Y}b_{A}\Big{(}C_{A}\cos(\omega_{\oplus}T_{\oplus})+S_{A}\sin(\omega_{\oplus}T_{\oplus})\Big{)}\right]\;,$
(15) $\displaystyle\Delta_{T}$ $\displaystyle=E_{t}+E_{\bar{t}}\;,$ (16)
$\displaystyle\Delta_{Z}$
$\displaystyle=-\sin\chi\cos\psi(p_{1}+\bar{p}_{1})+\cos\chi(p_{2}+\bar{p}_{2})-\sin\chi\sin\psi(p_{3}+\bar{p}_{3})\;,$
(17) $\displaystyle C_{X}$
$\displaystyle=\cos\chi\cos\psi(p_{1}+\bar{p}_{1})+\sin\chi(p_{2}+\bar{p}_{2})+\cos\chi\sin\psi(p_{3}+\bar{p}_{3})\;,$
(18) $\displaystyle S_{X}$
$\displaystyle=-\sin\psi(p_{1}+\bar{p}_{1})+\cos\psi(p_{3}+\bar{p}_{3})\;,$
(19) $\displaystyle C_{Y}$ $\displaystyle=-S_{X}\;,$ (20) $\displaystyle
S_{Y}$ $\displaystyle=C_{X}\;,$ (21)

where $p_{i}$ ($\bar{p}_{i}$) are the top (antitop) three-momentum components
in the laboratory frame. The symmetry of the $pp$ collisions guarantees
$\langle p_{i}+\bar{p}_{i}\rangle=0$, and hence $\langle\Delta
m_{t\bar{t}}^{\rm kin}\rangle=-\frac{b_{T}}{2m_{t}}\langle
E_{t}+E_{\bar{t}}\rangle$ when averaging over the whole phase space. This
means that the coefficients $b_{X,Y,Z}$ can only be constrained by dedicated
analyses. Moreover, note that the coefficients $b_{X}$ and $b_{Y}$ induce
effects that vanish also when averaged over time.

A strategy to build an observable sensitive exclusively to $b_{Z}$ is to
consider

$\displaystyle\langle{\Delta m_{t\bar{t}}^{\rm
kin}}^{\prime}\rangle=\langle\Delta m_{t\bar{t}}^{\rm kin}\;{\rm
sgn}[\Delta_{Z}]\rangle=-\frac{b_{Z}}{2m_{t}}\langle|\Delta_{Z}|\rangle.$ (22)

In fact, $\langle\Delta_{T}\;{\rm sgn}[\Delta_{Z}]\rangle=0$ and time
averaging eliminates $b_{X,Y}$. This quantity can be measured by first
calculating $\Delta m_{t\bar{t}}^{\rm kin}\;{\rm sgn}[\Delta_{Z}]$ on an
event-by-event basis ($\Delta_{Z}$ is a simple function of the top and antitop
three-momenta) and then studying its distribution. Since it is not possible to
simply convert the measured $\langle{\Delta m_{t\bar{t}}^{\rm kin}}\rangle$
distribution into the corresponding $\langle{\Delta m_{t\bar{t}}^{\rm
kin}}^{\prime}\rangle$ one, a constraint on $b_{Z}$ can only be obtained via a
dedicated re-analysis of the data. The $b_{Z}$ sensitivity of this analysis is
straightforward to assess because, in absence of a signal, we expect the
experimental $\langle{\Delta m_{t\bar{t}}^{\rm kin}}\rangle$ and
$\langle{\Delta m_{t\bar{t}}^{\rm kin}}^{\prime}\rangle$ distributions to be
very similar. In order to get an upper limit on $b_{Z}$, one should consider
the statistical uncertainty of the 8 TeV CMS analysis but adopt the larger
ATLAS systematic uncertainty, which has been calculated from a scheme where
the top and antitop masses are determined simultaneously on a per-event basis,
and center the distribution on zero: $\langle{\Delta m_{t\bar{t}}^{\rm
kin}}^{\prime}\rangle\sim[0\pm 0.19\pm 0.41]\;{\rm GeV}$. Combining these
sensitivities with the MC results for $\langle|\Delta_{Z}|\rangle$ presented
in Table 2 we obtain the expected sensitivity

$\displaystyle|b_{Z}|_{\rm expected}\lesssim 4.6\;{\rm GeV}.$ (23)

For the transverse components, $b_{X}$ and $b_{Y}$, the situation is more
complicated because a sidereal-time analysis is required. A possible strategy
is to divide the sidereal period in a number $N$ of bins. Focusing on the
coefficient $b_{X}$, the kinematical mass difference for events in the $n$-th
bin is given by:

$\displaystyle\Delta m_{t\bar{t}}^{\rm kin}$
$\displaystyle=-\frac{b_{X}}{2m_{t}}\Delta_{X}^{(n)},$ (24)
$\displaystyle\Delta_{X}^{(n)}$
$\displaystyle=C_{X}\langle\cos(\omega_{\oplus}T_{\oplus})\rangle_{n}+S_{X}\langle\sin(\omega_{\oplus}T_{\oplus})\rangle_{n},$
(25)

where $\langle\rangle_{n}$ indicates the time average over bin $n$. For each
bin we can then proceed exactly as for $b_{Z}$ and average over all events in
the bin weighting $\Delta m_{t\bar{t}}^{\rm kin}$ by ${\rm
sgn}[\Delta_{X}^{(n)}]$. In order to obtain the sensitivity to the
coefficients $b_{X}$ and $b_{Y}$ we simply calculate the phase space average
$[\langle|C_{A}|\rangle+\langle|S_{A}|\rangle]/\sqrt{2}$, where $A=X,Y$, and
the factor $1/\sqrt{2}$ is the root-mean-square of the $\sin$ and $\cos$
functions. Following the same procedure as for the $b_{Z}$ case (but using the
appropriate averages from Table 2), we find the following expected of the
sensitivity to $b_{X,Y}$:

$\displaystyle|b_{X,Y}|_{\rm expected}\lesssim 0.8\;{\rm GeV}.$ (26)

To summarize, we have established the first model-independent sensitivity to
CPT violation in the top-quark sector using LHC measurements of the top and
antitop kinematical mass difference. The constraint on $b_{T}$ (13) is about
two orders of magnitude stronger than those expected from single-top
production [21]. We have also suggested dedicated analyses which would lead to
constraints on $|b_{Z}|$ (23) and $|b_{X,Y}|$ (26) of the same order.

Regarding prospects, let us note that an analysis of the the entire Run-2
dataset (140 fb-1), assuming a factor of two improvement on the systematic
uncertainty and taking into account the larger $\sqrt{s}=13\;{\rm TeV}$
center-of-mass energy, is expected to yield a sensitivity to $b_{T}$ at the
0.05 GeV level (with similar fractional improvements for $b_{X,Y,Z}$). Further
accumulation of data will lead to a systematically dominated total
uncertainty.

## Acknowledgments

AB and SM are supported in part through the NExT Institute and STFC
Consolidated Grant No. ST/L000296/1. SM is also supported by the Knut and
Alice Wallenberg foundation under the grant KAW 2017.0100. NS is supported in
part by the Science and Technology Facilities Council (grant number
ST/T006048/1). We thank V. A. Kostelecký for useful comments and discussions.

## References

  * Bell [1955] J. S. Bell, Proc. Roy. Soc. Lond. A 231, 479 (1955).
  * Pauli [1955] W. Pauli, in  _Niels Bohr and the Development of Physics_ , edited by W. Pauli (McGraw-Hill, New York, 1955) p. 30.
  * Luders [1957] G. Luders, Annals Phys. 2, 1 (1957).
  * Kostelecký and Potting [1991] V. A. Kostelecký and R. Potting, Nucl. Phys. B 359, 545 (1991).
  * Kostelecký and Potting [1995] V. A. Kostelecký and R. Potting, Phys. Rev. D 51, 3923 (1995), arXiv:hep-ph/9501341 .
  * Colladay and Kostelecký [1997] D. Colladay and V. A. Kostelecký, Phys. Rev. D 55, 6760 (1997), arXiv:hep-ph/9703464 .
  * Colladay and Kostelecký [1998] D. Colladay and V. A. Kostelecký, Phys. Rev. D 58, 116002 (1998), arXiv:hep-ph/9809521 .
  * Coleman and Glashow [1997] S. R. Coleman and S. L. Glashow, Phys. Lett. B 405, 249 (1997), arXiv:hep-ph/9703240 .
  * Coleman and Glashow [1999] S. R. Coleman and S. L. Glashow, Phys. Rev. D 59, 116008 (1999), arXiv:hep-ph/9812418 .
  * Kostelecký and Russell [2011] V. A. Kostelecký and N. Russell, Rev. Mod. Phys. 83, 11 (2011), arXiv:0801.0287 [hep-ph] .
  * Aad  _et al._ [2012] G. Aad _et al._ (ATLAS), Phys. Lett. B 716, 1 (2012), arXiv:1207.7214 [hep-ex] .
  * Chatrchyan  _et al._ [2012a] S. Chatrchyan _et al._ (CMS), Phys. Lett. B 716, 30 (2012a), arXiv:1207.7235 [hep-ex] .
  * Abazov  _et al._ [2011] V. M. Abazov _et al._ (D0), Phys. Rev. D 84, 052005 (2011), arXiv:1106.2063 [hep-ex] .
  * Aaltonen  _et al._ [2013] T. Aaltonen _et al._ (CDF), Phys. Rev. D 87, 052013 (2013), arXiv:1210.6131 [hep-ex] .
  * Aad  _et al._ [2014] G. Aad _et al._ (ATLAS), Phys. Lett. B 728, 363 (2014), arXiv:1310.6527 [hep-ex] .
  * Chatrchyan  _et al._ [2012b] S. Chatrchyan _et al._ (CMS), JHEP 06, 109 (2012b), arXiv:1204.2807 [hep-ex] .
  * Chatrchyan  _et al._ [2017a] S. Chatrchyan _et al._ (CMS), Phys. Lett. B 770, 50 (2017a), arXiv:1610.09551 [hep-ex] .
  * Cembranos  _et al._ [2008] J. A. R. Cembranos, A. Rajaraman, and F. Takayama, EPL 82, 21001 (2008), arXiv:hep-ph/0609244 .
  * Abachi  _et al._ [1997] S. Abachi _et al._ (D0), Phys. Rev. Lett. 79, 1197 (1997), arXiv:hep-ex/9703008 .
  * Abe  _et al._ [1998] F. Abe _et al._ (CDF), Phys. Rev. Lett. 80, 2767 (1998), arXiv:hep-ex/9801014 .
  * Berger  _et al._ [2016] M. S. Berger, V. A. Kostelecký, and Z. Liu, Phys. Rev. D 93, 036005 (2016), arXiv:1509.08929 [hep-ph] .
  * Greenberg [2002] O. W. Greenberg, Phys. Rev. Lett. 89, 231602 (2002), arXiv:hep-ph/0201258 .
  * Kostelecký [1998] V. A. Kostelecký, Phys. Rev. Lett. 80, 1818 (1998), arXiv:hep-ph/9809572 .
  * Bluhm  _et al._ [2002] R. Bluhm, V. A. Kostelecký, C. D. Lane, and N. Russell, Phys. Rev. Lett. 88, 090801 (2002), arXiv:hep-ph/0111141 .
  * Bluhm  _et al._ [2003] R. Bluhm, V. A. Kostelecký, C. D. Lane, and N. Russell, Phys. Rev. D 68, 125008 (2003), arXiv:hep-ph/0306190 .
  * Kostelecký [2004] V. A. Kostelecký, Phys. Rev. D 69, 105009 (2004), arXiv:hep-th/0312310 .
  * Kostelecký and Mewes [2002] V. A. Kostelecký and M. Mewes, Phys. Rev. D 66, 056005 (2002), arXiv:hep-ph/0205211 .
  * Aad  _et al._ [2013] G. Aad _et al._ (ATLAS), Eur. Phys. J. C 73, 2518 (2013), arXiv:1302.4393 [hep-ex] .
  * Abdallah  _et al._ [2008] J. Abdallah _et al._ (DELPHI), Eur. Phys. J. C 55, 1 (2008), arXiv:0803.2534 [hep-ex] .
  * Chatrchyan  _et al._ [2013] S. Chatrchyan _et al._ (CMS), CMS-PAS-LUM-13-001 (2013).
  * Chatrchyan  _et al._ [2017b] S. Chatrchyan _et al._ (CMS), Phys. Lett. B 770, 50 (2017b), arXiv:1610.09551 [hep-ex] .
  * Belyaev  _et al._ [2013] A. Belyaev, N. D. Christensen, and A. Pukhov, Comput. Phys. Commun. 184, 1729 (2013), arXiv:1207.6082 [hep-ph] .
  * Alwall  _et al._ [2011] J. Alwall, M. Herquet, F. Maltoni, O. Mattelaer, and T. Stelzer, JHEP 06, 128 (2011), arXiv:1106.0522 [hep-ph] .
  * Sjöstrand  _et al._ [2015] T. Sjöstrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten, S. Mrenna, S. Prestel, C. O. Rasmussen, and P. Z. Skands, Comput. Phys. Commun. 191, 159 (2015), arXiv:1410.3012 [hep-ph] .
  * de Favereau  _et al._ [2014] J. de Favereau, C. Delaere, P. Demin, A. Giammanco, V. Lemaître, A. Mertens, and M. Selvaggi (DELPHES 3), JHEP 02, 057 (2014), arXiv:1307.6346 [hep-ex] .