# Inverse engineering of fast state transfer among coupled oscillators

Xiao-Jing Lu  School of Science, Xuchang University, Xuchang 461000, China
Departamento de Química Física, University of the Basque Country UPV/EHU,
Apdo. 644, 48080 Bilbao, Spain  Ion Lizuain  Department of Applied
Mathematics, University of the Basque Country UPV/EHU, Donostia-San Sebastián,
Spain  EHU Quantum Center, University of the Basque Country UPV/EHU  J. G.
Muga  Departamento de Química Física, University of the Basque Country
UPV/EHU, Apdo. 644, 48080 Bilbao, Spain  EHU Quantum Center, University of the
Basque Country UPV/EHU

###### Abstract

We design faster-than-adiabatic state transfers (switching of quantum numbers)
in time-dependent coupled-oscillator Hamiltonians. The manipulation to drive
the process is found using a two-dimensional invariant recently proposed in S.
Simsek and F. Mintert, Quantum 5 (2021) 409, and involves both rotation and
transient scaling of the principal axes of the potential in a Cartesian
representation. Importantly, this invariant is degenerate except for the
subspace spanned by its ground state. Such degeneracy, in general, allows for
infidelities of the final states with respect to ideal target eigenstates.
However, the value of a single control parameter can be chosen so that the
state switching is perfect for arbitrary (not necessarily known) initial
eigenstates. Additional 2D linear invariants are used to find easily the
parameter values needed and to provide generic expressions for the final
states and final energies. In particular we find time-dependent
transformations of a two-dimensional harmonic trap for a particle (such as an
ion or neutral atom) so that the final trap is rotated with respect to the
initial one, and eigenstates of the initial trap are converted into rotated
replicas at final time, in some chosen time and rotation angle.

##  I introduction

A goal of current Physics for fundamental science and technological
applications such as metrology, sensing, or information processing is to
achieve a better control of states and dynamics of single or few-body quantum
systems such as ions or neutral atoms. Both internal and motional states have
to be controlled. We shall focus here on setting fast control protocols for
two-dimensional (2D) systems described by Hamiltonians of the (dimensionless)
form

$H(t)=\frac{p^{2}_{x}}{2}+\frac{p^{2}_{y}}{2}+\frac{1}{2}M_{11}x^{2}+\frac{1}{2}M_{22}y^{2}+M_{12}xy,$
(1)

with $M$ real, symmetric ($M_{12}=M_{21}$), and controllable in time. More
specifically, we wish to consider a faster than adiabatic process to switch
eigenstates among the $x$ and $y$ oscillators, namely,
$|n,k{\rangle}_{i}\to|k,n{\rangle}_{f}$, up to a phase factor, (where $i$ and
$f$ stand for initial and final configurations) for arbitrary vibrational
quantum numbers $n$, $k$ and the uncoupling condition $M_{12}=0$ at boundary
times. This Hamiltonian may represent in particular two oscillators on a line
transiently coupled, or one particle in a time-dependent 2D trap, see other
possible interpretations in ref. Tobalina et al. (2020). Since the different
systems are mathematically equivalent and the single particle system is easier
to visualize we shall from now on use a language appropriate for the
manipulation of a Cartesian 2D potential that holds a single particle.
Specifically we shall look for a process where the final potential is rotated
by $\pi/2$ with respect to the initial one, see Fig. 1. This implies the
additional boundary relations $(M_{22})_{f}=(M_{11})_{i}$ and
$(M_{11})_{f}=(M_{22})_{i}$. Other natural boundary conditions such as
$(M_{11})_{f}=(M_{11})_{i}$ and $(M_{22})_{f}=(M_{22})_{i}$, or rotations by a
different rotation angle will be also discussed later.

Faster-than-adiabatic, excitationless rotations Masuda and Rice (2015); Masuda
and Nakamura (2011); Palmero et al. (2016); Lizuain et al. (2019, 2017), in
particular, form, together with transport Palmero et al. (2013, 2014); Lu et
al. (2018, 2014), expansions/compressions Chen et al. (2010); Palmero et al.
(2015) or separations/mergings Martínez-Garaot et al. (2018); Rodríguez-Prieto
et al. (2020), a basic set of motional operations. These operations are needed
to make the architectures based on moving the qubits encoded in atoms or
neutral atoms viable for practical implementations of information processing
Kielpinski et al. (2002). Fast rotations have been performed with increasing
accuracy experimentally but remain technically challenging Splatt et al.
(2009); Kaufmann et al. (2017); Urban et al. (2019); van Mourik et al. (2020).
From the theory side, they also imply some difficulties, even for the simple
setting of a rotating harmonic trap. In essence, a point transformation (i.e.,
one that does not mix positions and momenta) cannot separate the dynamics into
independent normal dynamical modes, so that simple inverse-engineering methods
known for the one-dimensional modes and the corresponding 1D invariants cannot
be used to find a fast, excitationless rotation protocol Lizuain et al.
(2017). Several wayouts have been investigated, such as purely numerical
optimizations Tobalina et al. (2021), or mode separations using more
complicated transformations mixing positions and momenta Lizuain et al.
(2019). However, numerical optimizations need specific calculations for each
particular state and process, and the approach based on non-point
transformations was limited by singularities in the protocol and conditions on
the rotation times and commensurate normal frequencies Lizuain et al. (2019).
Also, for a charged particle, if the trap manipulation is complemented with a
time-dependent magnetic field perpendicular to the rotation plane, and a time-
dependent scaling of the principal axes to compensate for non-inertial
effects, fast rotations can be done without final excitation Masuda and Rice
(2015); Masuda and Nakamura (2011); Lizuain et al. (2017). We assume here that
such compensation is not viable, either because the Hamiltonian does not
represent a charged particle or because of technical limitations. While these
mentioned wayouts might be effective solutions in some cases, the need for a
robust methodology to design smooth, fast trap and state rotations, remains.
Here “robustness” refers, in the spirit of simple shortcut-to-adiabaticity
(STA) processes Guéry-Odelin et al. (2019), to independence with respect to
the initial state. In this work we look for, and find, such a robust
methodology. We shall make use of recent work on two-dimensional invariants
Tobalina et al. (2020); Simsek and Mintert (2021a). First we adapt to the
rotation scenario a recently proposed two-dimensional quadratic invariant
commuting with initial and final trap Hamiltonians Simsek and Mintert (2021a)
to define the time-dependent protocols for the harmonic trap. Inverse
engineering will proceed from the 2D invariant to the evolution of the 2D
harmonic potential, which will include both rotations and scaling of the
instantaneous eigenfrequencies. Because of the degeneracy of this invariant,
an arbitrary initial eigenstate of the Hamiltonian leads to a combination of
states in the degenerate subspace which in general does not conserve the
initial energy at final time.111 The solutions of the Schrödinger equation for
a time dependent Hamiltonian may be written in terms of orthonormal
eigenstates $|\Phi_{{N,\alpha}};t{\rangle}$ of a Hermitian invariant, where
$N$ is the eigenvalue and $\alpha$ a label to distinguish orthogonal
eigenstates in the degenerate subspace, as Lewis and Riesenfeld (1969)
$|\psi(t){\rangle}=\sum_{N}\sum_{\alpha}a_{{}_{N\alpha}}e^{i\Theta_{{}_{\\!N\alpha}}\\!(t)}|\Phi_{N,\alpha};t{\rangle},$
where the $a_{{}_{N\alpha}}$ do not depend on time and
$\Theta_{{}_{\\!N\alpha}}\\!(t)$ is the Lewis-Riesenfeld phase, defined so
that $e^{i\Theta_{{}_{\\!N\alpha}}\\!(t)}|\Phi_{N,\alpha};t{\rangle}$ is a
solution of the Schrödinger equation.  (An exception is the non-degenerate
subspace of the ground state.)

To solve this problem we construct auxiliary linear invariants following
Tobalina et al. Tobalina et al. (2020). They will enable us to find the final
state, as well as the energy in terms of only a pair of complex-valued
classical trajectories. We can thus easily find parameters for which the
process performs as a perfect rotation for any (not necessarily known) initial
eigenstate.

All numerical calculations of the dynamics are performed in several ways in
this paper: The Wigner function provides a useful tool since it evolves
classically in a harmonic trap; we independently test the results by means of
a split-operator, fast-Fourier transform approach, and also using a finite
basis in a moving frame. Appendix A gives details on the Wigner function
method.

In Section II we introduce the notation and Hamiltonian; in Section III we
describe the two-dimensional invariant adapted from Simsek et al. Simsek and
Mintert (2021a); in Section IV, the potential is designed based on the
invariant, and in Section V the degeneracy of the invariant is discussed; In
Section VI, auxiliary linear invariants are introduced; Section VII uses the
different invariants to provide a generic expression of the final state, and
it describes how to achieve perfect state transfers solving the degeneracy
problem. The article ends with a Discussion and technical Appendices.

##  II Notation and model

As a relevant example of a physical setting that leads to the Hamiltonian (1)
we consider a two dimensional space of (laboratory) fixed axes
$\tilde{x},\tilde{y}$ and a particle of mass $m$ moving in an anisotropic
harmonic potential whose principal axes $\tilde{q}_{1},\tilde{q}_{2}$ rotate
by $\theta$ (counterclockwise polar angle measured from the $\tilde{x}$ axis
to the $\tilde{q}_{1}$ axis) and admits as well a time dependence of the
instantaneous eigenfrequencies $\tilde{\omega}_{1,2}(t)$. From now on we shall
consider dimensionless versions (without tilde) of the corresponding
coordinates, defined in terms of a unit of position
$\sqrt{\hbar/(m\tilde{\omega}_{1}(0))}$, unit of momentum $\sqrt{\hbar
m\tilde{\omega}_{1}(0)}$, unit of time $1/\tilde{\omega}_{1}(0)$, and unit of
energy $\hbar\tilde{\omega}_{1}(0)$. The coordinate and momentum
transformations between laboratory and rotating coordinates and momenta are

Figure 1: (Color online) Schematic representation of initial (long principal
axis in $x$ direction) and final (long principal axis in $y$ direction)
configurations of a 2D harmonic potential as well as an intermediate
configuration with principal axes along $q_{1}$ and $q_{2}$.

$\displaystyle\left(\begin{array}[]{cc}{q}_{1}\\\ {q}_{2}\\\
\end{array}\right)=W(t)\left(\begin{array}[]{cc}x\\\ y\\\
\end{array}\right),~{}~{}~{}~{}\left(\begin{array}[]{cc}{p}_{1}\\\ {p}_{2}\\\
\end{array}\right)=W(t)\left(\begin{array}[]{cc}p_{x}\\\ p_{y}\\\
\end{array}\right)$ (10)

with

$\displaystyle W(t)=\left(\begin{array}[]{cc}\cos\theta&\sin\theta\\\
-\sin\theta&\cos\theta\\\ \end{array}\right).$ (13)

Adding kinetic and potential terms the lab-frame Hamiltonian may simply be
written as

$\displaystyle H(t)$ $\displaystyle=$
$\displaystyle\frac{p^{2}_{x}}{2}+\frac{p^{2}_{y}}{2}+\frac{1}{2}{\omega}^{2}_{1}(t)q_{1}^{2}(x,y)$
(14) $\displaystyle+$
$\displaystyle\frac{1}{2}{\omega}^{2}_{2}(t)q_{2}^{2}(x,y)$ $\displaystyle=$
$\displaystyle\frac{p^{2}_{x}}{2}+\frac{p^{2}_{y}}{2}+\frac{1}{2}{\omega}^{2}_{1}(t)(x\cos\theta+y\sin\theta)^{2}$
$\displaystyle+$
$\displaystyle\frac{1}{2}{\omega}^{2}_{2}(t)(x\sin\theta-y\cos\theta)^{2},$

which is of the form (1) with

$\displaystyle M_{11}$ $\displaystyle=$
$\displaystyle{\omega}^{2}_{1}\cos^{2}\theta+{\omega}^{2}_{2}\sin^{2}\theta,$
(15) $\displaystyle M_{22}$ $\displaystyle=$
$\displaystyle{\omega}^{2}_{1}\sin^{2}\theta+{\omega}^{2}_{2}\cos^{2}\theta,$
(16) $\displaystyle M_{12}$ $\displaystyle=$ $\displaystyle
M_{21}=({\omega}^{2}_{1}-{\omega}^{2}_{2})\sin\theta\cos\theta,$ (17)

and

$\displaystyle\tan(2\theta)$ $\displaystyle=$
$\displaystyle\frac{M_{12}+M_{21}}{M_{11}-M_{22}},$
$\displaystyle{\omega}^{2}_{1}$ $\displaystyle=$
$\displaystyle\frac{\cos^{2}\theta}{\cos
2\theta}M_{11}-\frac{\sin^{2}\theta}{\cos 2\theta}M_{22},$
$\displaystyle{\omega}^{2}_{2}$ $\displaystyle=$
$\displaystyle-\frac{\sin^{2}\theta}{\cos
2\theta}M_{11}+\frac{\cos^{2}\theta}{\cos 2\theta}M_{22}.$ (18)

##  III The primary invariant to design the harmonic trap evolution

Figure 2: (Color online) Snapshots at equal time intervals of equipotential
lines (time increases from lighter to darker lines) of the trap for different
values of the control parameter $\lambda$. For $\lambda=0$ the trap simply
expands in $y$-direction and compresses in $x$-direction. For $\lambda\neq 0$
the trap rotates anticlockwise with some deformation of the principal axes. In
all cases $x$ and $y$ initial and final frequencies are switched. Parameters:
$t_{f}=5$, $w_{2}=5$.

Figure 3: (Color online) $\theta$ and $\dot{\theta}$ versus time for
$\lambda=1$ (purple solid line), $\lambda=5$ (black long dashed line),
$\lambda=10$ (red short dashed line), $\lambda=20$ (blue dot-dashed line) and
$\lambda=40$ (brown dotted line). The peak cut in (b) reaches $10$. $t_{f}=5$,
$w_{2}=5$.

Figure 4: (Color online) Coupling term $M_{12}(t)$ for different $\lambda$.
$t_{f}=5$, $w_{2}=5$.

Figure 5: (Color online) Trap frequencies squared ${\omega}^{2}_{1,2}$ versus
time for different $\lambda$, the dashed line is for ${\omega}_{1}^{2}$ and
the solid line for ${\omega}_{2}$. The lines do cross at $\lambda=0$ but the
crossing is avoided for $\lambda\neq 0$. $t_{f}=5$, $w_{2}=5$.

The Hamiltonian (1) can be written as a $4\times 4$ matrix

$\displaystyle
H=\frac{1}{2}\hat{X}^{T}\Omega\hat{X},~{}~{}~{}\Omega=\left(\begin{array}[]{cc}M&0\\\
0&\mathbf{1}\\\ \end{array}\right),$ (21)

where $\hat{X}^{T}=(x,y,p_{x},p_{y})$, $\mathbf{1}$ represents a $2\times 2$
unit matrix, and $M$ is a time-dependent (potential) matrix

$\displaystyle M=\left(\begin{array}[]{cc}M_{11}&M_{12}\\\ M_{21}&M_{22}\\\
\end{array}\right).$ (24)

So far the formal results are valid for classical or quantum particles. For a
quantum particle the components of $\hat{X}$ are to be interpreted as
operators with the standard position/momentum commutation relations. In the
following, we will construct the invariant to design the rotation following
ref. Simsek and Mintert (2021a), in a quantum frame. To find an invariant, a
quadratic ansatz may be assumed with the form

$I_{0}(t)=\frac{1}{2}\hat{X}^{T}{\Gamma}\hat{X},$ (25)

where the subscript $0$ will distinguish it from other invariants to be
defined later. The invariant $I_{0}$ has to satisfy the relation of all
invariants $\partial I(t)/\partial t=i[I(t),H(t)]$, which implies

$\frac{d{\Gamma}}{dt}=\Omega\mathcal{S}\Gamma-\Gamma\mathcal{S}\Omega,~{}~{}~{}\mathcal{S}=\left(\begin{array}[]{cc}0&{\bf
1}\\\ -{\bf 1}&0\\\ \end{array}\right).$ (26)

As invariants imply specific dynamics, $H(t)$ may be inverse engineered from
$I_{0}(t)$, but care must be exercised, as the generic form (25) might lead to
an $H(t)$ without the desired form (21). To impose the desired form, $\Gamma$
is designed with the parameterization Simsek and Mintert (2021a)

$\displaystyle\Gamma\\!=\Re\left(\begin{array}[]{cc}\dot{P}^{\dagger}\dot{P}&-\dot{P}^{\dagger}P\\\
-{P}^{\dagger}\dot{P}&P^{\dagger}{P}\\\ \end{array}\right),$ (29)

where $P$ is a $2\times 2$ complex matrix, $\Re$ stands for the real part of
the matrix, and the dots are time derivatives. Using Eq. (29) into Eq. (26),

$\displaystyle\left(\begin{array}[]{cc}\Re(\dot{P}^{\dagger}{\cal D}+{\cal
D}^{{\dagger}}\dot{P})&-\Re({\cal D}^{\dagger}P)\\\ -\Re(P^{\dagger}{\cal
D})&0\\\ \end{array}\right)=0,$ (32)

where ${\cal D}=\ddot{P}+PM$. Any $P$ such that ${\cal D}=0$ leads to an
invariant consistent with Eq. (21), but additional restrictions are necessary
to make $M$ real-symmetric. To that end it is useful to assume a polar
decomposition $P=UR$ with $U$ unitary and $R$ Hermitian and positive-
semidefinite. $R$ is left as a free matrix and $U$ is determined to make $M$
real-symmetric by defining $A=U^{\dagger}\dot{U}$, which gives $U$ from $A$
and the initial condition $U(0)=1$. Hermiticity of $M$ requires

$\displaystyle A=iR^{-2}+\frac{1}{2}[R^{-1},\dot{R}]+\frac{1}{2}R^{-1}{\cal
J}R^{-1},$ $\displaystyle\\{{\cal
J},R^{-2}\\}=[\dot{R},R^{-1}]+[R,R^{-2}]_{\dot{R}},$ (33)

where $[x,y]_{z}=xzy-yzx$ and $\\{x,y\\}=xy+yx$ is the anti-commutator. One
then finds the differential equation for $R$ Simsek and Mintert (2021a),

$\\{\ddot{R},R\\}+\\{R^{2},M\\}=2[\dot{R},R]_{A}-2RA^{2}R,$ (34)

which guarantees that $M$ is real-symmetric for all times if $R$ is real and
positive and $\dot{R}(0)=0$ Simsek and Mintert (2021a). The Hamiltonian is
inverse engineered by designing $R$ with proper boundary conditions, and then
using Eqs. (33) and (34). The boundary conditions are chosen so that

$[I_{0}(0),H(0)]=[I_{0}(t_{f}),H(t_{f})]=0,$ (35)

which implies

$\displaystyle\ddot{R}(0)=\dot{R}(0)=0,~{}~{}~{}\ddot{R}(t_{f})=\dot{R}(t_{f})=0,$
(36) $\displaystyle
R(0)=M(0)^{-\frac{1}{4}},~{}~{}~{}~{}R(t_{f})=M(t_{f})^{-\frac{1}{4}}.$ (37)

In simple scenarios without degeneracies of $H$ and $I$, the commutativity at
the boundaries (35) implies a one-to-one dynamical mapping from initial to
final eigenstates of the Hamiltonian, see e.g. Chen et al. (2010), carried out
by the eigenstates of the invariant. This is the basis for many inverse
engineering applications from $I(t)$ (representing the desired dynamics) to
$H(t)$ (the driving) Guéry-Odelin et al. (2019). The current $I_{0}$, however,
is a degenerate operator, see below, which makes the desired mapping possible,
as we shall see, but not necessary. A basic goal of this work is to find out
the way to impose that mapping.

##  IV Design of $H(t)$ from the invariant

We consider the following initial and final potentials

$\displaystyle V_{i}=\frac{1}{2}w^{2}_{1}x^{2}+\frac{1}{2}w^{2}_{2}y^{2},$
$\displaystyle V_{f}=\frac{1}{2}w^{2}_{2}x^{2}+\frac{1}{2}w^{2}_{1}y^{2},$
(38)

where the initial and final frequencies in $x$ and $y$ directions swap, and we
use a special notation for initial and final frequencies, namely,
$w_{1}\equiv\omega_{1}(0)=\omega_{2}(t_{f})$ and
$w_{2}=\omega_{2}(0)=\omega_{1}(t_{f})$. Let us recall that with our units
$w_{1}=1$ but in several equations we write $w_{1}$ (instead of its value 1)
to achieve a more clear physical interpretation. Our goal is to transform the
potential from $V_{i}$ to $V_{f}$ so that initial eigenstates
$|n,k{\rangle}_{i}$ become their rotated versions at final time, namely
$|k,n{\rangle}_{f}$, up to a phase factor. An adiabatic rotation of the
principal axes would accomplish such transformations but we want faster
drivings using the invariant as a guidance for inverse engineering $H(t)$.
From the boundary conditions in Eqs. (38) and (37), $R$ has to satisfy

$\displaystyle R(0)=\left(\begin{array}[]{cc}w^{-\frac{1}{2}}_{1}&0\\\
0&w^{-\frac{1}{2}}_{2}\\\
\end{array}\right)\\!,\,R(t_{f})=\left(\begin{array}[]{cc}w^{-\frac{1}{2}}_{2}&0\\\
0&w^{-\frac{1}{2}}_{1}\\\ \end{array}\right).$ (43)

(The exponents of $R$ are instead $+1/2$ in Simsek and Mintert (2021a), which
is a typo.)

Considering the conditions in Eqs. (36) and (37), $R(t)$ is interpolated with
a polynomial $p(t)$,

$R(t)=[1-p(t)]R(0)+p(t)R(t_{f})+(t/t_{f})^{3}(1-t/t_{f})^{3}R_{c},$ (44)

where $R_{c}$ is in principle an arbitrary real symmetric matrix and222This is
just a simple choice. Alternative forms may be used, for example, to optimize
robustness with respect to perturbations, or other variables such as transient
energies Espinós et al. (2022).

$p(t)=10(t/t_{f})^{3}-15(t/t_{f})^{4}+6(t/t_{f})^{5},$ (45)

so that the boundary conditions (36) are also satisfied. A purely diagonal
$R_{c}$ implies that the principal axes do not rotate. This would lead to a
mapping of the form $|n,k{\rangle}_{i}\to|n,k{\rangle}_{f}$ instead of the
desired $|n,k{\rangle}_{i}\to|k,n{\rangle}_{f}$ (up to phase factors). Thus we
choose instead

$\displaystyle R_{c}=\lambda(w_{1}w_{2})^{-1/4}\left(\begin{array}[]{cc}0&1\\\
1&0\\\ \end{array}\right)$ (48)

with $\lambda$ as a free parameter. Figure 2 shows snapshots of the trap (in
the laboratory frame) at equal time intervals between $t=0$ and $t_{f}$ for
different $\lambda$. For $\lambda=0$ there is no rotation but a pure
compression (in $x$ direction) and expansion (in $y$ direction) of the two
orthogonal oscillators with principal axes fixed in the laboratory frame. For
$\lambda\neq 0$ the principal axes rotate, more abruptly as $\lambda\to 0$,
see Fig. 3. For increasing values of $\lambda$ the rotation velocity shows
first a maximum at $t_{f}/2$ but eventually develops two maxima. The effect of
$\lambda$ on the coupling term $M_{12}$ is depicted in Fig. 4, and on
$\omega^{2}_{1,2}(t)$ in Fig. 5.

Initial state | $|0,0{\rangle}_{i}$ | $|1,1{\rangle}_{i}$ | $|1,0{\rangle}_{i}$ | $|0,1{\rangle}_{i}$ | $|2,0{\rangle}_{i}$ | $|0,2{\rangle}_{i}$ | $|3,0{\rangle}_{i}$ | $|0,3{\rangle}_{i}$ | $|2,1{\rangle}_{i}$ | $|1,2{\rangle}_{i}$ | $|5,2{\rangle}_{i}$ | $|2,5{\rangle}_{i}$
---|---|---|---|---|---|---|---|---|---|---|---|---
$E_{nk}(t=0)$ | 3 | 9 | 4 | 8 | 5 | 13 | 6 | 18 | 10 | 14 | 18 | 30
$\delta$ | 0 | 0 | 2.6134 | -2.6134 | 5.2268 | -5.2268 | 7.8402 | -7.8402 | 2.6134 | -2.6134 | 7.8402 | -7.8402
Table 1: For a given initial state $|n,k{\rangle}_{i}$, Initial energy
$E_{nk}(0)=(n+1/2)w_{1}+(k+1/2)w_{2}$ and final excess energy $\delta=\langle
H(t_{f})\rangle-E_{nk}(0)$. Parameters: $t_{f}=1$, $\lambda=20$, $w_{2}=5$.

##  V Degeneracy of $I_{0}$

The invariant may be calculated at initial and final times from the boundary
conditions. In terms of number operators for $x$ and $y$ oscillators it takes
the form

$I_{0}(t_{b})=n_{x}(t_{b})+n_{y}(t_{b})+1,$ (49)

where $t_{b}=0,t_{f}$. For simplicity it is useful to subtract the constant
and define $I_{+}=I_{0}-1$ as a pure sum of number operators at boundary
times,

$I_{+}(t_{b})=n_{x}(t_{b})+n_{y}(t_{b}).$ (50)

Being an invariant, the spectrum of $I_{+}(t)$ is constant, and degenerate at
all times, except for the ground state. The different subspaces can be
labelled by $N=0,1,2...$, where $N$ is the sum of quantum numbers, and their
dimension is $N+1$. The non-degenerate subspace $N=0$ is spanned by the ground
states of the Hamiltonian at boundary times. At any other time the state that
dynamically evolves from $|0,0{\rangle}_{i}$ is the nondegenerate ground
eigenstate of $I_{+}(t)$. In general it is not an eigenstate of $H(t)$ except
in the limit of very slow, adiabatic processes. The next subspace is $N=1$,
spanned by $|0,1{\rangle}_{i}$ and $|1,0{\rangle}_{i}$ at initial time and by
$|0,1{\rangle}_{f}$ and $|1,0{\rangle}_{f}$ at final time. The states that
evolve dynamically from $|0,1{\rangle}_{i}$ and $|1,0{\rangle}_{i}$ remain
orthogonal, are degenerate eigenstates of $I_{+}(t)$, and span at any time the
subspace $N=1$ of the invariant. However there is in general no guarantee that
the state $|0,1{\rangle}_{i}$ will become $|1,0{\rangle}_{f}$. A priori, all
we can say is that it will evolve into a linear combination of
$|0,1{\rangle}_{f}$ and $|1,0{\rangle}_{f}$ at $t_{f}$. As a further example,
$N=2$ is spanned at initial and final time by $|0,2{\rangle}_{i,f}$,
$|2,0{\rangle}_{i,f}$, and $|1,1{\rangle}_{i,f}$; and so on. Any state
$|n,k{\rangle}_{i}$ evolves in time within the invariant subspace $n+k$, which
is spanned by the states that evolve dynamically from the initial eigenstates
but, except for the ground state, the final state is in general a linear
combination of the eigenstates of $H(t_{f})$ that span the subspace, rather
than the desired $|k,n{\rangle}_{f}$. Numerical calculations, see an example
in Table I, demonstrate that the energy increment when the state starts at
$|n,k{\rangle}_{i}$ takes a remarkably simple form, namely,

${\langle}H(t_{f}){\rangle}-E_{nk}(0)=b(n-k)(w_{2}-w_{1}),$ (51)

where $E_{nk}(0)=E_{kn}(t_{f})=(n+1/2)w_{1}+(k+1/2)w_{2}$, and the
proportionality constant $b$ depends, for given $w_{1,2}$, only on the process
time $t_{f}$, and on $\lambda$, but not on $n$ or $k$. This might at first
sight seem to indicate that at least $|nn{\rangle}_{i},\,\,n>0$ becomes
$|nn{\rangle}_{f}$, but that is not the case in general as we shall see. In
fact we shall find from theory the physical meaning of the constant $b$,
which, for a given $t_{f}$, can be manipulated and made zero by choosing the
proper value of $\lambda$. We shall also see that this choice corresponds
indeed to the desired process $|n,k{\rangle}_{i}\to|k,n{\rangle}_{f}$, for all
$n,k$. To proceed further we need to introduce more invariants.

##  VI Construction of further invariants

To get a more explicit expression for the final state and final energy for the
processes designed from $I_{+}(t)$, we shall make use of further invariants.
The linear operators Castaños et al. (1994); Urzúa et al. (2019); Tobalina et
al. (2020)

$\displaystyle G(t)$ $\displaystyle=$ $\displaystyle
u_{x}(t)p_{x}-\dot{u}_{x}(t)x+u_{y}(t)p_{y}-\dot{u}_{y}(t)y$ (52)

are invariants provided

$\displaystyle\ddot{u}_{x}+M_{11}(t)u_{x}$ $\displaystyle=$ $\displaystyle-
M_{12}(t)u_{y},$ $\displaystyle\ddot{u}_{y}+M_{22}(t)u_{y}$ $\displaystyle=$
$\displaystyle-M_{21}(t)u_{x},$ (53)

which are classical (Newton) equations of motion driven by a Hamiltonian (1).
We may thus regard $u_{x}(t),u_{y}(t)$ as the Cartesian components of a
classical trajectory. Note that the $u_{x},u_{y}$ might be complex,
representing in that case two trajectories, since real and imaginary parts
evolve independently and solve the classical equations. It proves useful to
write $G(t_{b})$ at boundary times $t_{b}=0,t_{f}$ in terms of (initial and
final) creation and annihilation operators,

$\displaystyle G(t_{b})$ $\displaystyle=$
$\displaystyle\sum_{z=x,y}\frac{a^{{\dagger}}_{z}(t_{b})}{\sqrt{2}}\bigg{[}i\sqrt{\omega_{z}(t_{b})}u_{z}(t_{b})-\frac{\dot{u}_{z}(t_{b})}{\sqrt{\omega_{z}(t_{b})}}\bigg{]}$
(54) $\displaystyle-$
$\displaystyle\sum_{z=x,y}\frac{a_{z}(t_{b})}{\sqrt{2}}\bigg{[}i\sqrt{\omega_{z}(t_{b})}u_{z}(t_{b})+\frac{\dot{u}_{z}(t_{b})}{\sqrt{\omega_{z}(t_{b})}}\bigg{]},$

where
$a_{z}(t_{b})=\sqrt{\omega_{z}(t_{b})/2}z+ip_{z}/\sqrt{2\omega_{z}(t_{b})}$,
$z=x,y$. Different linear invariants may be constructed by choosing specific
boundary conditions for $u_{x,y}$ and $\dot{u}_{x,y}$. In particular, the
initial conditions Tobalina et al. (2020)

$\displaystyle u_{x}(0)$ $\displaystyle=$ $\displaystyle
i/\sqrt{2\omega_{x}(0)},~{}~{}\dot{u}_{x}(0)=-\sqrt{{\omega_{x}(0)}/{2}},$
$\displaystyle u_{y}(0)$ $\displaystyle=$ $\displaystyle
0,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\dot{u}_{y}(0)=0,$ (55)

define an invariant which is initially $G_{1}(0)=a_{x}(0)$. Instead, the
initial conditions

$\displaystyle u^{\prime}_{x}(0)$ $\displaystyle=$ $\displaystyle
0,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\dot{u}^{\prime}_{x}(0)=0,$
$\displaystyle u^{\prime}_{y}(0)$ $\displaystyle=$ $\displaystyle
i/\sqrt{2\omega_{y}(0)},~{}~{}\dot{u}^{\prime}_{y}(0)=-\sqrt{{\omega_{y}(0)}/{2}},$
(56)

define a different invariant which at time zero is $G_{2}(0)=a_{y}(0)$. Here
we use the prime ′ to distinguish the trajectories with initial conditions in
Eq. (56).

We may construct corresponding quadratic invariants as
$I_{1}(0)=G_{1}^{\dagger}(0)G_{1}(0)$ and
$I_{2}(0)=G_{2}^{\dagger}(0)G_{2}(0)$. Clearly
$I_{+}(0)=n_{x}(0)+n_{y}(0)=I_{1}(0)+I_{2}(0)$ so, for the $H(t)$ inverse
engineered from $I_{+}(t)$, we have that
$G_{1}^{\dagger}(t_{f})G_{1}(t_{f})+G_{2}^{\dagger}(t_{f})G_{2}(t_{f})=I_{+}(t_{f})=n_{x}(t_{f})+n_{y}(t_{f})$.
This implies that $G_{1}(t_{f})$ and $G_{2}(t_{f})$ have no
$a_{z}^{\dagger}(t_{f})$ components and the following relations hold for the
final values of the trajectories, see Eq. (54),

$\displaystyle\dot{u}_{x}(t_{f})=i\omega_{x}(t_{f})u_{x}(t_{f}),\,\,\dot{u}_{y}(t_{f})=i\omega_{y}(t_{f})u_{y}(t_{f}),$
$\displaystyle\dot{u}^{\prime}_{x}(t_{f})=i\omega_{x}(t_{f})u^{\prime}_{x}(t_{f}),\,\,\dot{u}^{\prime}_{y}(t_{f})=i\omega_{y}(t_{f})u^{\prime}_{y}(t_{f}).$
(57)

The corresponding $G_{1}(t_{f})$ and $G_{2}(t_{f})$ take the form

$\displaystyle G_{1}(t_{f})$ $\displaystyle=$ $\displaystyle
c_{x}a_{x}(t_{f})+c_{y}a_{y}(t_{f}),$ $\displaystyle G_{2}(t_{f})$
$\displaystyle=$ $\displaystyle
c^{\prime}_{x}a_{x}(t_{f})+c^{\prime}_{y}a_{y}(t_{f}),$ (58)

where, using Eq. (57),

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!c_{x}$ $\displaystyle=$
$\displaystyle-i\sqrt{2\omega_{x}(t_{f})}u_{x}(t_{f}),\,c_{y}=-i\sqrt{2\omega_{y}(t_{f})}u_{y}(t_{f}),$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!c^{\prime}_{x}$ $\displaystyle=$
$\displaystyle-i\sqrt{2\omega_{x}(t_{f})}u^{\prime}_{x}(t_{f}),\,c^{\prime}_{y}=-i\sqrt{2\omega_{y}(t_{f})}u^{\prime}_{y}(t_{f}),$
(59)

which, taking into account that $I_{+}(t_{f})=n_{x}(t_{f})+n_{y}(t_{f})$,
satisfy the relations

$c^{\ast}_{x}c_{x}+c^{\prime\ast}_{x}c^{\prime}_{x}=1,~{}~{}c^{\ast}_{y}c_{y}+c^{\prime\ast}_{y}c^{\prime}_{y}=1,~{}~{}c^{\ast}_{x}c_{y}+c^{\prime\ast}_{x}c^{\prime}_{y}=0.$
(60)

A property of invariants that we shall use repeatedly is that if
$|\psi(t){\rangle}$ is a solution of the Schrödinger equation, then
$I(t)|\psi(t){\rangle}$ is also a solution, as can be readily be checked from
the equation satisfied by $I(t)$. We shall also use that there is an
elementary (ground state) solution that begins at $|0,0{\rangle}_{i}$ and ends
at $|0,0{\rangle}_{f}$. If we now take
$|1,0\rangle_{i}=G^{{\dagger}}_{1}(0)|0,0\rangle_{i}$ as the initial state,
the final state will be

$\displaystyle|\psi(t_{f})\rangle$ $\displaystyle=$ $\displaystyle
G^{\dagger}_{1}(t_{f})|0,0\rangle_{f}=[c^{\ast}_{x}a^{{\dagger}}_{x}(t_{f})+c^{\ast}_{y}a^{{\dagger}}_{y}(t_{f})]|0,0\rangle_{f}$
(61) $\displaystyle=$ $\displaystyle
c^{\ast}_{x}|1,0\rangle_{f}+c^{\ast}_{y}|0,1\rangle_{f}.$

Using the normalization condition at final time, we thus find the additional
relation

$\langle\psi(t_{f})|\psi(t_{f})\rangle=c^{\ast}_{x}c_{x}+c^{\ast}_{y}c_{y}=1,$
(62)

and similarly we find, using $G_{2}$,

$c^{\prime\ast}_{x}c^{\prime}_{x}+c^{\prime\ast}_{y}c^{\prime}_{y}=1.$ (63)

From the Eq. (60), (62) and (63),

$c^{\ast}_{x}c_{x}=c^{\prime\ast}_{y}c^{\prime}_{y},~{}~{}c^{\ast}_{y}c_{y}=c^{\prime\ast}_{x}c^{\prime}_{x}.$
(64)

Thanks to the relations (60) and (64) several important results will depend
only on the parameter $b\equiv c^{\ast}_{x}c_{x}$ so we use a special notation
for it. We shall find shortly that this is exactly the parameter that appears
in the expression of the energy increment (51).

##  VII Final state

###  VII.1 Final state for a given initial eigenstate $|n,k{\rangle}_{i}$

Driven by the Hamiltonian designed using $I_{0}$, the evolution from the
ground state $|0,0\rangle_{i}$ will end at the ground state $|0,0\rangle_{f}$.
To see the fate of an arbitrary initial eigenstate of $H(0)$,
$|n,k\rangle_{i}$, we construct the invariant which is initially

$I_{n,k}(0)=\frac{[G^{{\dagger}}_{1}(0)]^{n}[G^{{\dagger}}_{2}(0)]^{k}}{\sqrt{n!k!}}=\frac{[a^{{\dagger}}_{x}(0)]^{n}[a^{{\dagger}}_{y}(0)]^{k}}{\sqrt{n!k!}},$
(65)

such that $I_{n,k}(0)|0,0\rangle_{i}=|n,k\rangle_{i}$. At final time it is

$\displaystyle\\!\\!I_{n,k}(t_{f})=\frac{1}{\sqrt{n!k!}}[G^{{\dagger}}_{1}(t_{f})]^{n}[G^{{\dagger}}_{2}(t_{f})]^{k}$
$\displaystyle\\!\\!=\frac{1}{\sqrt{n!k!}}[c^{\ast}_{x}a^{{\dagger}}_{x}(t_{f})+c^{\ast}_{y}a^{{\dagger}}_{y}(t_{f})]^{n}[c^{\prime\ast}_{x}a^{{\dagger}}_{x}(t_{f})+c^{\prime\ast}_{y}a^{{\dagger}}_{y}(t_{f})]^{k}$
$\displaystyle=\frac{1}{\sqrt{n!k!}}\sum^{n}_{i=0}C(n,i)[c^{\ast}_{x}a^{{\dagger}}_{x}(t_{f})]^{n-i}[c^{\ast}_{y}a^{{\dagger}}_{y}(t_{f})]^{i}$
$\displaystyle\times\sum^{k}_{j=0}C(k,j)[c^{\prime\ast}_{x}a^{{\dagger}}_{x}(t_{f})]^{k-j}[c^{\prime\ast}_{y}a^{{\dagger}}_{y}(t_{f})]^{j},$
(66)

where $C(n,k)\equiv\frac{n!}{k!(n-k)!}$. The final state will be

$\displaystyle\\!\\!\\!\\!\\!\\!|\psi(t_{f})\rangle=I_{n,k}(t_{f})|0,0\rangle_{f}$
$\displaystyle\\!\\!\\!\\!\\!\\!=\sum^{n}_{i=0}\\!\sum^{k}_{j=0}\\!\sqrt{\\!C(n,i)C(k,j)C(n\\!+\\!k\\!-\\!i-\\!j,k\\!-\\!j)C(i\\!+\\!j,i)}$
$\displaystyle\\!\\!\\!\\!\\!\\!\times\,(c^{\ast}_{1})^{n-i}(c^{\ast}_{2})^{i}(c^{\prime\ast}_{1})^{k-j}(c^{\prime\ast}_{2})^{j}|n+k-i-j,i+j\rangle_{f}.$
(67)

We can now compute the final energy
$\langle\psi(t_{f})|H(t_{f})|\psi(t_{f})\rangle$ for the given initial state
$|n,k\rangle_{i}$. An explicit computation for any $n,k$ pair using Eqs. (67),
(60), and (64) gives exactly Eq. (51).

###  VII.2 Perfect state transfers

In view of Eq. (51), the zeros of $b=c^{\ast}_{x}c_{x}$ (as a function of the
parameters $\lambda$, $t_{f}$, and $w_{2}$) play an important role. For $b=0$,
the initial and final energy are equal for all $n,k$. Moreover, since
$c_{x}=c_{y}^{\prime}=0$ for $b=0$, see Eq. (64) (equivalently
$u_{x}(t_{f})=0$ and $u_{y}^{\prime}(t_{f})=0$), $c_{y}$ and $c^{\prime}_{x}$
become complex numbers with unit modulus, see Eq. (60), and with phases $\phi$
and $\phi^{\prime}$ that depend on the final point of the trajectories, see
Eq. (59). The linear invariants $G_{1,2}$ at final time become

$\displaystyle G_{1}(t_{f})|_{b=0}$ $\displaystyle=$ $\displaystyle
e^{i\phi}a_{y}(t_{f}),$ $\displaystyle G_{2}(t_{f})|_{b=0}$ $\displaystyle=$
$\displaystyle e^{i\phi^{\prime}}a_{x}(t_{f}),$ (68)

so that the final state corresponding to $|n,k{\rangle}_{i}$ is, for $b=0$,

$\displaystyle|\psi(t_{f})\rangle|_{b=0}=I_{n,k}(t_{f})|0,0\rangle_{f}$
$\displaystyle=(c^{\ast}_{y})^{n}(c^{\prime\ast}_{x})^{k}|k,n\rangle_{f}=e^{-i(n\phi+k\phi^{\prime})}|k,n\rangle_{f}.$
(69)

In other words, the condition $b=0$ implies a perfect-transfer STA protocol.
Notice also the quadratic invariant evolutions when $b=0$,

$\displaystyle I_{1}(0)=n_{x}(0)\rightarrow I_{1}(t_{f})=n_{y}(t_{f}),$ (70)
$\displaystyle I_{2}(0)=n_{y}(0)\rightarrow I_{2}(t_{f})=n_{x}(t_{f}).$ (71)

Figure 6: (Color online) $b$ versus $\lambda$ (a) for different final times
and $w_{2}=5$; (b) for different $w_{2}$ and $t_{f}=5$.

As Fig. 6 demonstrates, for a given $w_{2}$ and $t_{f}$ one can always find
values of $\lambda$ that make $b=0$. Figure 7 shows the minimum value of such
$\lambda$ values as a function of $t_{f}$. We emphasize that $t_{f}$ is
arbitrary so that the process may be clearly non-adiabatic during the
transient protocol. Figure 8 gives some examples of transient excitations
consistent with a vanishing final excitation.

Figure 7: (Color online) For given $t_{f}$, $\lambda_{min}$ is the minimum
$\lambda$ to realize $b=0$.

Figure 8: (Color online) Instantaneous eigenenergies of the Hamiltonian (blue
solid lines) ${\cal
E}_{nk}(t)\equiv\omega_{1}(t)(n+1/2)+\omega_{2}(t)(k+1/2)$, and the
expectation value of the energy (red points) versus time for two different
initial states $|n,k{\rangle}_{i}$ and $w_{2}=5$. For (a) and (c): $t_{f}=5$,
$\lambda=18.81$; for (b) and (d): $t_{f}$=3, $\lambda=21.89$.

##  VIII Discussion

We have found that a combination of invariants, the 2D invariant $I_{0}$
proposed by Simsek and Mintert Simsek and Mintert (2021a), and the 2D linear
invariants used in Tobalina et al. Tobalina et al. (2020), allow for fast
inverse engineering state transfers, more specifically vibrational quantum
number swapping for any (not necessarily known) quantum numbers, in coupled
oscillators. The linear invariants and their powers serve to solve the
“degeneracy problem” of $I_{0}$, by which its commutation with the Hamiltonian
at boundary times does not guarantee perfect transfer except for the ground
state. They also provide explicit expression for the Lewis-Riesenfeld final
phase for each state, see Eq. (69). Different systems may be described as
coupled oscillators, and we have paid special attention to a particle in a
trap whose final form is rotated by $\pi/2$ with respect to the initial one.
In this case the method allows for rotations that produce final eigenstates of
the Hamiltonian rotated with respect to the initial eigenstates. The method
can be extended easily to arbitrary rotation angles, see Appendix B, as well
as to more complex operations such as a trap translation combined with its
rotation, i.e. to the process considered in Simsek and Mintert (2021a). The
methodology presented here may be extended into problems in higher dimensions
such as ion separations Simsek and Mintert (2021b). A further extension
concerns the boundary conditions for the potential. Suppose that a quantum
number exchange is desired, namely $|n,k{\rangle}_{i}\to|k,n{\rangle}_{f}$ up
to phase factors, with $M_{12}(0)=M_{12}(t_{f})=0$, as in the main text, but
such that the initial and final oscillators are equal, i. e., such that
$(M_{jj})_{f}=(M_{jj})_{i}$, $j=1,2$. The process may be designed in two STA
steps, the first one as in the main text, and the second one as an STA
expansion or compression of the oscillators Chen et al. (2010) keeping the
coupling $M_{12}$ zero. Quite different boundary conditions and Hamiltonian
parameter controllability have been set in ref. Villazon et al. (2019), where
an STA approach based on counter-diabatic driving and unitary transformations
to get rid of terms difficult to implement Ibáñez et al. (2012) was applied,
see also ref. Qi and Jing (2022). A detailed comparison will be carried out
elsewhere.

Finding perfect transfers requires fine tuning an invariant parameter
$\lambda$. This tuning is done by calculating a complex valued classical
trajectory, so that a final boundary condition $(u_{x}(t_{f})=0)$ is
satisfied. The resulting processes is defined by time dependent functions
$\theta(t)$, $\omega_{1,2}(t)$. As the process time is decreased, one of the
frequencies for the principal axes may become imaginary, which implies a
repeller in that direction. This effect is well known in fast STA expansion or
compression processes Chen et al. (2010) and may or may not represent a
problem to implement it in the laboratory depending on the physical setting.
The detailed discussion of specific implementations goes beyond the scope of
the present article, but we point out that in a trapped ion setting, the
rotation of the trap with varying frequencies and controlled angular speed are
possible with “point harmonic traps” making use of concentric rings to
implement a ponderomotive potentials and a rotating electrostatic quadrupole
potential Tobalina et al. (2021); Urban et al. (2019); Urban (2019). Different
(two particle) implementations may be based on two ions in controllable double
wells Sägesser et al. (2020).

##  IX Acknowledgments

We are very grateful to Florian Mintert, Selwyn Simsek, Sofía Martínez Garaot,
and Mikel Palmero for many discussions. This work was supported by National
Natural Science Foundation of China (Grant No. 12104390), the Natural Science
Foundation of Henan Province (Grants No. 212300410238), and by PGC2018-101355-
B-I00 (MCIU/AEI/FEDER,UE).

##  Appendix A Use of the Wigner function to calculate final energies

At initial time, the Hamiltonian is given by

$H(0)=\frac{p^{2}_{x}}{2}+\frac{p^{2}_{y}}{2}+\frac{1}{2}\omega^{2}_{1}(0)x^{2}+\frac{1}{2}\omega^{2}_{2}(0)y^{2}.$
(72)

Consider as initial states the product of eigenstates of the two harmonic
oscillators with vibrational quantum numbers $n$ and $k$,

$|\Psi(x,y;0)\rangle=|\phi^{(1)}_{n}(x;0)\rangle|\phi^{(2)}_{k}(y;0)\rangle$
(73)

The corresponding Wigner function is
$W(x,y,p_{x},p_{y};0)=W^{(1)}_{n}(x,p_{x};0)W^{(2)}_{k}(y,p_{y};0)$, where

$\displaystyle W^{(1)}_{n}(x,p_{x};0)\\!\\!$ $\displaystyle\\!=\\!$
$\displaystyle\\!\\!\\!\int\\!\frac{e^{-{ip_{x}z}}}{2\pi}\phi^{(1)}_{n}\\!\\!\left(x\\!+\\!\frac{z}{2};0\right)\\!\phi^{\ast(1)}_{n}\\!\\!\left(x\\!-\\!\frac{z}{2};0\right)\\!dz,$
$\displaystyle W^{(2)}_{k}(y,p_{y};0)\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!\int\\!\frac{e^{-{ip_{y}z}}}{2\pi}\phi^{(2)}_{k}\\!\left(y\\!+\\!\frac{z}{2};0\right)\\!\\!\phi^{\ast(2)}_{k}\\!\left(y\\!-\\!\frac{z}{2};0\right)\\!dz.$

Because the driving is harmonic we may use classical trajectories to propagate
the Wigner function, as well as Liouville’s theorem, namely,
$W(x_{f},y_{f},p_{x,f},p_{y,f};t_{f})=W(x,y,p_{x},p_{y};0)$, where
$x_{f},y_{f},p_{x,f},p_{y,f}$ is the classical trajectory with initial
conditions $x,y,p_{x},p_{y}$ at time $t=0$. Thus in the Weyl-Wigner
formulation, the final energy may be computed easily using the initial Wigner
function as weight function and classical trajectories to compute the Weyl
transform of the energy operator, $H_{W}(t_{f})$, which is simply the
classical expression of the energy, as $H(t_{f})$ does not involve any cross
terms between positions and momenta,

$\displaystyle\langle H(t_{f})\rangle$ $\displaystyle=$
$\displaystyle\int\\!\\!\\!...\\!\\!\int
dxdydp_{x}dp_{y}H_{W}(x_{f},y_{f},p_{x,f},p_{y,f};t_{f})$ (74)
$\displaystyle\times$ $\displaystyle W(x,y,p_{x},p_{y};0).$

The numerical computations proceed by discretizing the integral into a mesh of
initial phase points which serve as initial conditions to calculate classical
trajectories whose final values provide the corresponding final energies.

##  Appendix B Arbitrary final rotation angle

In the main text, we consider rotations of the trap from $0$ to $\pi/2$. In
fact the use of invariants to design fast state transfer processes may also be
applied to an arbitrary final rotation angle $\gamma$. Note that for an
arbitrary rotation angle $\gamma$ the Cartesian coordinates $x,y$ are still
useful as a fixed “laboratory frame” for calculations, but they are no more a
privileged frame to represent states and Hamiltonians, as $x$ and $y$
directions will be generally coupled in the final configuration. Instead it is
convenient to express initial and final states (similarly initial and final
number operators, or creation and annihilation operators) in terms of
uncoupled rotated coordinates $q_{1},q_{2}$ and momenta $p_{1},p_{2}$, see
Eqs. (10) to (14). In that representation perfect rotations take the form
$|n_{1},k_{2}{\rangle}_{i}\to|n_{1},k_{2}{\rangle}_{f}$, up to phase factors,
and at boundary times $t_{b}=0,t_{f}$, $H(t_{b})=H_{1}(t_{b})+H_{2}(t_{b})$,
where

$\displaystyle H_{1}(t_{b})$ $\displaystyle=$
$\displaystyle\frac{p^{2}_{1}}{2}+\frac{1}{2}w^{2}_{1}q_{1}^{2},$
$\displaystyle H_{2}(t_{b})$ $\displaystyle=$
$\displaystyle\frac{p^{2}_{2}}{2}+\frac{1}{2}w^{2}_{2}q_{2}^{2}.$ (75)

The corresponding invariant is (using the lab frame as in the main text)
$I_{0}(t_{b})=\frac{1}{2}X^{T}\Gamma(t_{b})X$, and

$\displaystyle\Gamma\\!=\\!\left(\\!\\!\begin{array}[]{cc}R^{-2}(t_{b})&0\\\
0&R^{2}(t_{b})\\\
\end{array}\\!\\!\right)\\!=\\!\left(\\!\\!\begin{array}[]{cc}M^{1/2}(t_{b})&0\\\
0&M^{-1/2}(t_{b})\\\ \end{array}\\!\\!\right)\\!,$ (80)

where now the boundary conditions for $R$ and $M$ at $t_{b}$ must be adapted
to the corresponding potential. Finally, we get

$I_{0}(t_{b})=n_{1}(t_{b})+n_{2}(t_{b})+1,$ (81)

similarly to the main text, but now the number operators correspond to the two
uncoupled oscillators in the orthogonal principal axes of the initial or final
trap. The auxiliary invariants can we worked out similarly to the main text to
achieve perfect state transfers.

## References

  * Tobalina et al. (2020) A. Tobalina, E. Torrontegui, I. Lizuain, M. Palmero, and J. G. Muga, Phys. Rev. A 102, 063112 (2020).
  * Masuda and Rice (2015) S. Masuda and S. A. Rice, The Journal of Physical Chemistry B 119, 11079 (2015).
  * Masuda and Nakamura (2011) S. Masuda and K. Nakamura, Phys. Rev. A 84, 043434 (2011).
  * Palmero et al. (2016) M. Palmero, S. Wang, D. Guéry-Odelin, J.-S. Li, and J. G. Muga, New J. Phys. 18, 043014 (2016).
  * Lizuain et al. (2019) I. Lizuain, A. Tobalina, A. Rodríguez-Prieto, and J. G. Muga, Journal of Physics A: Mathematical and Theoretical 52, 465301 (2019).
  * Lizuain et al. (2017) I. Lizuain, M. Palmero, and J. G. Muga, Phys. Rev. A 95, 022130 (2017).
  * Palmero et al. (2013) M. Palmero, E. Torrontegui, D. Guéry-Odelin, and J. G. Muga, Phys. Rev. A 88, 053423 (2013).
  * Palmero et al. (2014) M. Palmero, R. Bowler, J. P. Gaebler, D. Leibfried, and J. G. Muga, Phys. Rev. A 90, 053408 (2014).
  * Lu et al. (2018) X.-J. Lu, A. Ruschhaupt, and J. G. Muga, Phys. Rev. A 97, 053402 (2018).
  * Lu et al. (2014) X.-J. Lu, J. G. Muga, X. Chen, U. G. Poschinger, F. Schmidt-Kaler, and A. Ruschhaupt, Phys. Rev. A 89, 063414 (2014).
  * Chen et al. (2010) X. Chen, A. Ruschhaupt, S. Schmidt, A. del Campo, D. Guéry-Odelin, and J. G. Muga, Phys. Rev. Lett. 104, 063002 (2010).
  * Palmero et al. (2015) M. Palmero, S. Martínez-Garaot, U. G. Poschinger, A. Ruschhaupt, and J. G. Muga, New J. Phys. 17, 093031 (2015).
  * Martínez-Garaot et al. (2018) S. Martínez-Garaot, A. Rodríguez-Prieto, and J. G. Muga, Phys. Rev. A 98, 043622 (2018).
  * Rodríguez-Prieto et al. (2020) A. Rodríguez-Prieto, S. Martínez-Garaot, I. Lizuain, and J. G. Muga, Phys. Rev. Research 2, 023328 (2020).
  * Kielpinski et al. (2002) D. Kielpinski, C. Monroe, and D. J. Wineland, Nature 417, 709 (2002).
  * Splatt et al. (2009) F. Splatt, M. Harlander, M. Brownnutt, F. Zähringer, R. Blatt, and W. Hänsel, New J. Phys. 11, 103008 (2009).
  * Kaufmann et al. (2017) H. Kaufmann, T. Ruster, C. T. Schmiegelow, M. A. Luda, V. Kaushal, J. Schulz, D. von Lindenfels, F. Schmidt-Kaler, and U. G. Poschinger, Phys. Rev. A 95, 052319 (2017).
  * Urban et al. (2019) E. Urban, N. Glikin, S. Mouradian, K. Krimmel, B. Hemmerling, and H. Haeffner, Phys. Rev. Lett. 123, 133202 (2019).
  * van Mourik et al. (2020) M. W. van Mourik, E. A. Martinez, L. Gerster, P. Hrmo, T. Monz, P. Schindler, and R. Blatt, Phys. Rev. A 102, 022611 (2020).
  * Tobalina et al. (2021) A. Tobalina, J. G. Muga, I. Lizuain, and M. Palmero, Quantum Science and Technology 6, 045023 (2021).
  * Guéry-Odelin et al. (2019) D. Guéry-Odelin, A. Ruschhaupt, A. Kiely, E. Torrontegui, S. Martínez-Garaot, and J. G. Muga, Rev. Mod. Phys. 91, 045001 (2019).
  * Simsek and Mintert (2021a) S. Simsek and F. Mintert, Quantum 5, 409 (2021a).
  * Lewis and Riesenfeld (1969) H. R. Lewis and W. B. Riesenfeld, J. Math. Phys. 10, 1458 (1969).
  * Espinós et al. (2022) H. Espinós, J. Echanobe, X.-J. Lu, and J. G. Muga, arXiv preprint arXiv:2201.07555 (2022).
  * Castaños et al. (1994) O. Castaños, R. López-Peña, and V. I. Man’ko, Journal of Physics A: Mathematical and General 27, 1751 (1994).
  * Urzúa et al. (2019) A. R. Urzúa, I. Ramos-Prieto, M. Fernández-Guasti, and H. M. Moya-Cessa, Quantum Reports 1, 82 (2019).
  * Simsek and Mintert (2021b) S. Simsek and F. Mintert, _Quantum invariant-based control of interacting trapped ions_ (2021b), eprint arXiv 2112.13905.
  * Villazon et al. (2019) T. Villazon, A. Polkovnikov, and A. Chandran, Phys. Rev. A 100, 012126 (2019).
  * Ibáñez et al. (2012) S. Ibáñez, X. Chen, E. Torrontegui, J. G. Muga, and A. Ruschhaupt, Phys. Rev. Lett. 109, 100403 (2012).
  * Qi and Jing (2022) S.-f. Qi and J. Jing, arXiv preprint arXiv:2201.12536 (2022).
  * Urban (2019) E. Urban, Ph.D. thesis, University of California, Berkeley (2019).
  * Sägesser et al. (2020) T. Sägesser, R. Matt, R. Oswald, and J. P. Home, New J. Phys. 22, 073069 (2020).


Tracing snowlines and C/O ratio in a planet-hosting disk
ALMA molecular line observations towards the HD 169142 disk

Alice S. Booth 1 Charles J. Law2 Milou Temmink 1 Margot Leemker 1 Enrique Macías 3

Leiden Observatory, Leiden University, 2300 RA Leiden, the Netherlands Center for Astrophysics - Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138, USA  European Southern Observatory, Karl-Schwarzschild-Strasse 2, 85748 Garching bei München, Germany

<EMAIL_ADDRESS><EMAIL_ADDRESS>
Tracing snowlines and C/O ratio in a planet-hosting disk
A. S. Booth et al.

The composition of a forming planet is set by the material it accretes from its parent protoplanetary disk. Therefore, it is crucial to map the chemical make-up of the gas in disks to understand the chemical environment of planet formation. This paper presents molecular line observations taken with the Atacama Large Millimeter/submillimeter Array of the planet-hosting disk around the young star HD 169142. We detect N_2H^+, CH_3OH, [CI], DCN, CS, C^34S, ^13CS, H_2CS, H_2CO, HC_3N and c-C_3H_2 in this system for the first time. Combining these data with the recent detection of SO and previously published DCO^+ data, we estimate the location of H_2O and CO snowlines and investigate radial variations in the gas phase C/O ratio. We find that the HD 169142 disk has a relatively low N_2H^+ flux compared to the disks around Herbig stars HD 163296 and MWC 480 indicating less CO freeze-out and place the CO snowline beyond the millimetre disk at $\approx150$ au. The detection of CH_3OH from the inner disk is consistent with the H_2O snowline being located at the edge of the central dust cavity at $\approx20$ au. The radially varying CS/SO ratio across the proposed H_2O snowline location is consistent with this interpretation. Additionally, the detection of CH_3OH in such a warm disk adds to the growing evidence supporting the inheritance of complex ices in disks from the earlier, colder stages of star formation. Finally, we propose that the giant HD 169142 b located at 37 au is forming between the CO_2 and H_2O snowlines where the local elemental make of the gas is expected to have C/O$\approx$1.0.

§ INTRODUCTION

The elemental make-up of an exoplanet's atmosphere is set by the material it accretes locally from its parent protoplanetary disk [Öberg & Bergin, 2021].  With ground-based high-resolution spectroscopy and the James Webb Space Telescope (JWST), we can infer the C/H, O/H and C/O ratios in the atmospheres of giant exoplanets <cit.>. These observations now provide a direct link to unravelling the formation history of planets <cit.>. This is particularly important, as many of the gas-giant planets detected may not have formed at their current location and have undergone migration [Dawson & Johnson, 2018].
A key piece of this puzzle is to measure the elemental make-up of disks and in particular locating snowlines (CO and H_2O) at the time of planet formation.

The local gas and ice composition in a protoplanetary disk is regulated by various physical and chemical processes. The elemental make-up of the disk gas and ice varies as a function of radius and height due to the freeze-out of different molecules [Öberg et al., 2011]. The growth, settling, and inward drift of dust transports volatile elements leading to an elevated gas-phase C/O$>$1 in the outer regions of disks and enhanced gas-phase abundances of CO and H_2O within their respective snowlines [Du et al., 2015, Krijt et al., 2018, Zhang et al., 2019, Booth & Ilee, 2019, Banzatti et al., 2020]. Directly measuring the locations of snowlines is non-trivial. The most abundant and easily observable isotopologues of CO tend to be optically thick and trace gas in a layer above the disk midplane. The H_2O emission lines accessible at sub-millimetre wavelengths are weak and difficult to detect in Class II disks [Notsu et al., 2019, Bosman & Bergin, 2021].

Luckily we can use other molecules as chemical tracers of both the CO and H_2O snowlines. The cations DCO^+ and N_2H^+ have been used as tracers of the CO snowline and of the two species N_2H^+ is the more robust tracer [Qi et al., 2013, Qi et al., 2015]. N_2H^+ is destroyed by gas-phase CO and is thus abundant in the outer disk regions where CO is frozen-out [van 't Hoff et al., 2017]. Due to a similar chemical relationship, HCO^+ is chemical tracer of the H_2O snowline <cit.>. Additionally, the presence of methanol (CH_3OH) indicates warm enough temperatures for the sublimation of H_2O [Booth et al., 2021, van der Marel et al., 2021]. To complement these chemical tracers, the variations in C/O across disks due to snowlines can be traced by observing pairs of molecules, e.g, SO and CS or NO and CN [Le Gal et al., 2019, Booth et al., 2021, Le Gal et al., 2021, Leemker et al., 2023], or by detecting molecules that have enhanced formation in elevated C/O environments, e.g., C_2H and c-C_3H_2 [Miotello et al., 2019, Guzmán et al., 2021, Bosman et al., 2021, Calahan et al., 2023].

HD 169142 is a nearby (114.0 pc) F1 star hosting an almost face on (inclination of 13$^{\circ}$, position angle of 5$^{\circ}$) protoplanetary disk [Raman et al., 2006, Fedele et al., 2017, Gaia Collaboration et al., 2018, Vioque et al., 2018]. Multi-wavelength dust observations have revealed multiple dust rings in this transition disk that are likely carved by forming giant planets [Fedele et al., 2017, Pohl et al., 2017, Macías et al., 2017, Macías et al., 2019, Pérez et al., 2019]. CO gas detected within the central dust cavity points to the presence of a giant planet or low-mass companion within 10 au [Carmona, 2018, Leemker et al., 2022, Garg et al., 2022, Poblete et al., 2022]. In the outer disk, in a dust and gas gap, there is evidence supporting a giant planet located at 37 au [Gratton et al., 2019, Garg et al., 2022, Hammond et al., 2023]. Studies of the kinematics of the CO gas in the outer disk show strong deviations from Keplerian motion in the vicinity of the planet [Garg et al., 2022] and recently, Law et al., 2023 show evidence for shocks associated with this planet traced in SiS and SO. In this paper, we present Atacama Large Millimeter/submillimeter Array (ALMA) observations of the HD 169142 disk where we detect 13 different species. Using these data, we infer the elemental make-up of the gas and location of both the CO and H_2O snowlines in relation to these forming planets.

§ OBSERVATIONS

We have compiled ALMA observations towards the HD 169142 disk over Bands 6, 7 and 8. Full details, including, the observation dates, number of executions, baselines, integration times and maximum recoverable scales are listed in Table A.1.
All data sets were first calibrated by the ALMA staff and the subsequent data reduction, representative angular resolution and the molecules detected in each data set are summarised below:

Project 2012.1.00799.S (PI: Honda, M.) are Band 7 observations at $\approx$0.2" and in these data, we detect the CH_3OH $J=11_{(1,10)}-11_{(0,11)}$ and the SO $J=8_8-7_7$ where the SO was first presented in Law et al., 2023 along with CO isotopologues (^12CO and ^13CO $J=3-2$). The continuum data from this program were initially presented in Macías et al., 2019. For further details about the data reduction, please refer to Law et al., 2023.

Project 2015.1.00806.S (PI: Pineda, J.) are Band 7 observations at $\approx$0.03" and in these data, we detect CS $J=6-5$. These observations also provide long baselines of the CO isotopologues (^12CO and ^13CO $J=3-2$). For further details on the data reduction, please refer to Law et al., 2023.

Project 2016.1.00344.S (PI: Perez, S.) are Band 6 observations at $\approx$0.2". This project also has higher angular resolution data ($\approx$0.02") associated with it <cit.> but as we are searching for weak lines emission we only make use of the shorter baseline data, which are optimal for our science goals.
These observations consist of 5 spectral windows with 3 dedicated to the CO lines (^12CO, ^13CO, C^18O $J=2-1$) and two continuum windows centred at 218.004 GHz and 232.004 GHz. We obtained the calibrated measurement sets from the ALMA archive and performed one round of phase-only self-calibration (interval 60 seconds) on the data after flagging the line containing channels. This resulted in only a 20% improvement in the continuum signal-to-nose ratio. In these data, we detect lines of H_2CO $J=3_{(0,3)}-2_{(0,2)}$ and $J=3_{(2,2)}-3_{(2,1)}$, DCN $J=3-2$, HC_3N $J=24-23$ and c-C_3H_2 $J=6_{(0,6)}-5_{(1,5)}$ and $J=6_{(1,6)}-5_{(0,5)}$ in the 218.004 GHz continuum spectral window (see Figure C1).

Project 2016.1.00346.S (PI: Tsukagoshi, T.) are Band 8 observations at $\approx$1.0" and we detect CS $J=10-9$ and [CI] $^{3}P_1-^{3}P_0$. The dataset contains four spectral windows, which are centred at 478.651 GHz, 480.568 GHz, 489.785 GHz and 492.195 GHz with the latter two the line spectral windows at $\sim0.15$ km s$^{-1}$ spectral resolution. We obtained the calibrated measurement set from the ALMA archive and then performed self-calibration on the data. This consisted of three rounds of phase-only calibration using solution lengths (60, 20 and 6 seconds). The self-calibration yielded a signal-to-noise improvement of 40% for the continuum.

Project 2018.1.01237 (PI: Macías, E.) are Band 7 observations at $\approx$0.8" resolution. These data target and detect DCO+ $J=4-3$ and N_2H+ $J=3-2$ and also detect
DCN $J=4-3$, ^13CS $J=6-5$, C^34S $J=6-5$, H_2CS $J=8_{(1,7)}-7_{(1,6)}$, H_2CO $J=4_{(0,4)}-3_{(0,3)}$ and,
CH_3OH $J=9_{(-1,9)}-8_{(0,8)} $ and $J=6-5$ blend of lines.
Four spectral windows were used, centred at 279.515 GHz (resolution of $\sim0.15$ km s$^{-1}$), 277.990 GHz ($\sim1.2$ km s$^{-1}$), 288.147 GHz ($\sim0.29$ km s$^{-1}$) and 289.990 GHz ($\sim17$ km s$^{-1}$).
The data were then phase-only self-calibrated using the continuum emission. For this, the line emission channels were first flagged to create a continuum-only dataset.
Three rounds of phase-only self-calibration were applied, with solution intervals going down from inf (i.e., the scan duration) to 30 s. Overall, the SNR of the continuum improved from $\sim900$ to $\sim1600$. The solutions were then applied to the unflagged data.
A weak detection of N_2H^+ in the images was confirmed via both matched filtering and spectral stacking with GoFish <cit.>. The DCN, C^34S, H_2CO and one of the CH_3OH lines are observed in a continuum spectral window at $\approx$17 km s$^{-1}$ spectral resolution (see Figure C1).

The continuum was subtracted using uvcontsub with a fit order of 1 after excluding the line containing channels. The line containing channels of the continuum spectral windows were identified in an iterative process. First, the continuum was subtracted assuming no lines were present, we made a dirty image of the line, extracted spectra, identified any lines and then made a new continuum subtracted measurement set flagging these channels. We performed all of the imaging in CASA version 6.2 [McMullin et al., 2007] and used Keplerian masks for the tCLEAN mask (generates with the code from Teague, 2020) for lines at a spectral resolution $<$2 km s$^{-1}$ and hand-drawn elliptical tCLEAN masks for other lines that were observed at a coarser spectral resolution as total line width is only a few km s$^{-1}$ across the disk. Table 1 lists all of the properties of the molecular lines detected and the resulting images. All lines in the 2016.1.00344.S data aside from HC_3N were tapered to a 0.4" beam to increase signal-to-noise.

Molecular transitions and image properties for the detected lines in the HD 169142 disk.

Species      Transition     Frequency     E$_{\mathrm{up}}$      E$_{\mathrm{Aul}}$      g$_u$         Beam     $\delta$v     rms     Int. Flux

(GHz)     (K)     (s$^{-1}$)                    (km s$^{-1}$)      (mJy beam$^{-1}$)     (mJy km s${^{-1}}$)

[C I]     $^{3}P_1-^3P_0$                492.160651       24        7.994$\times10^{-4}$     3      2.0    1.07"$\times$0.51" (81)     0.15     71.9     (7.4$\pm$0.7)$\times10^{3}$

DCO^+    $J=4-3$                         288.143858       35        5.858$\times10^{-3}$     9      0.5     0.72"$\times$0.46" (89)     0.26      2.64     970$\pm$100

DCN    $J=3-2$                           217.238539       21        4.575$\times10^{-4}$     21     2.0, taper     0.41"$\times$0.37" (86)    2.7      0.67     93.0.0$\pm$10

DCN    $J=4-3$                           289.644917       35        1.124$\times10^{-3}$     27     0.5     0.73"$\times$0.46" (89)        16.2      0.24     240.0$\pm$20

N_2H+     $J=3-2$                        279.511701       27        1.259$\times10^{-3}$     63     2.0     0.91"$\times$0.58" (85)      0.14     2.95     55.0$\pm$6

CS$\dagger$     $J=7-6$                  342.882850       66        8.368$\times10^{-4}$     15     0.5     0.29"$\times$0.26"(37)     0.85     11.7     ${\gtrsim}$7$\times10^{3}$

CS     $J=10-9$                          489.750921       129       2.489$\times10^{-3}$     21     0.5     0.94"x0.40" (83)     0.15     130.7     (4.3$\pm$0.4)$\times10^{3}$

C^34S    $J=6-5$                       289.209066       38        4.812$\times10^{-4}$     13     0.5     0.73"$\times$0.46" (89)        16.2       0.23       261.0$\pm$30

^13CS    $J=6-5$                       277.455405       47        4.399$\times10^{-4}$     26     0.5     0.73"$\times$0.46" (89)         1.10       1.1       51$\pm$5

SO    $J=8_8-7_7$                        344.310612       88        5.186$\times10^{-4}$     17     2.0      0.19"$\times$0.14" (84)        0.46     1.50      120$\pm$10

H_2CS    $J=8_{(1,7)}-7_{1,6}$           278.887661       73        3.181$\times10^{-4}$     51     0.5     0.73"$\times$0.46" (89)      1.10       1.07     140$\pm$10

H_2CO    $J=3_{0,3}-2_{0,2}$         218.222192       21        2.818$\times10^{-4}$     7      2.0, taper     0.41"$\times$0.37" (86)      2.7       0.52      480$\pm$50

H_2CO    $J=3_{2,2}-3_{2,1}$         218.475632       68        1.571$\times10^{-4}$     7      2.0, taper     0.41"$\times$0.37" (86)      2.7      0.51        55$\pm$6

H_2CO    $J=4_{0,4}-3_{0,3}$         290.623405       35        6.902$\times10^{-4}$     9      0.5     0.73"$\times$0.46" (89)        16.2       0.26     900$\pm$90

CH_3OH     $J=4_{2,2}-3_{1,2}$     218.440050               45       4.687$\times10^{-5}$     36     2.0     0.20"$\times$0.15" (88)        0.46     1.44      65$\pm$7

CH_3OH    $J=9_{(-1,9)}-8_{(0,8)}$       278.304510       110       7.685$\times10^{-5}$     76     0.5      0.73"$\times$0.46" (89)           1.10       1.05      13$\pm$1

CH_3OH    blend$^*$                      290.209700       112       7.994$\times10^{-5}$     52     0.5     0.73"$\times$0.46" (89)        16.2       0.24       104$\pm$10

CH_3OH    $J=11_{(1,10)}-11_{(0,11)}$    331.502370       169       3.929$\times10^{-4}$     92     2.0     0.20"$\times$0.15" (88)        0.46     1.44     44$\pm$10

HC_3N      $J=24-23$                     218.324723       131       8.261$\times10^{-4}$     49     2.0     0.26"$\times$0.21" (74)    2.7      0.45        171$\pm$20

c-C_3H_2     $J=6_{(0,6)}-5_{(1,5)}$     217.822148       39        5.934$\times10^{-4}$     13     2.0, taper     0.41"$\times$0.37" (86)    2.7      0.59       59$\pm$6$\dagger$$\dagger$

c-C_3H_2     $J=6_{(1,6)}-5_{(0,5)}$     217.822148       39        5.934$\times10^{-4}$     39     2.0, taper     0.41"$\times$0.37" (86)    2.7      0.59       59$\pm$6$\dagger$$\dagger$

All of the molecular data are taken from the Cologne Database for Molecular Spectroscopy <cit.>. Since we are comparing molecular lines from a variety of different sets of observations we take the errors on the integrated fluxes to be the 10% absolute flux calibration uncertainty of ALMA. For reference, the C^18O $J=2-1$ disk integrated flux from Fedele et al., 2017 is 3.9$\pm$0.5 Jy km s$^{-1}$.
$^*$ Blended CH_3OH emission from up to 6 lines where the  E$_{\mathrm{up}}$, E$_{\mathrm{Aul}}$, and g$_u$ listed are of the strongest line $J=6_{(-3,4)}-5_{(-3,3)}$. Other lines are: $J=6_{(3,3)}-5_{(3,2)}$, $J=6_{(3,4)}-5_{(3,3)}$, $J=6_{(4,2)}-5_{(4,1)}$, $J=6_{(2,5)}-5_{(2,4)}$ and $J=6_{(3,3)}-5_{(3,2)}$.  $\dagger$ The CS $J=7-6$ image was generated with a Gaussian uv-taper of 0.25" and then further smoothed in the image plane with the CASA  task using a Gaussian kernel of 0.2". Due to
the lack of short-baseline data we consider this flux a lower limit. $\dagger \dagger$ flux is of two blended lines.

Integrated intensity maps of the millimetre dust and molecular line emission from the HD 169142 protoplanetary disk. Units of the colour bars are $\mathrm{mJy~beam^{-1}}$ for the continuum image and $\mathrm{mJy~beam^{-1}~km~s^{-1}}$ for the line images aside from ^12CO which is $\mathrm{Jy~beam^{-1}~km~s^{-1}}$ . The beam of each of the observations is shown in the bottom left-hand corner of the plots and the central dot denotes the position of the star.

Azimuthally averaged radial intensity profiles for the 0.9mm continuum emission and molecular lines from the HD169142 disk. The specific transitions are the same as in Figure 1. Vertical lines mark the peak locations of the two bright continuum rings. The bar in the top right of the images is the 0.5$\times$ the major axis of the beam as listed in Table 1. Units of the y-axis are
mJy beam$^{-1}$ for the continuum, Jy beam$^{-1}$ km s$^{-1}$ for the ^12CO and mJy beam$^{-1}$ km s$^{-1}$ for the rest of the lines.

§ RESULTS

All integrated intensity maps were generated with the same Keplerian or elliptical masks used in the imaging with no flux clipping to ensure accurate fluxes. Figure 1 presents the resulting maps for the different species detected in the HD 169142 disk. Where multiple transitions are detected, we show the "best image" and all of the detected transitions are shown in Figure D1.  For comparison, we also include the 0.9 mm continuum image and ^12CO map from Law et al., 2023. All of the newly reported molecules show varying ring-like morphologies with the SO and CH_3OH being the most compact. Figure 2 shows the azimuthally-averaged radial profiles for the same lines as presented in Figure 1. These profiles highlight a clear anti-correlation between SO and CS.

Table <ref> lists the disk-integrated fluxes for all of the detected lines.

We take the errors on these fluxes to be $\pm$10% due to the uncertainty on the absolute flux calibration of ALMA. We estimate disk-averaged column densities for all the molecules detected using these disk-integrated fluxes following the now standard methods <cit.>.
This assumes the gas is in local thermodynamic equilibrium and that the line emission is optically thin.
For most species, we only detect one transition therefore we calculate all of the column densities at 30 and 60 K for the molecules we expect to be tracing cooler gas and 100 K for the SO and CH_3OH which are emitting from the warmer inner disk. Where multiple lines are detected we take an average of the calculated column density.
We exclude the H_2CO $J=3_{(2,2)}-3_{(2,1)}$ and CH_3OH  $J=4_{(2,2)}-_{(1,2)}$ lines from this analysis as they are low signal-to-noise detections in the spectra (see Figure C1). Additionally, the H_2CO $J=3_{(2,2)}$-$3_{(2,1)}$ line does not share the same emission morphology as the other two strongly detected H_2CO lines (see Figure D1).

We also compute radial column densities from the emission profiles presented in Figure 2. This is done for SO, CS (using the optically thinner C^34S and assuming ^32S/^34S = 22, e.g.; Wilson, 1999), H_2CO and CH_3OH as we aim to investigate any radial variations in N(CS)/N(SO) and N(CH_3OH)/N(H_2CO) ratios. Since we do not have any excitation information we assume a fixed temperature of 60 K. The SO and CH_3OH are likely warmer than this and the CS and H_2CO may be cooler but given the data in hand this intermediate temperature is the simplest assumption to make. Figure 3 shows the radial column densities at the native angular resolution for each of the four molecules alongside a version where the SO and CH_3OH images are convolved to the same angular resolution of the ^34CS and H_2CO data using CASA tool imsmooth. We find that the N(CS)/N(SO) ratio is increasing as a function of radius and the N(CH_3OH)/N(H_2CO) ratio is decreasing as a function of radius (see Figure 3). These ratios in the context of the disk chemistry and forming planets in the HD 169142 disk will be discussed in Section 4.

Disk-averaged column densities assuming a circular emitting area of 2" for all species aside from CH_3OH and SO where we use 0.3". A range of T$\mathrm{_{ex}}$ are explored for all species and from DCN, DCO^+, H_2CO and CH_3OH.
Errors in column densities are propagated from the assumed 10% uncertainty in the disk-integrated fluxes.

Species          N (cm$^{-2}$)

$\mathrm{T_{ex}:}$ 30 K     60 K      100 K

[C I]     (5.1$\pm$0.5)$\times10^{16}$      (5.0$\pm$0.5)$\times10^{16}$     -

DCO^+     (3.0$\pm$0.3)$\times10^{11}$      (3.7$\pm$0.4)$\times10^{11}$     -

DCN     (2.6$\pm$0.3)$\times10^{11}$      (3.1$\pm$0.3)$\times10^{11}$     -

N_2H^+     (4.5$\pm$0.5)$\times10^{10}$     (5.7$\pm$0.6)$\times10^{10}$      -

CS     (3.8$\pm$0.4)$\times10^{13}$     (8.2$\pm$0.8)$\times10^{12}$     -

C^34S       (8.5$\pm$0.9)$\times10^{11}$     (8.9$\pm$0.9)$\times10^{11}$      -

^13CS       (2.5$\pm$0.3)$\times10^{11}$     (2.3$\pm$0.3)$\times10^{11}$      -

SO      -     -     (2.1$\pm$0.3)$\times10^{13*}$

H_2CS      (4.2$\pm$0.4)$\times10^{12}$     (3.5$\pm$0.4)$\times10^{12}$      -

H_2CO     (1.0$\pm$0.1)$\times10^{13}$     (1.6$\pm$0.2)$\times10^{13}$      -

CH_3OH      -     -     (3.3$\pm$0.3)$\times10^{14}$

HC_3N       (1.0$\pm$0.1)$\times10^{13}$     (2.2$\pm$0.2)$\times10^{12}$      -

c-C_3H_2       (5.6$\pm$0.6)$\times10^{11}$     (8.4$\pm$0.8)$\times10^{11}$      -

$^{*}$ taken from Law et al., 2023.

Radial profiles of column densities and column density ratios for SO, CS, CH_3OH and H_2CO in the HD 169142 disk. These are computed from the radial emission profiles in Figure 2 assuming a temperature of 60 K. Left: at the native resolution of the data, middle: convolved to the same beam, and right: ratios of the convolved profiles. For SO and CS ratio we show where N(CS)/N(SO)=1 with a horizontal line and for H_2CO and CH_3OH we show the upper-limit on the column density ratio where CH_3OH is not detected. The vertical dashed line highlights the radial location of the giant planet HD 169142 b.

Comparisons of line fluxes and line flux ratios for the TW Hya, HD 163296, MWC 480 and HD 169142 disks. Data are compiled from this work, Qi et al., 2015, Carney et al., 2018, Pegues et al., 2020, Loomis et al., 2020, Calahan et al., 2021, Öberg et al., 2021, Muñoz-Romero et al., 2023, Pegues et al., 2023. In the top panel line fluxes and their associated errors are scaled to 100 pc for each disk taking TW Hya at 60 pc, HD 169142 at 114 pc, MWC 480 at 160 pc, and HD 163296 at 101 pc [Gaia Collaboration et al., 2018, Vioque et al., 2018].

§ DISCUSSION

We have detected a range of molecules from the HD 169142 disk that can be used to unravel the physical and chemical conditions in the system. In this section, we discuss the locations of the CO and H_2O snowlines, the evidence for radial variations in C/O, and how this connects to the composition of the forming planet(s) in this disk.

§.§ Locating the CO and H_2O snowlines

The CO snowline in the HD 169142 disk has been previously investigated using observations of DCO^+ [Macías et al., 2017, Carney et al., 2018]. The dust temperature in the fiducial model of this disk from Carney et al., 2018 only reaches the freeze-out temperature of CO ($\approx$20 K) in the midplane of the very outer disk at $\approx$150 au and there is no complete freeze-out of CO indicating a CO snowline. In this work, we present the detection of N_2H^+ which is a more robust tracer of the CO snowline <cit.>. The N_2H^+ emission extends to $\approx$200 au. N_2H^+ does not solely trace the CO snowline in disks and can be present in the warm molecular layer of the disk where CO is photo-dissociated and N_2 is not yet photo-dissociated [van 't Hoff et al., 2017].
Comparing to the models presented in Anderson et al., 2019 and Trapman et al., 2022 the low N_2H^+ line flux observed towards HD 169142 relative to C^18O is consistent with minimal CO depletion, i.e., a gas-phase CO/H_2 abundance of $\approx10^{-4}$. Therefore, the weak N_2H^+ emission detected may be tracing the gas-phase chemistry in the upper layers of the disk atmosphere inside the CO snowline.
Figure 4 compares the C^18O $J=2-1$, N_2H^+ $J=3-2$, DCO^+ $J=3-2$, DCN $J=3-2$ and H_2CO $J=3_{(0,3)}-2_{(0,2)}$ line fluxes observed towards the HD 169142, TW Hya, MWC 480, and HD 163296 disks.  Scaled to a distance of 100 pc the N_2H^+ emission from HD 169142 is $\approx$10$\times$ weaker than in both of the Herbig disks and the T-Tauri disk TW Hya. Relative to the C^18O $J=2-1$ line flux the N_2H^+ is $\approx$5$\times$ weaker than in HD 163296 and MWC 480. The depletion of CO beyond the midplane CO snowline has been measured in both the HD 163296 and MWC 480 disks where Zhang et al., 2021 found a drop in CO abundance of at least a factor of 10 at a radius of 65 au and 100 au for each disk, respectively. In comparison, a similar drop in CO abundance is not required in the HD 169142 disk to reproduce the DCO^+ data [Carney et al., 2018].
That being said, if there is a CO snowline it is expected to be beyond the mm-dust rings and current known giant planet orbits [Macías et al., 2017, Carney et al., 2018].

We can trace the H_2O snowline indirectly under the assumption that the CH_3OH is thermally desorbed from the dust grains where $\mathrm{T_{dust}}>100$ K. The temperature of optically thick CO gas in the mm-dust cavity is also $>100$ K, which supports this hypothesis <cit.>. This is also consistent with the model from Carney et al., 2018 which places the H_2O midplane snowline at $\approx$20 au, which is at the edge of the inner dust ring (shown in Figure E1).

§.§ Volatile sulphur and C/O

In the HD 169142 disk, we have detections of SO, CS, C^34S, ^13CS and H_2CS. Using the ratios between the different CS isotopologues we can assess if the main isotopologue is optically thick. The inferred column density ratios of CS/C^34S (44$\pm$10 at 30 K, 10$\pm$5 at 60 K) and CS/^13CS (150$\pm$30 at 30 K, 36$\pm$8 at 60 K) compared to the ISM isotope ratios (^32S/^34S = 22, ^12C/^13C=69, e.g.; Wilson, 1999) indicate optically thick CS if the emitting layer of the gas is warm $\approx$60 K. If the gas is cooler (30 K) then the isotopes appear to be less abundant than in the ISM, but they may emit from different layers of the disk and therefore at different temperatures. Observation of CS $J=6-5$ rather than $J=10-9$ to match the C^34S and ^13CS lines, respectively, will help break this degeneracy.
The CS is ringed and the H_2CS follows this same spatial distribution with a H_2CS/CS ratio of 0.22$\pm$0.04 at 30 K and 0.18$\pm$0.04 at 60 K (assuming N(CS) =  N(C^34S)$\times$22). This value is consistent with what was observed for disks by Le Gal et al., 2019 who report that the H_2CS/CS ratio in disks is more similar to pre-stellar cores than protostellar envelopes and photo-dissociation regions where the H_2CS/CS values are higher. We also detect SO which lies interior to the CS emission. This change in volatile S carrier across the disk is likely tracing the change in C/O at the H_2O snowline. From Figure 3 the N(CS)/N(SO) ratio is $1.0\pm$0.2 at 10 au and then increases radially to $>10$ by 100 au. Since the SO data are convolved to a lower spatial resolution (0.7") the minimum N(CS)/N(SO) is likely lower than 1.0 since the emission is more compact at 0.2" and the CS already shows a gas cavity in the 0.7" data.
Compared to other Herbig disks, this ratio sits in the middle with Oph-IRS 48 having a disk averaged N(CS)/N(SO) $<$0.01 and HD 142527 a disk averaged ratio $>$1 [Booth et al., 2021, Temmink et al., 2023].
An elevated C/O in the outer disk gas where the CS is present is also supported by the rings of HC_3N and c-C_3H_2. Calahan et al., 2023 show with chemical models that in a C/O$>1$ UV-driven chemistry leads to the efficient formation of HC_3N in the gas-phase.

HD 169142 0.9 mm continuum emission (grey-scale) with coloured contours showing the 60, 70, 80 and 90% integrated flux values for C^34S and SO in the left panel and HC_3N and CH_3OH in the right panel. The green dot highlights the location of the giant planet at  37 au [Gratton et al., 2019, Garg et al., 2022, Hammond et al., 2023].
The different beams for each molecular line are shown in the same colours as the contours.

§.§ Organic content

We detect a number of organic molecules in the HD 169142 disk: H_2CO, CH_3OH, HC_3N and c-C_3H_3. It is interesting to compare the distribution of H_2CO and CH_3OH as they both share a formation route on the grains [Fuchs et al., 2009].
HD 169142 hosts a warm disk with $\mathrm{T_{dust}}>20$ K and therefore limited CO freeze-out. This implies the inheritance of CH_3OH ices rather than in-situ formation [Booth et al., 2021].
The disk temperature also favours the formation of the observed H_2CO in the gas-phase rather than on the grains [Terwisscha van Scheltinga et al., 2021]. We find that the two strongest H_2CO lines we detect have central cavities with a drop in emission where the CH_3OH is detected. This difference may be a result of excitation since these lines have upper energy levels of 21 and 35 K whereas we expect the CH_3OH to be present at temperatures $>$100 K. To potentially support this, the third H_2CO line we (weakly) detect has an upper energy level of 68 K and shows compact emission like the CH_3OH (see Figure D1).
The difference between H_2CO and CH_3OH can also be due to the gas-phase formation of the H_2CO in the outer disk which does not occur for CH_3OH.
In Figure 3 we show the radial ratio of N(CH_3OH)/N(H_2CO). Within the H_2O snowline line at 10 au this ratio is 140$\pm$40 which decreases sharply to $<$3 beyond $\approx$60 au. This is similar to the trend shown for the HD 100546 disk [Booth et al., 2021] although, the peak of the ratio is lower in HD 100546. A N(CH_3OH)/N(H_2CO) ratio $>$1 is also consistent with thermal ice desorption observed in the Class II disk Oph-IRS 48, observations of Class I protostars, and cometary ice abundances [Jørgensen et al., 2016, Bockelée-Morvan & Biver, 2017, van der Marel et al., 2021].
In the colder outer disk, the data are not sensitive enough to detect the same non-thermally desorbed reservoir of CH_3OH as seen for HD 100546 and TW Hya where N(CH_3OH)/N(H_2CO)$\approx$1-2 [Carney et al., 2019, Booth et al., 2021].
The H_2CO extends out to 200 au where a component of the H_2CO could originate from non-thermal desorption of H_2CO rich ices at the edge of the millimetre dust disk.

The H_2CO, DCO^+ and DCN line fluxes for HD 169142 are consistent with the trends derived in Pegues et al., 2023, where the Herbig disks have lower line fluxes relative to C^18O (a proxy for total gas mass) than the T-Tauri disks indicating less active cold chemistry (see Figure 4).
We take the HD169142 C^18O $J=2-1$ line flux of 3.9 Jy km s$^{-1}$ from Fedele et al., 2017.
The similarities in the DCO^+, DCN and H_2CO in the three Herbig disks may indicate that CO freeze-out is not significantly impacting the primary formation/destruction routes of these species in Herbig disks.

At the edge of mm-dust disk is also where the rings of HC_3N and c-C_3H_2 are peaking. Ilee et al., 2021 report clear rings of HC_3N and c-C_3H_2 (and CH_3CN) in the AS 209, HD 163296 and MWC 480 protoplanetary disks. In comparison to HD 169142, the HC_3N and c-C_3H_2 emission rings from these sources, all fall well within the radial extent of the mm-dust disks. However, there are some similarities between these disks. For the HD 163296 disk the peak of the H_2CO emission is just beyond the peaks of the HC_3N and c-C_3H_2 rings [Law et al., 2021, Guzmán et al., 2021, Ilee et al., 2021] and for HD 169142 we find that the H_2CO is more radially extended in than the HC_3N and c-C_3H_2. Higher angular resolution H_2CO observations are required to properly determine the relationship between the dust, H_2CO and the large organics CH_3OH, HC_3N and c-C_3H_2.

§.§ Connections to planet composition

Multiple complementary studies indicate that there are at least two giant planets in the HD 169142 system: one within the central mm-dust cavity (likely a massive companion rather than a planet) and another located in the gap between mm-dust rings at 37 au - HD 169142 b [Gratton et al., 2019, Leemker et al., 2022, Garg et al., 2022, Hammond et al., 2023]. Using the molecular line observations of the HD 169142 disk we can infer the chemical make-up of the gas in the proximity of these planets.
Figure 5 shows mm-dust emission, the expected location of HD 169142 b and contours tracing the emission from ^34CS, SO, CH_3OH and HC_3N. Under the assumption that SO and CH_3OH trace gas with C/O$<$1 any planet forming within the central dust cavity will be forming within the H_2O snowline and therefore, able accrete gas with a bulk C/O$<$1. In contrast, HD 169142 b is forming beyond the water snowline and given the temperature structure of the disk (See Figure E1), beyond the CO_2 snowline too.
Since the CO snowline lies well beyond the dust disk, HD 169142 b should be located within the CO snowline. From this, the most simple interpretation is that HD 169142 b is accreting gas with an elemental C/O ratio of 1.0 [Öberg et al., 2011]. From Figure 3 the N(CS)/N(SO) ratio at the location of HD 169142 b is 2.0$\pm$0.4. Comparing to astrochemical models, this is broadly consistent with our interpretation of an overall C/O$\approx$1 in this region of
the disk [Semenov et al., 2018, Fedele & Favre, 2020].
The [CI] also traces the gas in the outer disk and has been detected in two other Herbig disks - HD 163296 and HD 142527 [Alarcón et al., 2022, Temmink et al., 2023]. Chemical models show that [CI] is expected to trace the upper layers of the disk and can be used to determine the carbon abundance in the disk atmosphere [Kama et al., 2016]. Further work to specifically model [CI] in disks will allow for both comparisons of disk C/H and C/O to giant exoplanet atmospheres. The gas composition in the vicinity of the planet may be altered if the planet is locally heating the disk, leading to a C/O$<1$ in the gas as suggested for HD 100546 [Booth et al., 2023] and tentatively traced by the SO in HD 169142 <cit.>.

§ CONCLUSION

We have analysed ALMA data of the HD 169142 disk and presented observations of 13 different species. Our main conclusions are as follows:

* We present the first detections of N_2H^+, CH_3OH, [CI], DCN, CS, C^34S, ^13CS, H_2CS, H_2CO, HC_3N and c-C_3H_2 in the HD 169142 disk.

* The line ratios and abundance ratios of DCO^+, DCN, H_2CO, CS and H_2CS detected in HD 169142 are consistent with observations of the other Herbig disks HD 163296 and MWC 480.

* The N_2H^+ line flux from the HD 169142 disk is low relative to HD 163296 and MWC 480. This comparison and the temperature structure in the disk model from Carney et al., 2018 indicates the lack of CO freeze-out in the HD 169142 disk but, if there is a CO freeze-out in this disk it is likely beyond the mm-dust disk edge ($\approx$100 au) and out at $\approx$150 au.

* The detection of CH_3OH and SO from the inner disk can be explained by the sublimation of ices putting the H_2O snowline at the edge of the inner mm-dust $\approx$20 au. This is consistent with the brightness temperature of CO exceeding 100 K in this region of the disk [Leemker et al., 2022].

* HD 169142 is now the third warm Herbig disk where the signature of CH_3OH ice sublimation has been detected. Since CH_3OH needs to form in the cold phase, this further supports the hypothesis that complex ices can survive the disk formation process.

* The variation in the radial column density ratio of CS to SO supports gas-phase C/O variations in the HD 169142 disk across the H_2O snowline.
The large organic molecules HC_3N and c-C_3H_2 may be tracing gas with a C/O$\approx$1 and appear in rings at the edge of the mm-dust disk.

* We estimate that the giant planet HD 169142 b is forming between the CO and H_2O snowlines and given the temperature structure of the disks likely between the CO_2 and CO snowlines. This is a region of the disk where C/O$\approx$1.0. However, this ratio may be lowered if the planet is heating the disk locally resulting in the sublimation of O-rich ices. Future higher angular resolution observations will be able to further characterise the gas composition in the vicinity of this planet.  If there is another giant planet in the central cavity $<$20 au this is within the H_2O snowline and therefore the planet will be accreting gas with C/O$<$1.

Many of the line detections presented in this work are serendipitous detections in continuum spectral windows in archival ALMA data. As continuum observations with ALMA reach higher sensitivities the likelihood of picking up such lines increases. We urge the community to be mindful of the frequency ranges and spectral resolution (FDM mode) used for continuum observations so that if lines are detected there is the best opportunity possible to use them for science.

We thank Karin Öberg and Ewine van Dishoeck for their useful comments on this paper. This work makes use of the following ALMA data:  2012.1.00799.S, 2015.1.00806.S, 2016.1.00344.S, 2016.1.00346.S, 2018.1.01237. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. We acknowledge assistance from Allegro, the European ALMA Regional Centre node in the Netherlands. Astrochemistry in Leiden is supported by funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 101019751 MOLDISK). This work has used the following additional software packages that have not been referred to in the main text: Astropy, IPython, Jupyter, Matplotlib and NumPy [Astropy Collaboration et al., 2022, Pérez & Granger, 2007, Kluyver et al., 2016, Hunter, 2007, Harris et al., 2020].

[Alarcón et al., 2022]
Alarcón, F., Bergin, E. A., & Teague, R. 2022, , 941, L24

[Anderson et al., 2019]
Anderson, D. E., Blake, G. A., Bergin, E. A., et al. 2019, , 881,

[Astropy Collaboration et al., 2022]
Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022,
, 935, 167

[Banzatti et al., 2020]
Banzatti, A., Pascucci, I., Bosman, A. D., et al. 2020, , 903, 124

[Bean et al., 2023]
Bean, J. L., Xue, Q., August, P. C., et al. 2023, , 618, 43

[Bockelée-Morvan & Biver, 2017]
Bockelée-Morvan, D. & Biver, N. 2017, Philosophical Transactions of
the Royal Society of London Series A, 375, 20160252

[Booth et al., 2023]
Booth, A. S., Ilee, J. D., Walsh, C., et al. 2023, , 669, A53

[Booth et al., 2021]
Booth, A. S., van der Marel, N., Leemker, M., van Dishoeck, E. F., &
Ohashi, S. 2021a, , 651, L6

[Booth et al., 2021]
Booth, A. S., Walsh, C., Terwisscha van Scheltinga, J., et al.
2021b, Nature Astronomy, 5, 684

[Booth & Ilee, 2019]
Booth, R. A. & Ilee, J. D. 2019, , 487, 3998

[Bosman et al., 2021]
Bosman, A. D., Alarcón, F., Bergin, E. A., et al. 2021, , 257,

[Bosman & Bergin, 2021]
Bosman, A. D. & Bergin, E. A. 2021, , 918, L10

[Calahan et al., 2021]
Calahan, J. K., Bergin, E., Zhang, K., et al. 2021, , 908, 8

[Calahan et al., 2023]
Calahan, J. K., Bergin, E. A., Bosman, A. D., et al. 2023, Nature
Astronomy, 7, 49

[Carmona, 2018]
Carmona, A. 2018, in Take a Closer Look, 64

[Carney et al., 2018]
Carney, M. T., Fedele, D., Hogerheijde, M. R., et al. 2018, , 614,

[Carney et al., 2019]
Carney, M. T., Hogerheijde, M. R., Guzmán, V. V., et al. 2019,
, 623, A124

[Dawson & Johnson, 2018]
Dawson, R. I. & Johnson, J. A. 2018, , 56, 175

[Du et al., 2015]
Du, F., Bergin, E. A., & Hogerheijde, M. R. 2015, , 807, L32

[Eistrup, 2023]
Eistrup, C. 2023, ACS Earth and Space Chemistry, 7, 260

[Fedele et al., 2017]
Fedele, D., Carney, M., Hogerheijde, M. R., et al. 2017, , 600, A72

[Fedele & Favre, 2020]
Fedele, D. & Favre, C. 2020, , 638, A110

[Fuchs et al., 2009]
Fuchs, G. W., Cuppen, H. M., Ioppolo, S., et al. 2009, , 505, 629

[Gaia Collaboration et al., 2018]
Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, ,
616, A1

[Garg et al., 2022]
Garg, H., Pinte, C., Hammond, I., et al. 2022, , 517, 5942

[Gratton et al., 2019]
Gratton, R., Ligi, R., Sissa, E., et al. 2019, , 623, A140

[GRAVITY Collaboration et al., 2020]
GRAVITY Collaboration, Nowak, M., Lacour, S., et al. 2020, , 633,

[Guzmán et al., 2021]
Guzmán, V. V., Bergner, J. B., Law, C. J., et al. 2021, , 257,

[Hammond et al., 2023]
Hammond, I., Christiaens, V., Price, D. J., et al. 2023, , 522,

[Harris et al., 2020]
Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585,

[Hunter, 2007]
Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90

[Ilee et al., 2021]
Ilee, J. D., Walsh, C., Booth, A. S., et al. 2021, , 257, 9

[Jørgensen et al., 2016]
Jørgensen, J. K., van der Wiel, M. H. D., Coutens, A., et al. 2016,
, 595, A117

[Jørgensen et al., 2013]
Jørgensen, J. K., Visser, R., Sakai, N., et al. 2013, , 779,

[Kama et al., 2016]
Kama, M., Bruderer, S., van Dishoeck, E. F., et al. 2016, , 592,

[Kluyver et al., 2016]
Kluyver, T., Ragan-Kelley, B., Pérez, F., et al. 2016, in Positioning and
Power in Academic Publishing: Players, Agents and Agendas, ed. F. Loizides &
B. Schmidt, IOS Press, 87 – 90

[Krijt et al., 2018]
Krijt, S., Schwarz, K. R., Bergin, E. A., & Ciesla, F. J. 2018, ,
864, 78

[Law et al., 2023]
Law, C. J., Booth, A. S., & Öberg, K. I. 2023, , 952, L19

[Law et al., 2021]
Law, C. J., Loomis, R. A., Teague, R., et al. 2021, , 257, 3

[Le Gal et al., 2019]
Le Gal, R., Brady, M. T., Öberg, K. I., Roueff, E., & Le Petit,
F. 2019a, , 886, 86

[Le Gal et al., 2019]
Le Gal, R., Öberg, K. I., Loomis, R. A., Pegues, J., & Bergner,
J. B. 2019b, , 876, 72

[Le Gal et al., 2021]
Le Gal, R., Öberg, K. I., Teague, R., et al. 2021, , 257, 12

[Leemker et al., 2022]
Leemker, M., Booth, A. S., van Dishoeck, E. F., et al. 2022, , 663,

[Leemker et al., 2023]
Leemker, M., Booth, A. S., van Dishoeck, E. F., et al. 2023, , 673,

[Leemker et al., 2021]
Leemker, M., van't Hoff, M. L. R., Trapman, L., et al. 2021, , 646,

[Loomis et al., 2018]
Loomis, R. A., Cleeves, L. I., Öberg, K. I., et al.
2018a, , 859, 131

[Loomis et al., 2020]
Loomis, R. A., Öberg, K. I., Andrews, S. M., et al. 2020, ,
893, 101

[Loomis et al., 2018]
Loomis, R. A., Öberg, K. I., Andrews, S. M., et al.
2018b, , 155, 182

[Macías et al., 2017]
Macías, E., Anglada, G., Osorio, M., et al. 2017, , 838, 97

[Macías et al., 2019]
Macías, E., Espaillat, C. C., Osorio, M., et al. 2019, , 881,

[McMullin et al., 2007]
McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K.
2007, in Astronomical Society of the Pacific Conference Series, Vol. 376,
Astronomical Data Analysis Software and Systems XVI, ed. R. A. Shaw,
F. Hill, & D. J. Bell, 127

[Miotello et al., 2019]
Miotello, A., Facchini, S., van Dishoeck, E. F., et al. 2019, ,
631, A69

[Muñoz-Romero et al., 2023]
Muñoz-Romero, C. E., Öberg, K. I., Law, C. J., et al. 2023,
, 943, 35

[Müller et al., 2005]
Müller, H. S. P., Schlöder, F., Stutzki, J., & Winnewisser, G.
2005, Journal of Molecular Structure, 742, 215

[Notsu et al., 2019]
Notsu, S., Akiyama, E., Booth, A., et al. 2019, , 875, 96

[Öberg & Bergin, 2021]
Öberg, K. I. & Bergin, E. A. 2021, , 893, 1

[Öberg et al., 2021]
Öberg, K. I., Guzmán, V. V., Walsh, C., et al. 2021, ,
257, 1

[Öberg et al., 2011]
Öberg, K. I., Murray-Clay, R., & Bergin, E. A. 2011, , 743, L16

[Pegues et al., 2020]
Pegues, J., Öberg, K. I., Bergner, J. B., et al. 2020, , 890,

[Pegues et al., 2023]
Pegues, J., Öberg, K. I., Qi, C., et al. 2023, arXiv e-prints,

[Pérez & Granger, 2007]
Pérez, F. & Granger, B. E. 2007, Computing in Science and Engineering, 9, 21

[Pérez et al., 2019]
Pérez, S., Casassus, S., Baruteau, C., et al. 2019, , 158, 15

[Poblete et al., 2022]
Poblete, P. P., Cuello, N., Pérez, S., et al. 2022, , 510,

[Pohl et al., 2017]
Pohl, A., Benisty, M., Pinilla, P., et al. 2017, , 850, 52

[Qi et al., 2015]
Qi, C., Öberg, K. I., Andrews, S. M., et al. 2015, , 813, 128

[Qi et al., 2013]
Qi, C., Öberg, K. I., Wilner, D. J., et al. 2013, Science, 341, 630

[Raman et al., 2006]
Raman, A., Lisanti, M., Wilner, D. J., Qi, C., & Hogerheijde, M.
2006, , 131, 2290

[Reggiani et al., 2022]
Reggiani, H., Schlaufman, K. C., Healy, B. F., Lothringer, J. D., &
Sing, D. K. 2022, , 163, 159

[Semenov et al., 2018]
Semenov, D., Favre, C., Fedele, D., et al. 2018, , 617, A28

[Teague, 2019]
Teague, R. 2019, The Journal of Open Source Software, 4, 1632

[Teague, 2020]
Teague, R. 2020, richteague/keplerian_mask: Initial Release

[Temmink et al., 2023]
Temmink, M., Booth, A. S., van der Marel, N., & van Dishoeck, E. F.
2023, , 675, A131

[Terwisscha van Scheltinga et al., 2021]
Terwisscha van Scheltinga, J., Hogerheijde, M. R., Cleeves, L. I.,
et al. 2021, , 906, 111

[Trapman et al., 2022]
Trapman, L., Zhang, K., van't Hoff, M. L. R., Hogerheijde, M. R., &
Bergin, E. A. 2022, , 926, L2

[van der Marel et al., 2021]
van der Marel, N., Booth, A. S., Leemker, M., van Dishoeck, E. F., &
Ohashi, S. 2021, , 651, L5

[van 't Hoff et al., 2017]
van 't Hoff, M. L. R., Walsh, C., Kama, M., Facchini, S., & van
Dishoeck, E. F. 2017, , 599, A101

[Vioque et al., 2018]
Vioque, M., Oudmaijer, R. D., Baines, D., Mendigutía, I., &
Pérez-Martínez, R. 2018, , 620, A128

[Wang et al., 2020]
Wang, J., Wang, J. J., Ma, B., et al. 2020, , 160, 150

[Wilson, 1999]
Wilson, T. L. 1999, Reports on Progress in Physics, 62, 143

[Zhang et al., 2019]
Zhang, K., Bergin, E. A., Schwarz, K., Krijt, S., & Ciesla, F. 2019,
, 883, 98

[Zhang et al., 2021]
Zhang, K., Booth, A. S., Law, C. J., et al. 2021, , 257, 5

§ ALMA OBSERVATIONS

Summary of ALMA observations used in this work.

Project ID      PI      Band     Date     No. Ant.     Baselines     Integration Time     M.R.S.

[m]     [min]     ["]
2012.1.00799.S     Honda, M.    7     2015-07-26     41     15–1574     41.4     1.3

7     2015-07-27     41     15–1574     21.2     1.4

7     2015-08-08     43     35–1574     41.4     1.3

2015.1.00806.S     Pineda, J.     7     2015-12-06     32     19–7716     25.3     0.4

2016.1.00344.S    Pérez, S.     6     2016-10-04     40      19-3144     33.3     8.5

2017-07-05     44      21-2647      33.2     7.7

2016.1.00346.S     Tsukagoshi, T.     8     2017-04-05     40     15-279     2.92     5.9

2018.1.01237.S     Macías, E.     7     2019-04-10     43     15-500     81.0     5.7

§ COMPARISONS OF MATCHED FILTER AND GOFISH DETECTIONS OF N_2H^+

Matched filter response (left) and GoFish spectra (right) for N_2H+ in the HD 169142 disk [Loomis et al., 2018, Teague, 2019]. The matched filter response is for a 0-200 au Keplerian mask.

§ SPECTRA EXTRACTED FROM CONTINUUM SPECTRAL WINDOWS

Spectra extracted from the continuum spectral windows using a 10 circular aperture for the two data sets noted in the plot titles. The molecules detected are highlighted with vertical lines at their respective frequencies as listed in Table 1. The dashed horizontal line marks the 3 $\sigma$ level where $\sigma$ is calculated as the rms in the line free channels of the spectrum. We note that to encompass the total flux of most of the lines a 20 aperture is required but the 10 aperture shows both the detection of compact CH_3OH and the more extended molecules, e.g., H_2CO.

§ SUPPLEMENTARY LINE IMAGES

Integrated intensity maps of additional molecular transitions detected in the HD 169142 disk.

§ HD 169142 MODEL DUST TEMPERATURE

Dust temperature in the fiducial DALI model from Carney et al., 2018.