# Circularly polarized RABBIT applied to a Rabi-cycling atom

Yijie Liao  School of Physics and Wuhan National Laboratory for
Optoelectronics,
Huazhong University of Science and Technology, Wuhan, 430074,China  Edvin
Olofsson  Department of Physics, Lund University, Box 118, SE-221 00 Lund,
Sweden  Jan Marcus Dahlström<EMAIL_ADDRESS>Department of
Physics, Lund University, Box 118, SE-221 00 Lund, Sweden  Liang-Wen Pi
Center for Attosecond Science and Technology,
Xi’an Institute of Optics and Precision Mechanics of the Chinese Academy of
Sciences, Xi’an 710119, China  Yueming Zhou<EMAIL_ADDRESS>School of
Physics and Wuhan National Laboratory for Optoelectronics,
Huazhong University of Science and Technology, Wuhan, 430074,China  Peixiang
Lu  School of Physics and Wuhan National Laboratory for Optoelectronics,
Huazhong University of Science and Technology, Wuhan, 430074,China  Optics
Valley Laboratory, Hubei 430074, China

###### Abstract

We utilize the reconstruction of attosecond beating by interference of two-
photon transitions (RABBIT) technique to study the phase of a Rabi-cycling
atom using circularly polarized extreme ultraviolet and infrared (IR) fields,
where the IR field induces Rabi oscillations between the 2s and 2p states of
lithium. Autler-Townes splittings are observed in sidebands of the
photoelectron spectra and the relative phases of outgoing electron wave
packets are retrieved from the azimuthal angle. In this RABBIT scheme, more
ionization pathways beyond the usual two-photon pathways are required. Our
results show that the polar-angle-integrated and polar-angle-resolved RABBIT
phases have different behaviors when the XUV and IR fields have co- and
counter-rotating circular polarizations, which are traced back to the
different ionization channels according to the selection rules in these two
cases and their competition relying on the propensity rule in laser-assisted
photoionization.

††preprint: APS/123-QED

##  I INTRODUCTION

Rabi oscillations Rabi (1937) are ubiquitous in nuclear physics Pileio  _et
al._ (2009), atomic physics Dudin  _et al._ (2012), condensed-matter physics
Batalov  _et al._ (2008) and quantum information science Stievater  _et al._
(2001). In the semiclassical interpretation, when a single atom is radiated by
a laser field that is nearly resonant with the transition from the ground
state to an excited bound state Merlin (2021), the population is periodically
transferred via absorption and stimulated emission of a photon Knight and
Milonni (1980). Furthermore, when Rabi oscillations are measured by light-
induced perturbative or non-perturbative transitions to another bound or
continuum state, Autler-Townes (AT) splittings Autler and Townes (1955) are
manifested in the energy domain Jiang  _et al._ (2021); Nandi  _et al._
(2022); Liao  _et al._ (2022); Zhang  _et al._ (2022); Olofsson and Dahlström
(2023); Chu and Lin (2013); Wu  _et al._ (2013); Ott  _et al._ (2014); Harkema
_et al._ (2021), in analogy with Ramsey interference Ramsey and Hellwarth
(1956); Xu and Heinzen (1999) from energy domain to time domain. The formation
of the AT splitting separated by the Rabi frequency $\Omega_{\rm R}$ can be
understood by the quasi-eigenenergies in Floquet theory Floquet (1883); Chu
and Telnov (2004). In the Jaynes-Cummings model, the spectral AT doublet is
observed when the atom-field interaction splits a pair of near-degenerated
“unperturbed” states Jaynes and Cummings (1963); Shore and Knight (1993).

With the advent of attosecond extreme ultraviolet (XUV) pulses Hentschel  _et
al._ (2001); Paul  _et al._ (2001), attosecond pump-probe spectroscopy enables
monitoring and controlling the electronic motion on its natural time scale
Pazourek  _et al._ (2015). For instance, the RABBIT technique is widely
employed to investigate the electronic dynamics in laser-assisted
photoionization of atoms Isinger  _et al._ (2017), molecules Beaulieu  _et
al._ (2017), solids Lucchini  _et al._ (2015) and liquids Jordan  _et al._
(2020). In conventional RABBIT measurements, an XUV attosecond pulse train
(APT) ionizes a target synchronized with a weak IR probe field Dahlström  _et
al._ (2012, 2013). Owing to the interference between the outgoing electron
wave packets, the photoelectron yield of sidebands (SBs) located between the
main peaks oscillates with the time delay between the XUV and IR fields Paul
_et al._ (2001); Toma and Muller (2002); Véniard  _et al._ (1996). From the
modulations of the SB signals, the relative phase of the ionized electron wave
packets can be retrieved, which contains two contributions, the phase of the
laser fields Paul  _et al._ (2001); Géneaux  _et al._ (2016) and the phase
related to the electronic dynamics in laser-matter interaction Dahlström  _et
al._ (2013); Gruson  _et al._ (2016); Zhong  _et al._ (2020); Biswas  _et al._
(2020); Gong  _et al._ (2022).

In the conventional RABBIT scheme with resonant two-photon ionization of
atoms, photoelectrons are ionized from the solely populated ground state and
an abrupt phase variation around the resonance is observed in SB modulations
Swoboda  _et al._ (2010); Kheifets (2021); Autuori  _et al._ (2022). Recently,
the RABBIT simulations on a lithium atom have revealed AT splittings in the SB
signals of the photoelectron spectra when the IR field induces several cycles
of Rabi oscillations in the atom Liao  _et al._ (2022). During each Rabi cycle
on femtosecond timescale, electrons are periodically transferred between the
ground (2s) and excited (2p) states of lithium Liao  _et al._ (2022).
Meanwhile, a train of attosecond XUV bursts assisted with the IR field samples
this Rabi process by simultaneously releasing the electrons from both Rabi
states every half IR period. Hence, the out-of-phase feature of the temporal
Rabi flopping is imprinted in the spectral interference pattern of the emitted
electron wave packets in this self-referenced interferometry. Most
interestingly, it is revealed that this Rabi process leads to a near $\pi$
phase difference in SB modulations between the AT doublet in both angle-
integrated and angle-resolved Heuser  _et al._ (2016) spectroscopies Liao  _et
al._ (2022). This real-time evolution of the structured Rabi wave packets
below the ionization threshold is out of reach in the previous pump-probe
schemes where only one of the two Rabi states is measured Jiang  _et al._
(2021); Chu and Lin (2013); Wu  _et al._ (2013). In addition, the influence of
the ongoing Rabi oscillations on the phase of the outgoing electron wave
packets can not be uncovered by a one-color scheme where only an asymmetric AT
doublet is observed at a resonant laser frequency Nandi  _et al._ (2022).

During Rabi oscillations, the atom is polarized by the resonant IR field and
its dipole moment oscillates at the Rabi frequency $\Omega_{\rm R}$ and beats
at the energy spacing of the two Rabi states. In the previous streaking
Baggesen and Madsen (2010); Pazourek  _et al._ (2012); Ossiander  _et al._
(2017) and RABBIT measurements Donsa  _et al._ (2020), it was revealed that a
permanent dipole moment of a polar target (either in its initial or ionic
states) leads to an additional term in the phase of the outgoing electron wave
packet. Meanwhile, the induced dipole moment related to the target’s
polarizability only distorts the shape of the modulation in the streaking
spectrum Baggesen and Madsen (2010). However, because there are no conceptual
counterparts like permanent dipole moment nor polarizability in the case of
Rabi oscillations Sussman (2011), the influence of the Rabi process on the
phase of the ejected electron wave packets is investigated by comparing the
two peaks in the AT doublet.

In this work, we revisit the RABBIT process of a Rabi-cycling lithium atom and
adopt circularly polarized XUV Kfir  _et al._ (2015) and IR fields to avoid
scanning the time delay between the two laser fields Donsa  _et al._ (2019);
Sörngård  _et al._ (2020); Han  _et al._ (2023). In this single time-delay
RABBIT measurement, the relative phases of outgoing electron wave packets are
retrieved from the photoemission anisotropy along the azimuthal direction.
This relative phase encodes the information about the temporal Rabi dynamics
in spite of the contamination related to the bound-free and free-free
transitions in RABBIT measurements Dahlström  _et al._ (2013). In accordance
with our previous work Liao  _et al._ (2022), a near $\pi$ phase difference
related to the Rabi dynamics is observed between each pair of splitting SBs in
AT doublets for both polar-angle-integrated and polar-angle-resolved
photoelectron spectroscopies. Furthermore, the ionization channels are
controlled by using co- or counter-rotating XUV and IR fields according to the
dipole selection rules, as an advantage of utilizing circularly polarized
fields. Correspondingly, the polar-angle-integrated and polar-angle-resolved
RABBIT phases exhibit different behaviors due to different competition among
partial waves. To give a quantitative interpretation of this complex RABBIT
process, more ionization pathways with important roles are considered in
perturbation theory than the previous work Liao  _et al._ (2022), which
elaborates the modification on the Rabi-related phase by the measurement
process.

This paper is structured as follows. In Sec. II, we describe the numerical
methods in our calculations, including the atomic parameters [Sec. II.1] and
the numerical details [Sec. II.2]. In Sec. III, we introduce the theoretical
model, including the perturbative treatment of the ionization amplitudes on
the top of the Rabi model [Sec. III.1], the description of the ionization
pathways [Sec. III.2], and the review of the RABBIT scheme using circularly
polarized laser fields [Sec. III.3]. In Sec. IV, we present the numerical
results, including the co-rotating case in Sec. IV.1 and the counter-rotating
case in Sec. IV.2. We finish with a summary in Sec. V. The paper is completed
by Appendix A, which details the derivation of the wavefunction for the Rabi
process; Appendix B, which gives the derivation of the Dyson series for the
RABBIT measurement on the Rabi process; Appendix C, which supplements the
formulas of the amplitudes for less contributing ionization pathways; and
Appendix D, which selects some numerical results of the dipole transition
matrix elements for readers’ reference. Atomic units are used throughout this
paper unless otherwise stated.

##  II Methods

###  II.1 Atomic parameters

In this work, we consider a lithium atom to study the Rabi-RABBIT scheme,
where the IR field is resonantly tuned to the transition between the 2s
($m=0$) ($|\psi_{2s}\rangle$) and 2p ($m=1$ or $m=-1$)($|\psi_{2p}\rangle$)
states of lithium and thus induces the Rabi oscillations between these two
states. The lithium atom has a single electron outside a closed shell and thus
the single-active-electron approximation is reasonable for describing the
photoionization of the lithium atom. In this study, we adopt the one-electron
effective potential of the lithium atom from Ref. Sarsa  _et al._ (2004) for
the time-dependent Schrödinger equation [Sec. II.2] and for perturbative
calculations [Sec. III.1]. Using this model potential, the calculated energies
of the 2s, 2p, 3s, and 3d states of the lithium atom are given in Tab. 1. Note
that in this model, the energy spacing between the 2s and 2p levels is 1.6707
eV. It deviates from the experimental value given in NIST by 0.1771 eV Kramida
_et al._ (2022).

Bound states |  Energies (eV)
---|---
2s | -5.3821
2p | -3.7114
3s | -2.1822
3d | -1.7960
Table 1: The energies of the 2s, 2p, 3s, and 3d states of the lithium atom
calculated using the effective potential from Ref. Sarsa  _et al._ (2004).

Figure 1 shows the photoionization cross sections as a function of the
photoelectron energy for the 2s and 2p ($m=0$) states of the lithium atom
ionized by a linearly polarized laser field. As shown, the photoionization
cross sections both decrease monotonously with the photoelectron energy, but
the cross-section of the 2p state decreases faster than that of the 2s state.
The photoionization cross-section of the 2p state is larger (smaller) than
that of the 2s state below (above) the photoelectron energy of around 4 eV.
For the laser parameters in our Rabi-RABBIT scheme, we will focus on the SBs
of the photoelectron energies ranging from 10 eV to 30 eV, where single-photon
ionization from the 2s state dominates over that of the 2p state by less than
an order of magnitude.

Figure 1: Photoionization cross sections from the 2s and 2p ($m=0$) states of
the lithium atom as the functions of the photoelectron energy.

###  II.2 Time-dependent Schrödinger equation

To uncover the effects of Rabi oscillations on the phase of the ejected wave
packet, we solve the time-dependent Schrödinger equation (TDSE) for the
lithium atom. The TDSE is written as

$i\frac{\partial\Psi(\textbf{r},t)}{\partial t}=[H_{0}+H^{\rm V}_{\rm
int}(t)]\Psi(\textbf{r},t),$ (1)

where the atomic Hamiltonian is $H_{0}=\textbf{p}^{2}/2+V(\textbf{r})$. Here
the canonical momentum operator of the electron is $\textbf{p}=-i\nabla$ and
$V(\textbf{r})$ is the effective one-electron potential of the lithium atom
Sarsa  _et al._ (2004). The atom-field interaction term is written in velocity
gauge as $H^{\rm V}_{\rm int}(t;\tau)=\textbf{A}(t;\tau)\cdot\textbf{p}$,
where the external laser field is described as

$\textbf{A}(t;\tau)=\textbf{A}_{\textrm{IR}}(t;\tau)+\textbf{A}_{\textrm{XUV}}(t).$
(2)

Here $\textbf{A}_{\textrm{IR}}(t;\tau)$ and $\textbf{A}_{\textrm{XUV}}(t)$ are
the vector potentials of the IR and XUV fields, respectively, and $\tau$ is
the time delay between the two fields. In our coordinates, the polarizations
of both fields are within the $zOx$ plane and they propagate in the $y$
direction. The vector potential of the time-delayed IR field is expressed as

$\displaystyle\textbf{A}_{\textrm{IR}}(t;\tau)=$
$\displaystyle\frac{A_{\omega}\Lambda^{\rm
A}_{\textrm{IR}}(t;\tau)}{\sqrt{1+\xi^{2}}}\times$ (3a)
$\displaystyle\\{\sin[\omega(t-\tau)]\hat{e}_{z}-\xi\cos[\omega(t-\tau)]\hat{e}_{x}\\},$
with the envelope function $\displaystyle\Lambda^{\rm
A}_{\textrm{IR}}(t;\tau)=$ (3b)
$\displaystyle\begin{cases}\frac{1}{2}-\frac{1}{2}\cos[\frac{\pi(t-\tau+\tau_{\textrm{IR}})}{2T}],&-\tau_{\rm
IR}<t-\tau<-\tau_{\rm IR}+2T\\\ 1,&-\tau_{\rm IR}+2T\leqslant
t-\tau\leqslant\tau_{\rm IR}-2T\\\
\frac{1}{2}-\frac{1}{2}\cos[\frac{\pi(t-\tau-\tau_{\textrm{IR}})}{2T}],&\tau_{\rm
IR}-2T<t-\tau\leqslant\tau_{\rm IR}\\\ 0,&\rm otherwise.\end{cases}$

Here $A_{\omega}$ is the amplitude of the field, and $\hat{e}_{x}$ and
$\hat{e}_{z}$ are the unit vectors along the $x$ and $z$ axes, respectively.
In our numerical simulations, the central frequency $\omega=1.6707$ eV is used
for the IR field, which corresponds to a period of $T=2\pi/\omega=2.4963$ fs.
The duration of the IR field is $2\tau_{\textrm{IR}}=160T=399.4080$ fs,
corresponding to a spectral width of 0.023 eV. The vector potential of the XUV
field is modeled as an APT Kheifets (2021)

$\displaystyle\textbf{A}_{\textrm{XUV}}(t)=\sum_{n=-120}^{120}\frac{A_{0}\Lambda^{\rm
A,(n)}_{\rm APT}\Lambda^{\rm A,(n)}_{\rm XUV}(t)}{\sqrt{1+\xi^{2}}}\times$
(4a)
$\displaystyle\\{\sin[\omega_{\textrm{XUV}}(t-nT/2)]\hat{e}_{z}-\xi\cos[\omega_{\textrm{XUV}}(t-nT/2)]\hat{e}_{x}\\},$
where the relative amplitude of the $n$th XUV pulse in the pulse train is
$\Lambda^{\rm A,(n)}_{\rm APT}=\exp\left[-2\ln
2\frac{(nT/2)^{2}}{\tau_{\textrm{APT}}^{2}}\right],$ (4b) and the envelope
function of the $n$th XUV pulse is $\Lambda^{\rm A,(n)}_{\rm
XUV}(t)=(-1)^{n}\exp\left[-2\ln
2\frac{(t-nT/2)^{2}}{\tau_{\textrm{XUV}}^{2}}\right].$ (4c)

Here $A_{0}$ is the amplitude of the field. In our simulations, the central
frequency $\omega_{\textrm{XUV}}=15\omega$ is used for the XUV field.
$\tau_{\textrm{XUV}}=0.08T=0.1997$ fs and $\tau_{\textrm{APT}}=50T=124.8150$
fs are the durations of the XUV attosecond pulses and the APT, respectively.
The spectral width of the XUV field is 0.021 eV. In our calculations, both the
IR and XUV fields are circularly polarized, with $\xi=-1$ and $\xi=+1$
corresponding to right- and left-hand circularly polarized laser fields,
respectively.

The wavefunction of the TDSE is expanded as a partial wave series

$\Psi(\textbf{r},t)=\sum_{l=0}^{l_{\rm max}}\sum_{m=-m_{\rm max}}^{m_{\rm
max}}\frac{R_{l,m}(r,t)}{r}Y_{l,m}(\theta,\varphi),$ (5)

where $R_{l,m}(r,t)$ is the radial part of the wavefunction and
$Y_{l,m}(\theta,\varphi)$ are spherical harmonics with polar angle $\theta$
and azimuthal angle $\varphi$. The angular momentum quantum number and the
magnetic quantum number are denoted as $l$ and $m$, respectively.

In our calculations, $R_{l,m}(r,t)$ is discretized by the finite-element
discrete variable representation method Rescigno and McCurdy (2000), where the
box size is $r_{\rm max}=400$ a.u.. The numerical convergence is guaranteed
with $l_{\rm max}=15$ and $m_{\rm max}=5$. The initial state of the lithium
atom is obtained by imaginary-time propagation. The time propagation of the
wavefunction $\Psi(\textbf{r},t)$ is implemented by the split-Lanczos method
Jiang and Tian (2017); Liang  _et al._ (2022) with the time step $\Delta
t=0.02$ a.u.. In each propagation step, we apply an absorbing mask function,
$F(r)=1-1/[1+e^{(200.0-r)/2.0}]$, which splits the wavefunction
$\Psi(\textbf{r},t)$ into the inner part $\Psi_{\rm
in}(\textbf{r},t)=F(r)\Psi(\textbf{r},t)$ and the outer part $\Psi_{\rm
out}(\textbf{r},t)=[1-F(r)]\Psi(\textbf{r},t)$. The inner part $\Psi_{\rm
in}(\textbf{r},t)$ is kept in the propagation governed by the full Hamiltonian
$H_{0}+H^{V}_{\rm int}(t)$ and the outer part $\Psi_{\rm out}(\textbf{r},t)$
is approximately propagated by Coulomb-Volkov propagator Arbó  _et al._
(2008). Specifically, the ionization amplitude of the photoelectrons with the
momentum k at time $t_{i}$ is obtained by projecting the outer part $\Psi_{\rm
out}(\textbf{r},t_{i})$ on the set of Volkov states as

$f(\textbf{k},t_{i})=\langle\psi_{\textbf{k}}^{\rm
V}(\textbf{r},t_{i})|\Psi_{\rm out}(\textbf{r},t_{i})\rangle.$ (6a) Then the
total ionization amplitude at the final time $t_{f}$ is expressed as
$f(\textbf{k})=\sum_{i=1}^{N_{\rm
step}}U_{\textbf{k}}(t_{i},t_{f})f(\textbf{k},t_{i}),$ (6b) where $N_{\rm
step}$ is the number of propagation steps and the time evolution factor,
$U_{\textbf{k}}(t_{i},t_{f})=e^{-i\int_{t_{i}}^{t_{f}}[\frac{\textbf{k}^{2}}{2}+\textbf{A}(\tau)\cdot\textbf{k}]d\tau},$
(6c) is expressed in terms of the Volkov phase of the photoelectrons
accumulated from $t_{i}$ to $t_{f}$. Finally, the ionization probability
distributions are obtained as $P(\textbf{k})=|f(\textbf{k})|^{2}.$ (6d)

##  III Theory

###  III.1 Perturbative treatment of the ionization amplitudes on the top of
the Rabi model

In our perturbative calculations, we employ the length gauge to describe the
atom-field interaction term,

$H^{\rm L}_{\rm int}(t;\tau)=\textbf{E}(t;\tau)\cdot\textbf{r},$ (7)

where the external electric field
$\textbf{E}(t;\tau)\coloneqq-d\textbf{A}(t;\tau)/dt$ contains the two
contributions from the IR and XUV fields, i.e.,
$\textbf{E}(t;\tau)=\textbf{E}_{\textrm{IR}}(t;\tau)+\textbf{E}_{\textrm{XUV}}(t)$.
To describe the physical process in the Rabi-RABBIT scheme, we preferentially
deal with the Rabi oscillations within the two-level Rabi subspace
$\mathcal{R}=\\{|\psi_{2s}\rangle,|\psi_{2p}\rangle\\}$, followed by a
perturbative treatment of the transitions from Rabi subspace $\mathcal{R}$ to
its orthogonal complement $\mathcal{S}$ in Hilbert space
$\mathcal{H}=\mathcal{R}\bigoplus\mathcal{S}$. In doing so, the full
Hamiltonian in $\mathcal{H}$ space is repartitioned as

$H\coloneqq H_{0}+H^{\rm L}_{\rm int}(t;\tau)=H_{\mathcal{R}}(t)+H_{\rm
int}^{\perp\mathcal{R}}(t;\tau).$ (8)

Here the unperturbed Hamiltonian $H_{\mathcal{R}}$ governs both Rabi dynamics
within $\mathcal{R}$ subspace and the field-free evolution dynamics within
$\mathcal{S}$ subspace and it is written as

$H_{\mathcal{R}}\coloneqq RH_{0}R+H_{\rm int}^{\mathcal{R}}(t;\tau)+SH_{0}S,$
(9)

where
$R\coloneqq|\psi_{2s}\rangle\langle\psi_{2s}|+|\psi_{2p}\rangle\langle\psi_{2p}|$
and $S\coloneqq 1-R$ are the two projectors corresponding to $\mathcal{R}$ and
$\mathcal{S}$ subspaces, respectively. The interaction term $H_{\rm
int}^{\mathcal{R}}(t;\tau)$ in the Rabi Hamiltonian $H_{\mathcal{R}}$ only
contains the two rotating-wave terms related to the excitation and stimulated
emission processes induced by the IR field within Rabi subspace $\mathcal{R}$
and it reads

$\displaystyle H_{\rm int}^{\mathcal{R}}(t;\tau)=$
$\displaystyle\frac{1}{2}E_{\omega}\Lambda^{\rm E}_{\rm IR}(t;\tau)\times$
(10) $\displaystyle[$ $\displaystyle
e^{-i\omega(t-\tau)}|\psi_{2p}\rangle\langle\psi_{2p}|(\hat{\epsilon}_{\rm
IR}\cdot\textbf{r})|\psi_{2s}\rangle\langle\psi_{2s}|$ $\displaystyle+$
$\displaystyle
e^{i\omega(t-\tau)}|\psi_{2s}\rangle\langle\psi_{2s}|(\hat{\epsilon}_{\rm
IR}\cdot\textbf{r})|\psi_{2p}\rangle\langle\psi_{2p}|],$

where $E_{\omega}$ and $\Lambda^{\rm E}_{\rm IR}(t;\tau)$ are the amplitude
and the envelope of the IR electric field, respectively. The polarization
vector of the IR field is denoted as $\hat{\epsilon}_{\rm IR}$. The
interaction term $H_{\rm int}^{\perp\mathcal{R}}(t;\tau)$ of the full
Hamiltonian $H$ in Eq. (8) is treated as a perturbation to the ongoing Rabi
oscillations within $\mathcal{R}$ subspace. Using the projection operators $R$
and $S$, it can be further separated as

$\displaystyle H_{\rm int}^{\perp\mathcal{R}}(t;\tau)$ $\displaystyle\coloneqq
H^{\rm L}_{\rm int}(t;\tau)-H_{\rm int}^{\mathcal{R}}(t;\tau)$ (11)
$\displaystyle=[RH^{\rm L}_{\rm int}(t;\tau)R-H_{\rm
int}^{\mathcal{R}}(t;\tau)]$ $\displaystyle+SH^{\rm L}_{\rm int}(t;\tau)S$
$\displaystyle+[RH^{\rm L}_{\rm int}(t;\tau)S+SH^{\rm L}_{\rm int}(t;\tau)R],$

where the role of $H_{\rm int}^{\perp\mathcal{R}}$ has three aspects: within
$\mathcal{R}$, $H_{\rm int}^{\perp\mathcal{R}}$ is responsible for the two
counter-rotating-wave transitions induced by the IR field, and all possible
non-resonant transitions induced by the XUV field; within $\mathcal{S}$,
$H_{\rm int}^{\perp\mathcal{R}}$ is responsible for all possible transitions
induced by both the IR and XUV fields; crossing between $\mathcal{R}$ and
$\mathcal{S}$, $H_{\rm int}^{\perp\mathcal{R}}$ is responsible for all the
possible transitions induced by both the IR and XUV fields.

According to the Rabi model Rabi (1937), the Rabi oscillations within
$\mathcal{R}$ subspace is described by the two-level wavefunction as

$|\Psi_{\mathcal{R}}(t)\rangle=C_{2s}(t;\tau)e^{-i\omega_{2s}t}|\psi_{2s}\rangle+C_{2p}(t;\tau)e^{-i\omega_{2p}t}|\psi_{2p}\rangle,$
(12)

where $H_{0}|\psi_{2s,2p}\rangle=\omega_{2s,2p}|\psi_{2s,2p}\rangle$ with
$\omega_{2s}$ and $\omega_{2p}$ the energies of the 2s and 2p states. Within
the subspace $\mathcal{R}$, the wavefunction $\Psi_{\mathcal{R}}$ obeys the
Schrödinger equation governed by $H_{\mathcal{R}}$, i.e.,
$i\frac{d}{dt}|\Psi_{\mathcal{R}}(t)\rangle=H_{\mathcal{R}}|\Psi_{\mathcal{R}}(t)\rangle$.
Symbolically, the propagator $U_{\mathcal{R}}$ corresponding to
$H_{\mathcal{R}}$ is defined as

$|\Psi_{\mathcal{R}}(t)\rangle\coloneqq
U_{\mathcal{R}}(t,t_{0})|\Psi_{\mathcal{R}}(t_{0})\rangle.$ (13)

Considering that the system is initially in the 2s state at time $t_{0}$,
i.e.,
$|\Psi_{\mathcal{R}}(t_{0})\rangle=e^{-i\omega_{2s}t_{0}}|\psi_{2s}\rangle$,
the time-dependent coefficients in Eq. (12) are solved as (see Appendix A)

$\begin{cases}&C_{2s}(t;\tau)=\cos\left[\frac{1}{2}\int_{t_{0}}^{t}\Omega_{\rm
R}(t^{\prime};\tau)dt^{\prime}\right]\\\
&C_{2p}(t;\tau)=-ie^{i\omega\tau}\sin\left[\frac{1}{2}\int_{t_{0}}^{t}\Omega_{\rm
R}(t^{\prime};\tau)dt^{\prime}\right],\end{cases}$ (14)

where $\Omega_{\textrm{R}}(t;\tau)=\Omega_{\textrm{R}}^{0}\Lambda^{\rm E}_{\rm
IR}(t;\tau)$ is the instantaneous Rabi frequency with its maximum value of
$\Omega_{\textrm{R}}^{0}=E_{\omega}\langle\psi_{2p}|(\hat{\epsilon}_{\rm
IR}\cdot\textbf{r})|\psi_{2s}\rangle$. The factor $e^{i\omega\tau}$ in Eq.
(14) comes from the time-delayed IR field. Here and hereafter the initial time
for the evolution of the system is taken from minus infinity
$t_{0}\rightarrow-\infty$.

Figure 2: The schematic of ionization pathways for $\rm P^{h}$ with the
photoelectron energy of $2q\omega-I^{2s}_{p}+\omega_{\rm R}$ ($\omega_{\rm
R}\approx\frac{\Omega_{\rm R}^{0}}{2}$). The purple and red arrows indicate
the transition by exchanging the XUV and IR photons, respectively. The
essential pathways are drawn in thick-line arrows. The less-contributed
pathways are represented by thin-line arrows.

In the Rabi-RABBIT scheme, assuming $H_{\rm int}^{\mathcal{R}}(t_{0};\tau)=0$
and $H_{\rm int}^{\rm L}(t_{0};\tau)=0$, the system is initially in the 2s
state within $\mathcal{R}$ subspace at time $t_{0}$ and thus the initial
wavefunction of Hilbert space $\mathcal{H}$ is
$\Psi(t_{0})=\Psi_{\mathcal{R}}(t_{0})=e^{-i\omega_{2s}t_{0}}|\psi_{2s}\rangle$.
Then the wavefunction $\Psi(t)$ describing the dynamics of the system at time
$t$ in $\mathcal{H}$ space can be written in the form of the Dyson series as
Nandi  _et al._ (2022) (see Appendix B)

$\displaystyle|\Psi(t)\rangle\coloneqq U(t,t_{0})|\Psi(t_{0})\rangle$ (15)
$\displaystyle=U_{\mathcal{R}}(t,t_{0})|\Psi(t_{0})\rangle$
$\displaystyle-i\int_{t_{0}}^{t}dt_{1}U(t,t_{1})H_{\rm
int}^{\perp\mathcal{R}}(t_{1})U_{\mathcal{R}}(t_{1},t_{0})|\Psi(t_{0})\rangle$
$\displaystyle=|\Psi_{\mathcal{R}}(t)\rangle-i\int_{t_{0}}^{t}dt_{1}U(t,t_{1})H_{\rm
int}^{\perp\mathcal{R}}(t_{1})|\Psi_{\mathcal{R}}(t_{1})\rangle,$

where $U$ is the propagator related to the full Hamiltonian $H$. In the last
line of Eq. (15), the first term describes the unperturbed Rabi oscillations
within the two-level subspace $\mathcal{R}$; and the second term describes the
transition from the Rabi wavefunction $\Psi_{\mathcal{R}}$ of the subspace
$\mathcal{R}$ outwards to its orthogonal complement $\mathcal{S}$ through the
interaction $H_{\rm int}^{\perp\mathcal{R}}$ followed by its full propagation
under $H$, which contains the ionization part of the wavefunction.

According to Eq. (8), the full propagator $U(t,t_{1})$ can be alternatively
written as

$U(t,t_{1})=U_{0}(t,t_{1})-i\int_{t_{1}}^{t}dt_{2}U(t,t_{2})H^{\rm L}_{\rm
int}(t_{2})U_{0}(t_{2},t_{1})$ (16)

where $U_{0}$ is the propagator related to the field-free Hamiltonian $H_{0}$.
Here $U_{0}(t,t^{\prime})$ can be expanded on the eigenstates of $H_{0}$ as
$U_{0}(t,t^{\prime})=\sum_{\nu}\kern-15.00002pt\int\,\,\,|\psi_{\nu}\rangle\langle\psi_{\nu}|\exp[-i\omega_{\nu}(t-t^{\prime})]$
with $H_{0}|\psi_{n}\rangle=\omega_{n}|\psi_{n}\rangle$. Substituting Eq. (16)
into Eq. (15), the wavefunction can be expanded iteratively to arbitrary order
in the interaction term as

$|\Psi(t)\rangle=\sum_{N=0}^{\infty}|\bar{\Psi}^{(N)}(t)\rangle,$ (17)

where the zeroth-order wave packet is the Rabi wavefunction
$|\bar{\Psi}^{(0)}(t)\rangle=|\Psi_{\mathcal{R}}(t)\rangle$ in Eq. (12) and
where the first-order wave packet is
$|\bar{\Psi}^{(1)}(t)\rangle=-i\int_{t_{0}}^{t}dt_{1}U_{0}(t,t_{1})H_{\rm
int}^{\perp\mathcal{R}}(t_{1})|\Psi_{\mathcal{R}}(t_{1})\rangle$, which
describes the transition from the Rabi wavefunction $\Psi_{\mathcal{R}}$
through the first-order interaction $H_{\rm int}^{\perp\mathcal{R}}$.
Likewise, for the wave packets to higher orders $N>1$, $|\Psi^{(N)}(t)\rangle$
originates from the $N$th-order interaction and it is derived from Eq. (15) as

$\displaystyle|\bar{\Psi}^{(N>1)}(t)\rangle=(-i)^{N}\int_{t_{0}}^{t}dt_{1}\int_{t_{1}}^{t}dt_{2}\cdots\int_{t_{N-1}}^{t}dt_{N}\times$
(18) $\displaystyle U_{0}(t,t_{N})H^{\rm L}_{\rm
int}(t_{N})U_{0}(t_{N},t_{N-1})H^{\rm L}_{\rm int}(t_{N-1})\cdots$
$\displaystyle H^{\rm L}_{\rm int}(t_{2})U_{0}(t_{2},t_{1})H_{\rm
int}^{\perp\mathcal{R}}(t_{1})|\Psi_{\mathcal{R}}(t_{1})\rangle,$

with the time order $t>t_{N}>t_{N-1}>\cdots>t_{2}>t_{1}>t_{0}$. Here in our
perturbative treatment of the ionization process on the top of the Rabi model,
we assume the ongoing Rabi oscillations with $\mathcal{R}$ are perfect: the
energy shifts of the 2s and 2p states, the population leakage to $\mathcal{S}$
subspace, and the dampening in Rabi oscillations are all negligible. Moreover,
we assume the time scale of the process we are investigating is much shorter
than the lifetime of the 2s and 2p states.

The $N$th-order transition amplitude $a_{f}^{(N)}$ is obtained by projecting
the $N$th-order wavefunction $|\Psi^{(N)}(t)\rangle$ to the final state
$|\phi_{f}\rangle$ of a specific energy $\omega_{f}$ as Joachain  _et al._
(2011)

$a_{f}^{(N)}=\lim_{t\rightarrow\infty}\langle\phi_{f}(t)|\Psi^{(N)}(t)\rangle,$
(19)

where $|\phi_{f}(t)\rangle=e^{-i\omega_{f}t}|\psi_{f}\rangle$, with
$|\psi_{f}\rangle$ satisfying
$H_{0}|\psi_{f}\rangle=\omega_{f}|\psi_{f}\rangle$. Here we take
$t\rightarrow\infty$ because the photoelectrons are measured long after the
interaction with the field is over. In the case of $\omega_{f}>0$,
$a_{f}^{(N>0)}$ represents the $N$-photon ionization amplitude from the Rabi
wavefunction $\Psi_{\mathcal{R}}$ with all possible combinations of photons in
different time orders. To further specify different ionization processes in
the total $N$th-order transition amplitude $a_{f}^{(N)}$, the interaction
terms $H^{\rm L}_{\rm int}$ and $H_{\rm int}^{\perp\mathcal{R}}$ are separated
into different laser components in Eq. (18). For example, the interaction term
is $E_{\rm IR}(t;\tau)\hat{O}_{\omega}$ for the IR photon and it is
$E_{\Omega_{2q+1}}(t)\hat{O}_{\Omega}$ for the XUV harmonics of the frequency
$\Omega_{2q+1}=(2q+1)\omega$, which contain both absorption and emission
terms. Here the transition operators are written in the length gauge as
$\hat{O}_{\omega}\coloneqq\hat{\epsilon}_{\rm IR}\cdot\textbf{r}$ and
$\hat{O}_{\Omega}\coloneqq\hat{\epsilon}_{\rm XUV}\cdot\textbf{r}$, with
$\hat{\epsilon}_{\rm XUV}$ the polarization vector of the XUV field assuming
all XUV harmonics have the same polarization. In particular, when considering
the XUV, XUV+IR, and XUV+IR+IR absorption processes, the corresponding
ionization amplitudes are derived as

$\displaystyle a_{\rm XUV}^{(1)}=$
$\displaystyle-i\int_{-\infty}^{\infty}d\omega_{\rm
R}\tilde{C}_{2s}(\omega_{\rm
R};\tau)\tilde{E}_{\Omega_{2q+1}}(\omega_{f}-\omega_{2s}-\omega_{\rm
R})\langle\psi_{f}|\hat{O}_{\Omega}|\psi_{2s}\rangle$ (20a)
$\displaystyle-i\int_{-\infty}^{\infty}d\omega_{\rm
R}\tilde{C}_{2p}(\omega_{\rm
R};\tau)\tilde{E}_{\Omega_{2q+1}}(\omega_{f}-\omega_{2p}-\omega_{\rm
R})\langle\psi_{f}|\hat{O}_{\Omega}|\psi_{2p}\rangle,$ $\displaystyle a_{\rm
XUV+IR}^{(2)}=$
$\displaystyle-\frac{i}{\sqrt{2\pi}}\int_{-\infty}^{\infty}d\omega_{\rm
R}\tilde{C}_{2s}(\omega_{\rm R};\tau)\int_{0}^{\infty}d\Omega
e^{i(\omega_{f}-\Omega-\omega_{2s}-\omega_{\rm R})\tau}\times$ (20b)
$\displaystyle\tilde{E}_{\rm IR}(\omega_{f}-\Omega-\omega_{2s}-\omega_{\rm
R})\tilde{E}_{\Omega_{2q+1}}(\Omega)\sum_{\nu_{1}}\kern-15.00002pt\int\frac{\langle\psi_{f}|\hat{O}_{\omega}|\psi_{\nu_{1}}\rangle\langle\psi_{\nu_{1}}|\hat{O}_{\Omega}|\psi_{2s}\rangle}{\omega_{2s}+\omega_{\rm
R}+\Omega-\omega_{\nu_{1}}}$
$\displaystyle-\frac{i}{\sqrt{2\pi}}\int_{-\infty}^{\infty}d\omega_{\rm
R}\tilde{C}_{2p}(\omega_{\rm R};\tau)\int_{0}^{\infty}d\Omega
e^{i(\omega_{f}-\Omega-\omega_{2p}-\omega_{\rm R})\tau}\times$
$\displaystyle\tilde{E}_{\rm IR}(\omega_{f}-\Omega-\omega_{2p}-\omega_{\rm
R})\tilde{E}_{\Omega_{2q+1}}(\Omega)\sum_{\nu_{1}}\kern-15.00002pt\int\frac{\langle\psi_{f}|\hat{O}_{\omega}|\psi_{\nu_{1}}\rangle\langle\psi_{\nu_{1}}|\hat{O}_{\Omega}|\psi_{2p}\rangle}{\omega_{2p}+\omega_{\rm
R}+\Omega-\omega_{\nu_{1}}},$ $\displaystyle a_{\rm XUV+IR+IR}^{(3)}=$
$\displaystyle-\frac{i}{2\pi}\int_{-\infty}^{\infty}d\omega_{\rm
R}\tilde{C}_{2s}(\omega_{\rm
R};\tau)\int_{0}^{\infty}d\omega\int_{0}^{\infty}d\Omega
e^{i(\omega_{f}-\Omega-\omega_{2s}-\omega_{\rm R})\tau}\times$ (20c)
$\displaystyle\tilde{E}_{\rm
IR}(\omega_{f}-\omega-\Omega-\omega_{2s}-\omega_{\rm R})\tilde{E}_{\rm
IR}(\omega)\tilde{E}_{\Omega_{2q+1}}(\Omega)\sum_{\nu_{1},\nu_{2}}\kern-15.00002pt\int\frac{\langle\psi_{f}|\hat{O}_{\omega}|\psi_{\nu_{2}}\rangle\langle\psi_{\nu_{2}}|\hat{O}_{\omega}|\psi_{\nu_{1}}\rangle\langle\psi_{\nu_{1}}|\hat{O}_{\Omega}|\psi_{2s}\rangle}{(\omega_{2s}+\omega_{\rm
R}+\Omega+\omega-\omega_{\nu_{2}})(\omega_{2s}+\omega_{\rm
R}+\Omega-\omega_{\nu_{1}})}$
$\displaystyle-\frac{i}{2\pi}\int_{-\infty}^{\infty}d\omega_{\rm
R}\tilde{C}_{2p}(\omega_{\rm
R};\tau)\int_{0}^{\infty}d\omega\int_{0}^{\infty}d\Omega
e^{i(\omega_{f}-\Omega-\omega_{2p}-\omega_{\rm R})\tau}\times$
$\displaystyle\tilde{E}_{\rm
IR}(\omega_{f}-\omega-\Omega-\omega_{2p}-\omega_{\rm R})\tilde{E}_{\rm
IR}(\omega)\tilde{E}_{\Omega_{2q+1}}(\Omega)\sum_{\nu_{1},\nu_{2}}\kern-15.00002pt\int\frac{\langle\psi_{f}|\hat{O}_{\omega}|\psi_{\nu_{2}}\rangle\langle\psi_{\nu_{2}}|\hat{O}_{\omega}|\psi_{\nu_{1}}\rangle\langle\psi_{\nu_{1}}|\hat{O}_{\Omega}|\psi_{2p}\rangle}{(\omega_{2p}+\omega_{\rm
R}+\Omega+\omega-\omega_{\nu_{2}})(\omega_{2p}+\omega_{\rm
R}+\Omega-\omega_{\nu_{1}})},$

where $\nu_{i}$ identifies the intermediate unperturbed states
$|\psi_{\nu_{i}}\rangle$ of the energy $\omega_{\nu_{i}}$ as
$H_{0}|\psi_{\nu_{i}}\rangle=\omega_{\nu_{i}}|\psi_{\nu_{i}}\rangle$; $\omega$
and $\Omega$ are the frequencies of the IR photon and the $(2q+1)$th-order XUV
harmonics, respectively. Here, we use the relation of Delta-function
$\delta(\omega)=\int_{-\infty}^{\infty}e^{i\omega t}dt/(2\pi)$ and the Fourier
transform (FT):
$\tilde{f}(\omega)=\frac{1}{\sqrt{2\pi}}\int^{\infty}_{-\infty}f(t)e^{i\omega
t}dt$,
$f(t)=\frac{1}{\sqrt{2\pi}}\int^{\infty}_{-\infty}\tilde{f}(\omega)e^{-i\omega
t}d\omega$, and
$\tilde{f}(\omega;\tau)\coloneqq\frac{1}{\sqrt{2\pi}}\int^{\infty}_{-\infty}f(t;\tau)e^{i\omega
t}dt=\tilde{f}(\omega)e^{i\omega\tau}$. The FT of the Rabi amplitudes [Eq.
(14)] are $\tilde{C}_{2s,2p}(\omega_{\rm
R};\tau)\coloneqq\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}C_{2s,2p}(t;\tau)e^{-i\omega_{\rm
R}t}d\omega_{\rm R}=\tilde{C}_{2s,2p}(\omega_{\rm R})e^{i\omega_{\rm R}\tau}$
(see Appendix A), where $\omega_{\rm R}$ is the frequency component of the
Rabi amplitudes. For the long IR field of flat-top envelope [Eq. (3b)], the
Rabi amplitudes mainly distribute at $\omega_{\rm
R}\approx\pm\frac{\Omega_{\rm R}^{0}}{2}$ in the frequency domain. Eq. (20a)
describes the one-photon transition from the zeroth-order Rabi wavefunction
$|\Psi_{\mathcal{R}}\rangle$ to the continuum state $|\phi_{f}\rangle$ by
absorbing an XUV photon of the frequency
$\Omega=\omega_{f}-\omega_{2s,2p}-\omega_{\rm R}$; Eq. (20b) describes the
two-photon transition through the absorption of an XUV photon followed by the
absorption of an IR photon, where the energy-preserving condition
$\omega_{f}=\omega_{2s,2p}+\omega_{\rm R}+\Omega+\omega$ is satisfied; and Eq.
(20c) describes the three-photon transition through the absorption of an XUV
photon followed by the absorption of two IR photons, where the energy-
preserving condition $\omega_{f}=\omega_{2s,2p}+\omega_{\rm R}+\Omega+2\omega$
is satisfied.

As indicated by the term $\omega_{\rm R}$ in the denominators in Eqs. (20c),
the ionization from the Rabi state $|\Psi_{\mathcal{R}}\rangle$ can be
understood to start from the dressed states of $\rm|\psi_{2s}\rangle$ and
$\rm|\psi_{2p}\rangle$, which spilt into two quasienergies considering the
sign of $\omega_{\rm R}$. More specifically, the first (second) term in Eqs.
(20c) corresponds to the ionization from the dressed 2s (2p) states with the
quasienergies of $\omega_{2s}+\omega_{\rm R}$ ($\omega_{2p}+\omega_{\rm R}$),
where $\omega_{\rm R}\approx-\frac{\Omega_{\rm R}^{0}}{2}$ for
$\rm|\phi_{2s}^{l}\rangle$ ($\rm|\phi_{2p}^{l}\rangle$) and $\omega_{\rm
R}\approx\frac{\Omega_{\rm R}^{0}}{2}$ for $\rm|\phi_{2s}^{h}\rangle$
($\rm|\phi_{2p}^{h}\rangle$), respectively. In the photoelectron spectra, the
peak of a lower (higher) photoelectron energy in the AT doublet, namely $\rm
P^{l}$ ($\rm P^{h}$), is formed by the interference among all the ionization
pathways initiated from $\rm|\phi_{2s}^{l}\rangle$ and
$\rm|\phi_{2p}^{l}\rangle$ ($\rm|\phi_{2s}^{h}\rangle$ and
$\rm|\phi_{2p}^{h}\rangle$).

###  III.2 Ionization pathways

As an example, we display the ionization pathways contributing to $\rm P^{h}$
with the photoelectron energy of $2q\omega-I^{2s}_{p}+\omega_{\rm R}$
($\omega_{\rm R}\approx\frac{\Omega_{\rm R}^{0}}{2}$) in Fig. 2. Pathways 1 to
4 are initiated from $\rm|\phi_{2s}^{h}\rangle$, among which pathways 1 and 3
are reminiscent of the well-known two-photon transitions in the traditional
RABBIT cases Klünder  _et al._ (2011). Pathways 1 and 2 denote the absorption
of one $(2q-1)$th-order XUV harmonic $\rm\Omega_{2q-1}$ and one IR photon
$\omega$ in different time orders. Note that, as indicated by the term $H_{\rm
int}^{\perp\mathcal{R}}$ in Eq. (18), $\psi_{2p}$ needs to be excluded from
the intermediate summation in calculating the transition amplitude of pathway
2 since this transition has been handled in the non-perturbative treatment of
the Rabi oscillations. Pathways 3 and 4 denote the absorption of one
$(2q+1)$th-order XUV harmonic $\rm\Omega_{2q+1}$ and the emission of one IR
photon $\omega$ in different time orders. In our case, because the 2p state is
also populated during Rabi oscillations, there are also pathways initiated
from $\rm|\phi_{2p}^{h}\rangle$, as denoted by pathways 5 to 11 in Fig. 2.
Pathway 5 indicates the absorption of one $(2q-1)$th-order XUV harmonic
$\rm\Omega_{2q-1}$. Pathways 6 to 8 refer to the absorption of one $(2q-3)$th-
order XUV harmonic $\rm\Omega_{2q-3}$ and two IR photons in different time
orders. Pathways 9 to 11 denote the absorption of one XUV photon
$\rm\Omega_{2q+1}$ and the emission of two IR photons in different time
orders. Similarly as pathway 2, in calculating the transition amplitude of
pathways 10 and 11, $\psi_{2s}$ needs to be excluded from the first
intermediate summation. Pathways of essential contribution are pathways 1, 3,
5, 6, and 9, whose numbers are highlighted in red in Fig. (2).

In our calculations, the XUV and flat-top IR fields are both spectrally narrow
enough and display no spectral overlap. Therefore, the convolution can be
reasonably dismissed and both laser fields are approximated as monochromatic:
let $E_{\rm IR}(t;\tau)=E_{\omega}\cos[\omega(t-\tau)]$ and $E_{\Omega_{\rm
2q+1}}(t)=E_{\rm 2q+1}\cos(\Omega_{\rm 2q+1}t-\phi_{\rm 2q+1})$, with the
amplitude $E_{\rm 2q+1}$ and the phase $\phi_{\rm 2q+1}$ associated to the
$(2q+1)$th-order XUV harmonic; and let $\tilde{C}_{2s}(\omega_{\rm
R};\tau)=\sqrt{\frac{\pi}{2}}e^{i\omega_{\rm R}\tau}[\delta(\omega_{\rm
R}+\frac{|\Omega_{\rm R}^{0}|}{2})+\delta(\omega_{\rm R}-\frac{|\Omega_{\rm
R}^{0}|}{2})]$ and $\tilde{C}_{2p}(\omega_{\rm
R};\tau)=-\sqrt{\frac{\pi}{2}}\frac{\Omega_{\rm R}^{0}}{|\Omega_{\rm
R}^{0}|}e^{i(\omega+\omega_{\rm R})\tau}[\delta(\omega_{\rm
R}+\frac{|\Omega_{\rm R}^{0}|}{2})-\delta(\omega_{\rm R}-\frac{|\Omega_{\rm
R}^{0}|}{2})]$. Then each term in Eqs. (20c) can be approximately calculated
by multiplying the corresponding IR ($E_{\omega}$) and XUV ($E_{\rm 2q+1}$)
electric field amplitudes to the residual integration of wavefunctions. For
example, the ionization amplitudes $\mathcal{A}_{\mathcal{P}}$ for SB2q, where
$\mathcal{P}$ characterizes the specific pathway, i.e.,
$\mathcal{P}\in\\{\textrm{pathway}\,\,(i)|i=1,\,\cdots,11\\}$, are given as
(see Appendix C for the other pathways)

$\mathcal{A}_{(1)}(\vec{k}_{\pm},\tau)=-\frac{i\pi}{4}e^{i\left[\left(\omega\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}\right)\tau+\phi_{2q-1}\right]}E_{\omega}E_{\rm
2q-1}\sum_{\nu_{1}}\kern-15.00002pt\int\frac{\langle\psi_{f\pm}|\hat{O}_{\omega}|\psi_{\nu_{1}}\rangle\langle\psi_{\nu_{1}}|\hat{O}_{\Omega}|\psi_{2s}\rangle}{\omega_{2s}\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}+\Omega_{2q-1}-\omega_{\nu_{1}}},$
(21a)
$\mathcal{A}_{(3)}(\vec{k}_{\pm},\tau)=-\frac{i\pi}{4}e^{-i\left[\left(\omega\mp\frac{|\Omega_{\textrm{R}}^{0}|}{2}\right)\tau-\phi_{2q+1}\right]}E_{\omega}E_{\rm
2q+1}\sum_{\nu_{1}}\kern-15.00002pt\int\frac{\langle\psi_{f\pm}|\hat{O}^{\dagger}_{\omega}|\psi_{\nu_{1}}\rangle\langle\psi_{\nu_{1}}|\hat{O}_{\Omega}|\psi_{2s}\rangle}{\omega_{2s}\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}+\Omega_{2q+1}-\omega_{\nu_{1}}},$
(21b) $\mathcal{A}_{(5)}(\vec{k}_{\pm},\tau)=\mp\frac{i\pi\Omega_{\rm
R}^{0}}{2|\Omega_{\rm
R}^{0}|}e^{i\left[\left(\omega\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}\right)\tau+\phi_{2q-1}\right]}E_{\rm
2q-1}\langle\psi_{f\pm}|\hat{O}_{\Omega}|\psi_{2p}\rangle,$ (21c)
$\displaystyle\mathcal{A}_{(6)}(\vec{k}_{\pm},\tau)=$
$\displaystyle\mp\frac{i\pi\Omega_{\rm R}^{0}}{8|\Omega_{\rm
R}^{0}|}e^{i\left[\left(3\omega\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}\right)\tau+\phi_{2q-3}\right]}E_{\omega}^{2}E_{\rm
2q-3}\times$ (21d)
$\displaystyle\sum_{\nu_{1},\nu_{2}}\kern-15.00002pt\int\frac{\langle\psi_{f\pm}|\hat{O}_{\omega}|\psi_{\nu_{2}}\rangle\langle\psi_{\nu_{2}}|\hat{O}_{\omega}|\psi_{\nu_{1}}\rangle\langle\psi_{\nu_{1}}|\hat{O}_{\Omega}|\psi_{2p}\rangle}{(\omega_{2p}\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}+\Omega_{2q-3}+\omega-\omega_{\nu_{2}})(\omega_{2p}\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}+\Omega_{2q-3}-\omega_{\nu_{1}})},$
$\displaystyle\mathcal{A}_{(9)}(\vec{k}_{\pm},\tau)=$
$\displaystyle\mp\frac{i\pi\Omega_{\rm R}^{0}}{8|\Omega_{\rm
R}^{0}|}e^{-i\left[\left(\omega\mp\frac{|\Omega_{\textrm{R}}^{0}|}{2}\right)\tau-\phi_{2q+1}\right]}E_{\omega}^{2}E_{\rm
2q+1}\times$ (21e)
$\displaystyle\sum_{\nu_{1},\nu_{2}}\kern-15.00002pt\int\frac{\langle\psi_{f\pm}|\hat{O}^{\dagger}_{\omega}|\psi_{\nu_{2}}\rangle\langle\psi_{\nu_{2}}|\hat{O}^{\dagger}_{\omega}|\psi_{\nu_{1}}\rangle\langle\psi_{\nu_{1}}|\hat{O}_{\Omega}|\psi_{2p}\rangle}{(\omega_{2p}\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}+\Omega_{2q+1}-\omega-\omega_{\nu_{2}})(\omega_{2p}\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}+\Omega_{2q+1}-\omega_{\nu_{1}})},$

where $\psi_{f+}$ ($\psi_{f-}$) is the final continuum state related to the
peak $\rm P^{h}$ ($\rm P^{l}$) of the AT doublet. The energy and the
asymptotic momentum of $\psi_{f\pm}$ are $E_{\pm}=2q\omega-
I_{p}^{2s}\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}=k_{\pm}^{2}/2$ and
$\vec{k}_{\pm}\coloneqq k_{\pm}\hat{k}_{\pm}$, respectively. Here, as
indicated by Eqs. (12) and (14), the zeroth-order Rabi wavefunction
$\Psi_{\mathcal{R}}$ evolving with the time-delayed IR field leads to the
factors $e^{\pm i\frac{|\Omega_{\textrm{R}}^{0}|}{2}\tau}$ and
$e^{i(\omega\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2})\tau}$ in Eqs. (21e),
respectively for the ionization pathways from the 2s and 2p states. In
addition, the initial phase of the electron wave packet generated from
$\rm|\phi_{2p}^{l}\rangle$ has an additional term of $\pi$ compared to that
for $\rm|\phi_{2p}^{h}\rangle$, as implied by the beginning sign $\mp$ of the
ionization amplitudes for pathways 5 to 11. These different initial phases
between the dressed 2p states stem from the sine-like Rabi amplitude of the 2p
state, as indicated by Eq. (14).

Case |  XUV field |  IR field |  P1 |  P2 |  P3 |  P4 |  P5 |  P6 |  P7 |  P8 |  P9 |  P10 |  P11 |  $\varphi-$oscillations |  Phase matched
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
1 | $\Delta m=+1$ | $\Delta m=+1$ | 2 | 2 | 0 | 0 | 2 | 4 | 4 | 4 | 0 | 0 | 0 |  $2\varphi$,$4\varphi$ | Yes
2 | $\Delta m=+1$ | $\Delta m=-1$ | 0 | 0 | 2 | 2 | 0 | -2 | -2 | -2 | 2 | 2 | 2 |  $2\varphi$,$4\varphi$ | Yes
3 | $\Delta m=-1$ | $\Delta m=+1$ | 0 | 0 | -2 | -2 | 0 | 2 | 2 | 2 | -2 | -2 | -2 |  $2\varphi$,$4\varphi$ | Yes
4 | $\Delta m=-1$ | $\Delta m=-1$ | -2 | -2 | 0 | 0 | -2 | -4 | -4 | -4 | 0 | 0 | 0 |  $2\varphi$,$4\varphi$ | Yes
5 | $\Delta m=0$ | $\Delta m=+1$ | 1 | 1 | -1 | -1 | 1 | 3 | 3 | 3 | -1 | -1 | -1 |  $2\varphi$,$4\varphi$ | No
6 | $\Delta m=0$ | $\Delta m=-1$ | -1 | -1 | 1 | 1 | -1 | -3 | -3 | -3 | 1 | 1 | 1 |  $2\varphi$,$4\varphi$ | No
7 | $\Delta m=+1$ | $\Delta m=0$ | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 | No | No
8 | $\Delta m=-1$ | $\Delta m=0$ | -1 | -1 | -1 | -1 | -1 | -1 | -1 | -1 | -1 | -1 | -1 | No | No
9 | $\Delta m=0$ | $\Delta m=0$ | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | No | Yes
Table 2: The magnetic quantum number of the final continuum state ($M$) for
pathways 1 to 11 (denoted as P1 to P11), the oscillation components in SBs as
a function of azimuthal angle $\varphi$, and the phase matching situation in
experiments, under nine combinations of the laser field polarizations (denoted
as Cases 1 to 9). The polarization of the fields is characterized by the
dipole selection rules for the magnetic quantum number ($\Delta m$): $\Delta
m=0$, $\Delta m=+1$, and $\Delta m=-1$ denote the linearly, left-hand
circularly, and right-hand circularly polarized fields, respectively.

In calculation, the incoming final continuum state
$\psi_{\vec{k}_{\pm}}^{-}(\vec{r})\coloneqq\langle\vec{r}|\psi_{f\pm}\rangle$
Breit and Bethe (1954); Altshuler (1956) can be further expanded on the
partial wave series as Starace (1982)

$\displaystyle\psi_{\vec{k}_{\pm}}^{-}(\vec{r})$
$\displaystyle=\frac{1}{k_{\pm}^{1/2}}\sum_{L=0}^{\infty}\sum_{M=-L}^{L}i^{L}e^{-i(\sigma_{L\pm}+\delta_{L\pm})}$
(22) $\displaystyle\times
Y^{\ast}_{L,M}(\hat{k}_{\pm})R_{E_{\pm},L}(r)Y_{L,M}(\hat{r}),$

where $L$ and $M$ are the angular momentum quantum number and the magnetic
quantum number of the partial wave, respectively. Here the phase shift due to
the short-range potential is $\delta_{L\pm}$ in lithium atom Sarsa  _et al._
(2004). The Coulombic phase is
$\sigma_{L\pm}\coloneqq\arg[\Gamma(1+L-iZ/k_{\pm})]$ with the effective
nuclear charge $Z=1$. The energy-normalized radial wavefunction is
$R_{E_{\pm}L}(r)$ with its asymptotic behavior of $\sqrt{2/(\pi
k_{\pm}r)}\sin[k_{\pm}r-L\pi/2-Z\ln(2k_{\pm}r)/k_{\pm}+\sigma_{L\pm}+\delta_{L\pm}]$
when $r\rightarrow\infty$.

###  III.3 The Rabi-RABBIT scheme using circularly polarized laser fields

In the traditional RABBIT scheme with linearly polarized XUV and IR fields
Klünder  _et al._ (2011), absorbing (emitting) an IR photon from the retarded
($\tau>0$) IR field $\textbf{E}_{\rm IR}(t-\tau)$ with respect to the XUV
field $\textbf{E}_{\rm XUV}(t)$ contributes a phase like $e^{+i\omega\tau}$
($e^{-i\omega\tau}$) with $\omega>0$ (see footnote 111Alternatively, absorbing
an XUV photon from the XUV field $\textbf{E}_{\rm XUV}(t+\tau)$ in advance of
the IR field $\textbf{E}_{\rm IR}(t)$ corresponds to an interaction phase of
$e^{-i\Omega_{\rm 2q+1}\tau}$ with $\Omega_{\rm 2q+1}>0$. The equivalence of
these two perspectives is established upon the energy-preserving condition
during the ionization process.) . Therefore, the number of exchanged IR
photons in each ionization pathway is imprinted in the $\tau$-related phase of
the outgoing electron wave packet. Thus, the relative phase of the outgoing
electron wave packets at the same energy via different pathways, named as
RABBIT phase, can be retrieved from the $2\omega\tau$-modulation of SBs Toma
and Muller (2002). As discussed in Ref. Sörngård  _et al._ (2020), by using
different combinations of the XUV and IR fields with either linear or circular
polarizations (see footnote 222Note that in the coordinate system used here,
the $xOy$ plane is defined as the polarization plane of the circularly
polarized fields; and the $z$ axis is defined as the propagation direction of
the circularly polarized fields and the polarization axis of the linearly
polarized fields. This definition of the coordinate system for the circularly
polarized fields differs from that used in our TDSE calculations [Eqs. (3a)
and (4a)], however, the observed physical quantities are unchanged under the
rotation
[$\mathcal{C}:(\hat{e}_{x},\hat{e}_{y},\hat{e}_{z})\rightarrow\mathcal{C}^{\prime}:(\hat{e}_{z},\hat{e}_{x},\hat{e}_{y})$]
of the coordinate system fixed with an observer.) , it is possible to extract
the RABBIT phases from the modulation of SBs as a function of azimuthal angle
$\varphi$. The mechanism behind this alternative way (see footnote
333According to the dipole selection rules, the involved ionization channels
are different in the cases of linearly and circularly polarized fields and
thus the retrieved RABBIT phases are inequivalent.) to extract the RABBIT
phases is based on the dipole selection rules: exchanging one circularly
polarized IR photon will increase or decrease the magnetic quantum number
($M$). Finally, the number of exchanged circularly polarized IR photons is
encoded in the $\varphi$-related phase of the outgoing wave packet via
$e^{iM\varphi}$ originating from $Y_{L,M}^{\ast}(\hat{k}_{\pm})$ in Eq. (22).

In principle, we have nine possible combinations of the polarizations of the
XUV and IR fields, as listed by Cases 1 to 9 in Tab. 2. The corresponding
transition operators $\hat{O}_{\omega}$ and $\hat{O}_{\Omega}$ used in Eqs.
(21e) are $r\cos\theta$, $r\sin\theta e^{i\varphi}/\sqrt{2}$, and $r\sin\theta
e^{-i\varphi}/\sqrt{2}$ for linear, left-hand circular, and right-hand
circular polarizations, respectively Arfken and Weber (1999). Note that phase
matching is only achievable in experiments for Cases 1, 2, 3, 4, and 9. Tab. 2
gives that, if only the circularly polarized IR field is used (Cases 1 to 6),
the RABBIT phases can be successfully extracted from the
$2\varphi-$oscillations of SBs without scan of the time delay between the XUV
and IR fields. In addition, Cases 1 to 6 all give the same interference scheme
as using both linearly polarized XUV and IR fields (Case 9). More
specifically, $2\varphi$-signal ($4\varphi$-signal) results from the
interference of pathways 1, 2, 5 versus 3, 4, 6, 7, 8, 9, 10, 11 (pathways 3,
4, 9, 10, 11 versus 6, 7, 8). As a demonstration, we show as follows how to
extract the RABBIT phases from the $2\varphi-$oscillations in the Rabi-RABBIT
scheme for Case 1 (Sec. IV.1) and Case 2 (Sec. IV.2).

Substituting Eq. (22) into Eqs. (21e), then the ionization amplitude of
pathway $i$ can be written as

$\mathcal{A}_{(i)}(\vec{k}_{\pm},\tau)=\sum_{\mathcal{Q}}\mathcal{A}_{(i),\mathcal{Q}}(\vec{k}_{\pm},\tau).$
(23)

Here $\mathcal{A}_{(i),\mathcal{Q}}$ is the amplitude of a specified
ionization channel for pathway $i\in\\{1,2,\,\cdots,11\\}$, which is
unambiguously characterized by the quantum numbers throughout all the states
in transition. In $N-$photon transition, the ensemble of these quantum numbers
is defined as
$\mathcal{Q}=(l_{i},m_{i}),(\lambda_{1},\mu_{1}),\ldots,(\lambda_{N-1},\mu_{N-1}),(L,M)$,
where $l_{i}$, $\lambda_{N-1}$, and $L$ ($m_{i}$, $\mu_{N-1}$, and $M$)
respectively label the angular (magnetic) quantum numbers of the initial, the
$(N-1)$th intermediate, and the final states. The specific quantum numbers in
$\mathcal{Q}$ are determined by the selection rules, taking Case 1 as an
example, we have $\Delta m=+1$ ($\Delta m=-1$) for absorption (emission) with
$l\geq m$. For a single time-delay ($\tau=0$) RABBIT measurement in Case 1,
$\mathcal{A}_{(i),\mathcal{Q}}$ can be written as

$\mathcal{A}_{(i),\mathcal{Q}}(\vec{k}_{\pm},\tau=0)=\frac{\pi}{k_{\pm}^{1/2}}Y_{L,M}(\hat{k}_{\pm})\mathcal{M}_{(i),\mathcal{Q}}(E_{\pm}),$
(24)

where $\hat{k}_{\pm}=(\theta,\varphi)$ indicates the emission direction of
photoelectrons. The amplitudes $\mathcal{M}_{(i),\mathcal{Q}}(E_{\pm})$ can be
further separated into its angular and radial parts in coordinate
representation as (see Appendix C for the other pathways)

$\displaystyle\mathcal{M}_{(1),\mathcal{Q}}(E_{\pm})$
$\displaystyle=\frac{\pi}{3}E_{\omega}E_{\rm
2q-1}i^{-(L+1)}e^{i(\sigma_{L\pm}+\delta_{L\pm}+\phi_{2q-1})}\langle
Y_{L,M}|Y_{1,1}|Y_{\lambda_{1},\mu_{1}}\rangle\langle
Y_{\lambda_{1},\mu_{1}}|Y_{1,1}|Y_{0,0}\rangle$ (25a)
$\displaystyle\times\sum_{\nu_{1}}\kern-15.00002pt\int\frac{\langle
R_{E_{\pm},L}|r|R_{\nu_{1},\lambda_{1}}\rangle\langle
R_{\nu_{1},\lambda_{1}}|r|R_{2,0}\rangle}{\omega_{2s}\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}+\Omega_{2q-1}-\omega_{\nu_{1}}},$
$\displaystyle\mathcal{M}_{(3),\mathcal{Q}}(E_{\pm})$
$\displaystyle=-\frac{\pi}{3}E_{\omega}E_{\rm
2q+1}i^{-(L+1)}e^{i(\sigma_{L\pm}+\delta_{L\pm}+\phi_{2q+1})}\langle
Y_{L,M}|Y_{1,-1}|Y_{\lambda_{1},\mu_{1}}\rangle\langle
Y_{\lambda_{1},\mu_{1}}|Y_{1,1}|Y_{0,0}\rangle$ (25b)
$\displaystyle\times\sum_{\nu_{1}}\kern-15.00002pt\int\frac{\langle
R_{E_{\pm},L}|r|R_{\nu_{1},\lambda_{1}}\rangle\langle
R_{\nu_{1},\lambda_{1}}|r|R_{2,0}\rangle}{\omega_{2s}\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}+\Omega_{2q+1}-\omega_{\nu_{1}}},$
$\mathcal{M}_{(5),\mathcal{Q}}(E_{\pm})=-\left(\frac{\pi}{3}\right)^{\frac{1}{2}}\frac{\Omega_{\rm
R}^{0}}{|\Omega_{\rm R}^{0}|}E_{\rm 2q-1}i^{-(L\pm
1)}e^{i(\sigma_{L\pm}+\delta_{L\pm}+\phi_{2q-1})}\langle
Y_{L,M}|Y_{1,1}|Y_{1,1}\rangle\langle R_{E_{\pm},L}|r|R_{2,1}\rangle,$ (25c)
$\displaystyle\mathcal{M}_{(6),\mathcal{Q}}(E_{\pm})$
$\displaystyle=-\left(\frac{\pi}{3}\right)^{\frac{3}{2}}\frac{\Omega_{\rm
R}^{0}}{|\Omega_{\rm R}^{0}|}E_{\omega}^{2}E_{\rm 2q-3}i^{-(L\pm
1)}e^{i(\sigma_{L\pm}+\delta_{L\pm}+\phi_{2q-3})}\langle
Y_{L,M}|Y_{1,1}|Y_{\lambda_{2},\mu_{2}}\rangle\langle
Y_{\lambda_{2},\mu_{2}}|Y_{1,1}|Y_{\lambda_{1},\mu_{1}}\rangle\langle
Y_{\lambda_{1},\mu_{1}}|Y_{1,1}|Y_{1,1}\rangle$ (25d)
$\displaystyle\times\sum_{\nu_{1},\nu_{2}}\kern-15.00002pt\int\frac{\langle
R_{E_{\pm},L}|r|R_{\nu_{2},\lambda_{2}}\rangle\langle
R_{\nu_{2},\lambda_{2}}|r|R_{\nu_{1},\lambda_{1}}\rangle\langle
R_{\nu_{1},\lambda_{1}}|r|R_{2,1}\rangle}{(\omega_{2p}\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}+\Omega_{2q-3}+\omega-\omega_{\nu_{2}})(\omega_{2p}\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}+\Omega_{2q-3}-\omega_{\nu_{1}})},$
$\displaystyle\mathcal{M}_{(9),\mathcal{Q}}(E_{\pm})$
$\displaystyle=-\left(\frac{\pi}{3}\right)^{\frac{3}{2}}\frac{\Omega_{\rm
R}^{0}}{|\Omega_{\rm R}^{0}|}E_{\omega}^{2}E_{\rm 2q+1}i^{-(L\pm
1)}e^{i(\sigma_{L\pm}+\delta_{L\pm}+\phi_{2q+1})}\langle
Y_{L,M}|Y_{1,-1}|Y_{\lambda_{2},\mu_{2}}\rangle\langle
Y_{\lambda_{2},\mu_{2}}|Y_{1,-1}|Y_{\lambda_{1},\mu_{1}}\rangle\langle
Y_{\lambda_{1},\mu_{1}}|Y_{1,1}|Y_{1,1}\rangle$ (25e)
$\displaystyle\times\sum_{\nu_{1},\nu_{2}}\kern-15.00002pt\int\frac{\langle
R_{E_{\pm},L}|r|R_{\nu_{2},\lambda_{2}}\rangle\langle
R_{\nu_{2},\lambda_{2}}|r|R_{\nu_{1},\lambda_{1}}\rangle\langle
R_{\nu_{1},\lambda_{1}}|r|R_{2,1}\rangle}{(\omega_{2p}\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}+\Omega_{2q+1}-\omega-\omega_{\nu_{2}})(\omega_{2p}\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}+\Omega_{2q+1}-\omega_{\nu_{1}})},$

where $\hat{O}_{\Omega}$ and $\hat{O}_{\omega}$ [$\hat{O}^{\dagger}_{\Omega}$
and $\hat{O}^{\dagger}_{\omega}$] are replaced by
$-\sqrt{\frac{4\pi}{3}}rY_{1,1}(\hat{r})$
[$\sqrt{\frac{4\pi}{3}}rY_{1,-1}(\hat{r})$] for absorbing [emitting] each
photon. Here the wavefunctions in Eqs. (21e) are separated into their radial
part and spherical harmonics in coordinate presentation:
$\langle\vec{r}|\psi_{2s}\rangle=Y_{0,0}(\hat{r})R_{2,0}(r)$,
$\langle\vec{r}|\psi_{2p}\rangle=Y_{1,1}(\hat{r})R_{2,1}(r)$, and
$\langle\vec{r}|\psi_{\nu_{i}}\rangle=Y_{\lambda_{i},\mu_{i}}(\hat{r})R_{\nu_{i},\lambda_{i}}(r)$
with $\hat{r}=(\theta_{r},\varphi_{r})$.

In our calculation, for the two- and three-photon transition amplitudes in
Eqs. (25e), the infinite summation in the radial part is evaluated with the
Dalgarno-Lewis method Dalgarno  _et al._ (1955). Then the radial part can be
calculated using perturbed wavefunctions Pi and Starace (2014), which satisfy
the inhomogeneous equation with its boundary conditions described in Ref.
Aymar and Crance (1981). Particularly, the integration of two continuum
wavefunctions appearing in Eqs. (25e) is calculated by using the complex
coordinate rotation method Gao and Starace (1987). Additionally, verification
of the matrix elements has been performed by the extrapolation method Cormier
and Lambropoulos (1995). We select some numerical results of dipole transition
matrix elements in Appendix D, where their radial and angular parts are
separately shown.

Figure 3: (a) The schematic of the essential states of the lithium atom
coupled to the IR field. The double-side (single-side) red arrows denote the
strong (weak) coupling. The thick horizontal solid (dashed) lines indicate the
states that are allowed (forbidden) to be populated according to the dipole
selection rules. The detuning of the 2p-3d transition is
$\Delta\coloneqq\omega-E_{3d}+E_{2p}$, with $E_{2p}$ and $E_{3d}$ the energies
of the 2p and 3d states of the lithium atom, respectively. (b, c, d) The
populations of the 2s (red-solid lines), 2p (orange-dashed lines), 3d states
of lithium (bluish-grey-dotted lines) and their summation (thin-purple-solid
lines) as a function of the atom-field interaction time (the beginning time is
set as zero). The intensities of the left-hand circularly polarized IR field
are (b) $1\times 10^{10}$ (corresponding to the energy spacing of 0.0367 eV),
(c) $3\times 10^{10}$ (corresponding to the energy spacing of 0.0642 eV), and
(d) $1\times 10^{11}\rm W/cm^{2}$ (corresponding to the energy spacing of
0.1154 eV).

##  IV Results

###  IV.1 Case 1: Rabi-RABBIT using the left-hand circularly polarized XUV and
IR fields

Figure 4: From top to bottom: (a1, b1) The RABBIT phases extracted from the
photoelectron spectra within the polarization plane of the left-hand
circularly polarized XUV and IR fields for $\rm P^{l}$ as a function of the
photoelectron energy. (a2, b2) The same as (a1, b1), but for $\rm P^{h}$. (a3,
b3) The relative RABBIT phases between $\rm P^{h}$ and $\rm P^{l}$, as a
function of the photoelectron energy. The different columns correspond to the
results obtained (a1, a2, a3) by solving the TDSE and (b1, b2, b3) by
perturbation theory including all pathways, respectively. The purple circles,
blue rhombuses and orange squares correspond to the IR intensities of $1\times
10^{10}$, $3\times 10^{10}$ and $1\times 10^{11}\rm W/cm^{2}$, respectively.

Figure 3(b) shows the populations of the 2s ($m=0$), 2p ($m=1$), and 3d
($m=2$) states of the lithium atom as a function of the atom-field interaction
time, which are approximately obtained by projecting the time-dependent TDSE
solution to the eigenstates of the field-free Hamiltonian $H_{0}$ near the
zero values of the vector potential of the external field. Within the dipole
approximation, the time-dependent coefficients obtained by projecting the
wavefunction onto the unperturbed eigenstates can be interpreted as
probability amplitudes when using the full Hamiltonian written in the length
gauge Kobe and Smirl (1978). Moreover, when using the velocity gauge, the
calculated projection coefficients coincide with the ones employing the length
gauge near the zeros of the vector potential Vábek  _et al._ (2022). Here the
intensities of the left-hand circularly polarized XUV and IR fields are $\rm
1\times 10^{13}$ and $\rm 1\times 10^{10}W/cm^{2}$, respectively. As shown in
Fig. 3(b), the populations of the 2s and 2p states oscillate out of phase with
the interaction time, as a demonstration of Rabi oscillations Rabi (1937). As
the IR intensity increases, the Rabi oscillations between the 2s and 2p states
become faster, as shown in Figs. 3(c) and 3(d) for the IR intensities of $\rm
3\times 10^{10}$ and $\rm 1\times 10^{11}W/cm^{2}$, respectively. As displayed
in Fig. 3(a), the 3d ($m=2$) state [the 3s state ($m=0$)] is allowed
[forbidden] to be populated through absorbing two IR photons from the 2s state
($m=0$) according to the dipole selection rule $\Delta m=+1$ for the left-hand
circular polarization. Due to the large detuning in the 2p-3d transition,
adopting the circular polarized IR field prevents the population leakage from
the Rabi subspace $\mathcal{R}$, as shown in Figs. 3(b)-3(d).

During the entire interaction time, the Rabi oscillations shown in Figs.
3(b)-3(d) are smooth and uniform for the three IR intensities, respectively
validating the rotating wave approximation in solving the Rabi amplitudes [Eq.
(10)] and the monochromatic approximation for the IR field in Eqs. (21e).
Moreover, there is nearly no dampening in the Rabi oscillations for the three
IR intensities due to the relatively negligible transition (including
ionization) from $\mathcal{R}$ to $\mathcal{S}$, as indicated by the total
population of the 2s, 2p, and 3d states after the laser pulses end in Figs.
3(b)-3(d). Therefore, it is a good approximation to treat the ionization
perturbatively on the top of the Rabi wavefunction $\Psi_{\mathcal{R}}$ [Eq.
(12)] which describes perfectly ongoing Rabi oscillations within $\mathcal{R}$
subspace. In addition, the incomplete Rabi oscillations are owing to the
energy shifts of the 2s and 2p states in the presence of the external fields,
which are neglected in our analytical treatment owing to their small effects.

####  IV.1.1 Comparison of TDSE and perturbation theory (PT): RABBIT phases
extracted from the photoelectron spectra within the polarization plane

Figure 4 shows the RABBIT phases as a function of photoelectron energy, which
are extracted from the $2\varphi$-oscillation of SBs in the photoelectron
spectra within the polarization plane of the XUV and IR fields. Note that the
TDSE simulations and the PT calculations are performed in different
coordinates without changing the physical observables (see footnote
22footnotemark: 2). Here and hereafter the polarization plane refers to the
$xOy$ plane with $\theta=\pi/2$ as defined in the PT calculations. Figures
4(a1) and 4(a2) respectively show the RABBIT phases $\rm\phi^{l}(\pi/2)$ and
$\rm\phi^{h}(\pi/2)$ of the peaks with the lower ($\rm P^{l}$) and higher
($\rm P^{h}$) energies in the AT doublet for the three IR intensities, which
are calculated by solving the TDSE. Despite the tiny energy spacing between
the AT doublet $\rm P^{l}$ and $\rm P^{h}$, $\rm\phi^{l}(\pi/2)$ and
$\rm\phi^{h}(\pi/2)$ behave quite differently as a function of the
photoelectron energy, which reveals the influence of the Rabi dynamics on the
phase of the ionized electron wave packets. For both $\rm\phi^{l}(\pi/2)$ and
$\rm\phi^{h}(\pi/2)$, their specific variations with the photoelectron energy
rely on the IR intensity. In addition, $\rm\phi^{l}(\pi/2)$ and
$\rm\phi^{h}(\pi/2)$ obviously show a change in its slope before and after SB
16. Figure 4(a3) shows the relative RABBIT phase
$\rm\Delta\phi(\pi/2)=\phi^{h}(\pi/2)-\phi^{l}(\pi/2)$ as a function of the
photoelectron energy for the three IR intensities, which are obtained by
solving the TDSE. For all the three IR intensities here, $\Delta\phi(\pi/2)$
varies in a similar way with the photoelectron energy and it is close to
$\pi$. In addition, $\Delta\phi(\pi/2)$ deviates more from $\pi$ for the
higher IR intensity. Likewise, there is an obvious change in the slope of
$\Delta\phi(\pi/2)$ around SB 16. Figures 4(b1), 4(b2), and 4(b3) respectively
show the RABBIT phases $\rm\phi^{l}(\pi/2)$, $\rm\phi^{h}(\pi/2)$, and
$\Delta\phi(\pi/2)$ calculated by PT including all pathways shown in Fig. 2
and using the same laser intensities and polarizations as the TDSE
simulations. By comparing Figs. 4(b1), 4(b2), and 4(b3) separately with Figs.
4(a1), 4(a2), and 4(a3), their agreement validates the PT in Sec. III.1 on
describing the physics process in the Rabi-RABBIT scheme for all the three IR
intensities. In addition, the agreement between the PT and TDSE results are
better at lower IR intensities when the Rabi oscillations within $\mathcal{R}$
are more accurately described by the zeroth-order wavefunction
$\Psi_{\mathcal{R}}$ [Eq. (12)] due to smaller energy shifts of the 2s and 2p
states, less population leakage to $\mathcal{S}$, and less dampening in the
Rabi oscillations [Fig. 3].

To specifically reveal the underlying physics in the Rabi-RABBIT scheme, we
analyze the PT results based on the two reduced subsets containing important
ionization pathways in Fig. 2: Subset five consists of the essential pathways
1, 3, 5, 6, and 9; and Subset three is composed of pathways 1, 3, and 5. As a
demonstration, Figs. 5(a) and 5(b) show the possible ionization channels of
several essential pathways in Case 1 and their relative strengths, which are
determined by the dipole selection rules and the propensity rule in laser-
assisted photoionization Busto  _et al._ (2019); Bertolino  _et al._ (2020).
Pathway 1 has only one partial wave $\varepsilon d_{2}$ (here and hereafter
the continuum partial wave is denoted as $\varepsilon l_{m}$ with $l$ the
angular quantum number and $m$ the magnetic quantum number) via one channel
and pathway 3 has two partial waves $\varepsilon s_{0}$ and $\varepsilon
d_{0}$ via two possible channels, as shown in Fig. 5(a). Pathway 5 (not shown)
has the same partial wave as pathway 1. Pathway 6 has only one partial wave
$\varepsilon g_{4}$ via one channel and pathway 9 has three partial waves
$\varepsilon s_{0}$, $\varepsilon d_{0}$, and $\varepsilon g_{0}$ via four
possible channels, as shown in Fig. 5(b).

The two rightmost columns in Fig. 5 show the RABBIT phases extracted from the
$2\varphi$-oscillations of SBs in the photoelectron spectra within the
polarization plane of the XUV and IR fields ($xOy$ plane with $\theta=\pi/2$)
as a function of the photoelectron energy, which are calculated by PT
including part of pathways in Fig. 2 and using the same laser intensities and
polarizations as the TDSE simulations. Figures 5(c1), 5(c2), and 5(c3)
respectively show the RABBIT phases $\rm\phi^{l}(\pi/2)$,
$\rm\phi^{h}(\pi/2)$, and $\Delta\phi(\pi/2)$ calculated only including
pathways of Subset five, which exhibits an excellent agreement with the PT
results including all pathways plotted in Figs. 4(b1), 4(b2), and 4(b3). Their
further agreement with the TDSE results plotted in Figs. 4(a1), 4(a2), and
4(a3) indicates that the TDSE results can be accurately interpreted by only
analyzing the interference among pathways of Subset five. Our calculations
show that, in the formation of the 2$\varphi$-oscillation signal, the
interference of pathways 3 and 5 has the most dominant contribution throughout
all the photoelectron energy for the three IR intensities, followed by the
interference of the usual RABBIT pathways 3 and 1 which has the second most
important contribution. Furthermore, when comparing the phase of the
2$\varphi$-oscillation signal between the AT doublet $\rm P^{l}$ and $\rm
P^{h}$, the interference of pathways 3 and 5 (3 and 1) will give a $\pi$
(zero) difference between the AT doublet $\rm P^{l}$ and $\rm P^{h}$,
recalling that there is a $\pi$ difference between the initial phases of the
electrons ionized from the dressed 2p states $\rm|\phi_{2p}^{l}\rangle$ and
$\rm|\phi_{2p}^{h}\rangle$ due to the sine-like Rabi amplitude
$\sin(\Omega_{\rm R}^{0}t/2)\sim e^{i\Omega_{\rm R}^{0}t/2}-e^{-i\Omega_{\rm
R}^{0}t/2}$ of the 2p state [Eq. (14)]. Therefore, the dominance of the
interference of pathways 3 and 5 explains why $\Delta\phi(\pi/2)$ is very
close to $\pi$ throughout all photoelectron energies for all the three IR
intensities as shown in Figs. 4(a3), 4(b3), and 5(c3), which is in qualitative
accordance with the explanations only based on Subset three in Ref. Liao  _et
al._ (2022). In contrast to the Rabi-RABBIT scheme using the linearly
polarized fields in Ref. Liao  _et al._ (2022), adopting the circularly
polarized IR field prevents the 3s ($m=0$) from being populated and thus we
have a more perfect two-level Rabi system within $\mathcal{R}$, as mentioned
in the previous section. So here we aim to make a quantitative interpretation
of the energy- and polar-angle-dependent RABBIT phases (time-delays) including
more processes.

Figure 5: (a) The schematic of the ionization channels for pathways 1 and 3
when using the left-hand circularly polarized (denoted as $\sigma_{+}$) XUV
and IR fields. (b) The same as (a), but for pathways 6 and 9. The purple and
red arrows indicate the transition via exchanging the XUV and IR photons,
respectively. The solid (dashed) arrows denote the relatively more (less)
probable transition according to the propensity rule in laser-assisted
photoionization when comparing each step of absorption and emission from the
same state. The partial waves of each ionization channel are illustrated as
the real spherical harmonics represented on polar plots. Pathway 5, having the
same partial waves as pathway 1, is not shown. From top to bottom of the two
rightmost columns: (c1, d1) The RABBIT phases extracted from the photoelectron
spectra within the polarization plane of the circularly polarized fields for
$\rm P^{l}$ as a function of the photoelectron energy. (c2, d2) The same as
(c1, d1), but for $\rm P^{h}$. (c3, d3) The relative RABBIT phases between
$\rm P^{h}$ and $\rm P^{l}$, as a function of the photoelectron energy. The
different columns correspond to the results obtained by perturbation theory
only including (c1, c2, c3) pathways 1, 3, 5, 6, 9 and (d1, d2, d3) pathways
1, 3, 5, respectively. The purple circles, blue rhombuses and orange squares
correspond to the IR intensities of $1\times 10^{10}$, $3\times 10^{10}$ and
$1\times 10^{11}\rm W/cm^{2}$, respectively.

For comparison, Figs. 5(d1), 5(d2), and 5(d3) respectively show the RABBIT
phases $\rm\phi^{l}(\pi/2)$, $\rm\phi^{h}(\pi/2)$, and $\Delta\phi(\pi/2)$
calculated with only pathways of Subset three included. The RABBIT phases
$\rm\phi^{l}(\pi/2)$, $\rm\phi^{h}(\pi/2)$, and $\Delta\phi(\pi/2)$ all change
monotonously with the photoelectron energy for the three IR intensities. In
particular, the relative phase $\Delta\phi(\pi/2)$ deviates more from $\pi$ at
higher photoelectron energies and for higher IR intensities because the
relative contribution of the interference of pathways 3 and 1 with respect to
that between pathways 3 and 5 increases with the photoelectron energy and with
the IR intensity. This relatively increasing contribution of the interference
between pathways 3 and 1 with the photoelectron energy results from a slower
decline in the radial integral of the dipole moment matrix element associated
with pathway 1 starting from the 2s state compared to that of pathway 5
initiated from the 2p state, as implied by the one-photon ionization cross
sections from the 2s and 2p states for the photoelectron energies from 10 eV
to 30 eV shown in Fig. 1. The relatively increasing contribution of the
interference between pathways 3 and 1 with the IR intensity is attributed to a
higher order of exchanging IR photons in the continuum for pathway 1 than
pathway 5. Note that the photoelectrons with higher kinetic energies also
interact more strongly with the laser field.

Compared to the monotonous tendencies of $\rm\phi^{l}(\pi/2)$,
$\rm\phi^{h}(\pi/2)$, and $\Delta\phi(\pi/2)$ in the PT results for Subset
three, the PT results for Subset five show some extra “bending” structures
from SB 14 to SB 20 in the cases of $\rm\phi^{l}(\pi/2)$,
$\rm\phi^{h}(\pi/2)$, and $\Delta\phi(\pi/2)$, which implies the non-
negligible influence of the interference of pathways 5 with pathways 6 and 9.
Indeed, our calculations indicate that these ”bending” structures are
attributed to the competition between the interference of pathways 5 and 6
versus the interference of pathways 1 and 3. Similar to the interference of
pathways 1 and 3 both starting from the 2s state, the interference of pathways
5 and 6 both initiated from the 2p state also gives a zero difference in the
phase of the $2\varphi$-oscillations between the AT doublet $\rm P^{l}$ and
$\rm P^{h}$. However, the interference between the zeroth-IR-order pathway 5
and the third-IR-order pathway 6 has an additional $\sim\pi$ phase compared to
the interference between the two first-IR-order usual RABBIT pathways 1 and 3
due to the interaction phase shift accumulated during the absorption or
emission of IR photons in multiphoton above threshold ionization Bertolino and
Dahlström (2021). In this sense, the interference of pathways 5 and 6 will
cancel in strength with the interference of pathways 1 and 3. As the
photoelectron energy increases, the strength of the interference between
pathways 5 and 6 becomes more comparable with (but still weaker than) the
interference of pathways 1 and 3 from SB 16 to SB 20 for the three IR
intensities, which is mainly determined by the relative amplitudes of the
absorbed XUV harmonics in pathways 3 and 6 [e.g., $E_{2q-3}$, $E_{2q-1}$, and
$E_{2q+1}$ for SB2q in Eqs. (21e)]. Hence, the summed strength of the
interferences of pathways 1, 3 and of pathways 5, 6 increases from SB 10 to SB
14 and then decreases from SB 16 to SB 20. Recall that the most dominant
interference of pathways 3 and 5 contributes a Rabi $\pi$ phase in the
formation of the $2\varphi$-oscillations between the AT doublet $\rm P^{l}$
and $\rm P^{h}$. As a superposition effect of the three pairs of interferences
among pathways 1, 3, 5, and 6, $\Delta\phi(\pi/2)$ is near $\pi$ while it is
driven further from $\pi$ from SB 10 to SB 14 and then becomes closer to $\pi$
from SB 16 to SB 20, as shown in Figs. 4(a3), 4(b3), and 5(c3).

Figure 6: From top to bottom: (a1, b1, c1) The relative polar-angle-resolved
RABBIT phases for $\rm P^{l}$ as a function of polar emission angle of
photoelectrons, with the intensities of the left-hand circularly polarized XUV
and IR fields $\rm 1\times 10^{13}$ and $\rm 1\times 10^{10}W/cm^{2}$,
respectively. (a2, b2, c2) The same as (a1, b1, c1), but for $\rm P^{h}$. (a3,
b3, c3) The relative phases between $\rm P^{h}$ and $\rm P^{l}$, as a function
of the polar emission angle of photoelectrons. The different columns
respectively correspond to the results obtained by perturbation theory
including (a1, a2, a3) all pathways; (b1, b2, b3) only pathways 1, 3, 5, 6,
and 9; and (c1, c2, c3) only pathways 1, 3, and 5. The circles, pluses,
squares, rhombuses, crosses and triangles correspond to SBs 10, 12, 14, 16, 18
and 20, respectively.

####  IV.1.2 Detailed analysis of PT results: polar-angle-resolved RABBIT
phases

As shown in Sec. IV.1.1, our PT in Sec. III.1 is reliable to quantitatively
reproduce the TDSE results, especially for low IR intensities in Case 1.
Compared to the TDSE simulations, PT has the advantage of specifying the
amplitude and the phase of each ionization channel so that the essential
processes can be accurately identified in order to uncover the underlying
mechanism. In the following, we apply PT for more detailed examinations on the
polar-angle-integrated and polar-angle-resolved RABBIT phases in Case 1 (Sec.
IV.1.2) and Case 2 (Sec. IV.2). Our PT calculations for Case 1 show that, the
phases extracted from the $2\varphi$-oscillations in SBs of the polar-angle-
integrated photoelectron spectra (polar-angle-integrated RABBIT phases) have
similar behaviors as the RABBIT phases $\rm\phi^{l}(\pi/2)$,
$\rm\phi^{h}(\pi/2)$, and $\Delta\phi(\pi/2)$ due to the similar competition
among the essential pathways of Subset five. In this sense, only the polar-
angle-resolved RABBIT phases for the lowest IR intensity will be discussed in
this section.

Figure 7: (a) The schematic of the ionization channels for pathways 1 and 3
when using the left-hand circularly polarized (denoted as $\sigma_{+}$) XUV
field and the right-hand circularly polarized (denoted as $\sigma_{-}$) IR
field. (b) The same as (a), but for pathways 6 and 9. The purple and red
arrows indicate the transition via exchanging the XUV and IR photons,
respectively. The solid (dashed) arrows denote the relatively more (less)
probable transition according to Fano’s propensity rule and the propensity
rule in laser-assisted photoionization when comparing each step of absorption
and emission from the same state. The partial waves of each ionization channel
are illustrated as the real spherical harmonics represented on polar plots.
Pathway 5, having the same partial waves as pathway 1, is not shown. From top
to bottom of the two rightmost columns: (c1, d1) The polar-angle-resolved
RABBIT phases for $\rm P^{l}$ as a function of polar emission angle of
photoelectrons, with the intensities of the XUV and the IR fields $\rm 1\times
10^{13}$ and $\rm 1\times 10^{10}W/cm^{2}$, respectively. (c2, d2) The same as
(c1, d1), but for $\rm P^{h}$. (c3, d3) The relative RABBIT phases between
$\rm P^{h}$ and $\rm P^{l}$, as a function of the polar emission angle of
photoelectrons. The different columns correspond to the results obtained by
perturbation theory including (c1, c2, c3) all pathways and (d1, d2, d3) only
pathways 1, 3, 5, respectively. The circles, pluses, squares, rhombuses,
crosses and triangles correspond to SBs 10, 12, 14, 16, 18 and 20,
respectively.

Figure 6 shows the RABBIT phases as a function of polar emission angle
$\theta$ of photoelectrons, which are extracted from the
$2\varphi$-oscillation of SBs in the three-dimensional photoelectron spectra.
Here the intensities of the left-hand circularly polarized XUV and IR fields
are $\rm 1\times 10^{13}$ and $\rm 1\times 10^{10}W/cm^{2}$, respectively.
Figure 6(a1) [Fig. 6(a2)] shows the polar-angle-resolved RABBIT phases
$\rm\phi^{l}(\theta)$ [$\rm\phi^{h}(\theta)$] of the peak $\rm P^{l}$ [$\rm
P^{h}$] in the AT doublet for SBs 10 to 20, which are calculated with all
pathways in Fig. 2 included. The phases $\rm\phi^{l}(\theta)$ and
$\rm\phi^{h}(\theta)$ are given relative to a small polar angle
$\theta_{0}=0.25^{\circ}$. Note that only the partial waves associated with
the zero magnetic quantum number ($M=0$) contribute to the signal along the
propagation direction ($\theta=0^{\circ}$) in the photoelectron spectra and
thus there is no $2\varphi$-oscillation for $\theta=0^{\circ}$. The RABBIT
phases $\rm\phi^{l}(\theta)$ and $\rm\phi^{h}(\theta)$ both gently vary with
polar angle $\theta$ for all SBs. For each SB, the obvious and
counterintuitive discrepancy between the behaviors of $\rm\phi^{l}(\theta)$
and $\rm\phi^{h}(\theta)$ as a function of polar angle $\theta$ uncovers the
influence of the Rabi dynamics on the phase of the ionized electron wave
packets. Besides, the specific behaviors of both $\rm\phi^{l}(\theta)$ and
$\rm\phi^{h}(\theta)$ depend on the photoelectron energy. For SBs 10 to 14,
$\rm\phi^{l}(\theta)$ and $\rm\phi^{h}(\theta)$ both change smoothly and
monotonously with polar angle $\theta$. For SBs 16 to 20, however,
$\rm\phi^{l}(\theta)$ and $\rm\phi^{h}(\theta)$ both exhibit some fluctuations
and bendings in the vicinities of $\theta=20^{\circ}$ and $\theta=80^{\circ}$
as a function of polar angle $\theta$. In addition, the curvatures of their
bendings depend on the photoelectron energy. Figure 6(a3) shows the relative
RABBIT phases $\rm\Delta\phi(\theta)=\rm\phi^{h}(\theta)-\rm\phi^{l}(\theta)$
as a function of polar angle $\theta$ for SBs 10 to 20. The relative phase
$\Delta\phi(\theta)$ is close to $\pi$ at all polar emission angles $\theta$
for all SBs and it increases with polar angle $\theta$ for all SBs. Similarly,
$\Delta\phi(\theta)$ varies smoothly with polar angle $\theta$ for SBs 10 to
14 while $\Delta\phi(\theta)$ shows a downward and an upward bendings
separately near $\theta=20^{\circ}$ and $\theta=80^{\circ}$ for SBs 14 to 20.
In addition, the slope of the increase of $\Delta\phi(\theta)$ becomes sharper
and the curvatures of the twice bendings become bigger as the photoelectron
energy increases. Figures 6(b1), 6(b2), and 6(b3) respectively show the RABBIT
phases $\rm\phi^{l}(\theta)$, $\rm\phi^{h}(\theta)$, and $\Delta\phi(\theta)$
as a function of polar angle $\theta$ for SBs 10 to 20, which are calculated
with only pathways of Subset five included. By comparing the first two columns
in Fig. 6, it is seen that Subset five can already reproduce the results in
the first column (including all pathways in Fig. 2) well, for both the
location of the “bending” structures and their dependence on the photoelectron
energy.

For comparison, we calculate the RABBIT phases $\rm\phi^{l}(\theta)$,
$\rm\phi^{h}(\theta)$, and $\Delta\phi(\theta)$ for SBs 10 to 20 by only
including pathways in Subset three, which are respectively shown in Figs.
6(c1), 6(c2), and 6(c3). In Figs. 6(c1) and 6(c2), $\rm\phi^{l}(\theta)$ and
$\rm\phi^{h}(\theta)$ both decrease monotonously with polar angle $\theta$. In
Fig. 6(c3), $\Delta\phi(\theta)$ is flat as a function of polar angle $\theta$
for all SBs. Our calculations show that the interference between pathways 3
and 5 (3 and 1) has the most (second most) important contribution in the
formation of the $2\varphi$-oscillation signal throughout all emission polar
angles $\theta$ for all SBs. This dominance of the interference between
pathways 3 and 5 explains why $\Delta\phi(\theta)$ is close to a Rabi $\pi$
phase throughout all emission polar angles $\theta$ for all SBs in Figs.
6(a3), 6(b3), and 6(c3), rather than the usual near zero RABBIT phase at high
photoelectron energies in consideration of the interaction phase shift
Bertolino and Dahlström (2021). Furthermore, the relative strength of the
interference between pathways 3 and 1 compared to that between pathways 3 and
5 increases with the photoelectron energy and thus $\Delta\phi(\theta)$
deviates more from $\pi$ at higher photoelectron energies in Fig. 6(c3) [and
at higher IR intensities (not shown)]. For each SB, because pathways 1 and 5
have the same partial wave $\varepsilon d_{2}$, the relative strength of the
interference between pathways 3 and 1 with respect to that between pathways 3
and 5 keeps constant as a function of polar angle $\theta$ and thus
$\Delta\phi(\theta)$ is flat with polar angle $\theta$ in Fig. 6(c3).

The discrepancy between the two rightmost columns suggests that the
interferences between pathways 5, 6 and between pathways 5, 9 are responsible
for the observed “bending” structures in Figs. 6(a3) and 6(b3). Similarly to
the interference of pathways 5 and 6, the interference of pathways 5 and 9
also gives a zero difference in the phase of the $2\varphi$-oscillations
between the AT doublet $\rm P^{l}$ and $\rm P^{h}$ and it also cancels in
strength with the interference of pathways 1 and 3 due to the $\pi$
interaction phase existing in the third-IR-order pathway 9. At small polar
angles $\theta$, the amplitude of the partial wave $\varepsilon g_{4}$ of
pathway 6 is nearly zero [Fig. 5(b)] and so as the strength of the
interference of pathways 5 and 6. Therefore, the behavior of
$\Delta\phi(\theta)$ near $\theta=20^{\circ}$ in Figs. 6(a3) and 6(b3) is
mainly determined by the cancellation on the interference of pathways 1 and 3
by the interference of pathways 5 and 9, which depends on the relative
strength of these two pairs of interferences. For SBs 10 and 12, the relative
strength of the interference between pathways 5 and 9 monotonously increases
with polar angle $\theta$ and thus there is no bending for SBs 10 and 12. For
SBs 14 to 20, their relative strength shows a dip in the vicinity of
$\theta=20^{\circ}$ and this dip becomes deeper with the increasing
photoelectron energy. Correspondingly, the curvatures of the downward bending
around $\theta=20^{\circ}$ are bigger at higher energies for SBs 14 to 20. At
big polar angles $\theta$, the interference of pathways 5 and 6 dominates over
that of pathways 5 and 9 for SBs 16 to 20 (which is mainly determined by the
relative amplitudes of the absorbed XUV harmonics). Hence, the upward shape of
the bending structure near $\theta=80^{\circ}$ for SBs 16 to 20 in Figs. 6(a3)
and 6(b3) is mainly attributed to the cancellation on the interference of
pathways 1 and 3 by the interference of pathways 5 and 6. For SBs 16 to 20,
the relative strength of the interference between pathways 5 and 6
monotonously increases with polar angle $\theta$ and its increase becomes more
gentle as $\theta$ approaches $80^{\circ}$. Moreover, the slope of its
increase in their relative strength is sharper for the SBs at higher energies
and thus the curvature of the upward bending near $\theta=80^{\circ}$ is
bigger with increasing photoelectron energy, as shown in Figs. 6(a3) and
6(b3). Note that for the emission pathways 3 and 9, the partial wave
$\varepsilon d_{0}$ ($\varepsilon g_{0}$) has a node at $\theta\approx
54.58^{\circ}$ (two nodes at $\theta\approx 34.54^{\circ}$ and $\theta\approx
69.85^{\circ}$). Nevertheless, no abrupt phase jumps occur in the interference
of pathway 3 with pathway 1 nor 5 because the $\varepsilon s_{0}$ wave
dominates over the $\varepsilon d_{0}$ wave in pathway 3 according to the
propensity rule in laser-assisted photoionization Busto  _et al._ (2019).
Likewise, in pathway 9 where two IR photons are emitted in the continuum, each
step of emission favors decreasing the electron angular momentum Bertolino
_et al._ (2020) and thus there is also no phase jump in the interference
between pathways 5 and 9.

###  IV.2 Case 2: Rabi-RABBIT using the left-hand circularly polarized XUV
field and right-hand circularly polarized IR field

As shown in Tab. 2, the RABBIT phases can be retrieved from the
$2\varphi-$oscillations of SBs when the circularly polarized IR field is used.
According to the dipole selection rules, the ionization channels differ when
changing the polarizations of the laser fields in the Rabi-RABBIT scheme and
so do the behaviors of the RABBIT phases. As demonstrated in Sec. IV.1, when
the circularly polarized XUV and IR fields are co-rotating, the behaviors of
the RABBIT phases are mainly determined by the interferences among pathways of
Subset five. As a comparison, in this section, we apply PT to examine the
polar-angle-resolved (Sec. IV.2.1) and polar-angle-integrated (Sec. IV.2.2)
RABBIT phases in Case 2 of Tab. 2 where the XUV and IR fields have opposite
circular polarizations (counter-rotating).

####  IV.2.1 Polar-angle-resolved RABBIT phases

Figure 8: From top to bottom: (a1, b1) The polar-angle-integrated RABBIT
phases for $\rm P^{l}$ as a function of the photoelectron energy. (a2, b2) The
same as (a1, b1), but for $\rm P^{h}$. (a3, b3) The relative RABBIT phases
between $\rm P^{h}$ and $\rm P^{l}$, as a function of the photoelectron
energy. The different columns correspond to the results obtained by
perturbation theory including (a1, a2, a3) all pathways and (b1, b2, b3) only
pathways 1, 3, and 5, respectively. The purple circles, blue rhombuses and
orange squares correspond to the IR intensities of $1\times 10^{10}$, $3\times
10^{10}$ and $1\times 10^{11}\rm W/cm^{2}$, respectively.

As a demonstration, Figs. 7 (a) and 7(b) show the possible ionization channels
of several pathways in Case 2 and their relative strengths, which are
determined by the dipole selection rules, Fano’s propensity rule Fano (1985),
and the propensity rule in laser-assisted photoionization Busto  _et al._
(2019); Bertolino  _et al._ (2020). Pathway 1 has two partial waves
$\varepsilon s_{0}$ and $\varepsilon d_{0}$ via two possible channels and
pathway 3 only has one partial wave $\varepsilon d_{2}$ via one channel, as
shown in Fig. 7 (a). Pathway 5 (not shown) has the same two partial waves as
pathway 1. Pathway 6 has two partial waves $\varepsilon d_{-2}$ and
$\varepsilon g_{-2}$ via four possible channels and pathway 9 has two partial
waves $\varepsilon d_{2}$ and $\varepsilon g_{2}$ via four possible channels,
as shown in Fig. 7 (b).

The two rightmost columns in Fig. 7 show the RABBIT phases as a function of
polar emission angle $\theta$ of photoelectrons, which are extracted from the
$2\varphi$-oscillations of SBs in the three-dimensional photoelectron spectra.
Here the intensities of the left-hand circularly polarized XUV field and the
right-hand circularly polarized IR field are $1\times 10^{13}$ and $1\times
10^{10}\rm W/cm^{2}$, respectively. Figure 7(c1) [Fig. 7(c2)] shows the polar-
angle-resolved RABBIT phases $\rm\phi^{l}(\theta)$ [$\rm\phi^{h}(\theta)$] of
the peak $\rm P^{l}$ [$\rm P^{h}$] in the AT doublet for SBs 10 to 20, which
are calculated by including all pathways in Fig. 2. Here $\rm\phi^{l}(\theta)$
and $\rm\phi^{h}(\theta)$ are given relative to a small polar angle
$\theta_{0}=0.25^{\circ}$. In contrast to co-rotating Case 1, when adopting
the counter-rotating XUV and IR fields in Case 2, $\rm\phi^{l}(\theta)$ and
$\rm\phi^{h}(\theta)$ both exhibit an obvious phase jump around $\theta\approx
40^{\circ}$ for all SBs, as shown in Figs. 7(c1) and 7(c2). In addition, the
steepness of the phase jump in $\rm\phi^{h}(\theta)$ is slightly sharper than
that of $\rm\phi^{l}(\theta)$ for all SBs. For both $\rm\phi^{l}(\theta)$ and
$\rm\phi^{h}(\theta)$, the steepness of the phase jump becomes more gentle and
the amplitude of the phase jump becomes smaller with the increasing
photoelectron energy. Figure 7(c3) shows the relative RABBIT phases
$\rm\Delta\phi(\theta)=\rm\phi^{h}(\theta)-\rm\phi^{l}(\theta)$ as a function
of polar angle $\theta$ for SBs 10 to 20. For all SBs, the relative phase
$\Delta\phi(\theta)$ is close to $\pi$ and exhibits a phase jump of around
0.25$\pi$ near $\theta\approx 40^{\circ}$. As the photoelectron energy
increases, the amplitude of the phase jump in $\Delta\phi(\theta)$ becomes
bigger.

Figures 7(d1), 7(d2), and 7(d3) respectively show the RABBIT phases
$\rm\phi^{l}(\theta)$, $\rm\phi^{h}(\theta)$, and $\Delta\phi(\theta)$
calculated by only involving pathways in Subset three. The results agree
excellently with the PT results including all pathways. This agreement implies
that, in Case 2, only the interferences of pathway 3 with pathways 1 and 5
have the determinant influence on the behaviors of the RABBIT phases. Indeed,
different from Case 1 where pathways 5 and 6 only have a single partial wave
(Fig. 5), pathways 5, 6, and 9 all have two partial waves in Case 2. More
importantly, for each of pathways 5, 6, and 9 in Case 2, their two partial
waves will interfere destructively due to the opposite signs of their dipole
moment matrix elements originating from the angular integrals. Thus, the
interferences of pathway 5 with pathways 6 and 9, as coherent summations of
all the ionization channels, are negligible compared to that of the
interferences of pathway 3 with pathways 1 and 5. Focusing on the
interferences among pathways in Subset three, the interference of pathways 3
and 5 dominates over that of pathways 3 and 1 throughout all polar angles
$\theta$ for all SBs except in the vicinity of $\theta\approx 40^{\circ}$
where the interference of pathways 3 and 5 has a minimum in strength.
Correspondingly, the relative phase $\Delta\phi(\theta)$ is close to a Rabi
$\pi$ phase in Figs. 7(c3) and 7(d3). In addition, the abrupt jumps of
$\rm\phi^{l}(\theta)$ and $\rm\phi^{h}(\theta)$ respectively in Figs. 7(d1)
and 7(d2) for all SBs are attributed to the dominance of the $\varepsilon
d_{0}$ wave over the $\varepsilon s_{0}$ wave in pathway 5 according to Fano’s
propensity rule favoring to increase angular quantum number in absorption of
one photon from the dressed 2p states Fano (1985). As the $\varepsilon d_{0}$
wave changes its sign as a function of polar angle around its node
$\theta\approx 54.5^{\circ}$, the $\varepsilon d_{0}$ wave destructively
(constructively) interferes with the nodeless $\varepsilon s_{0}$ wave in
pathway 5 at the angles smaller (bigger) than $54.5^{\circ}$, in reminiscent
of the opposite signs of the dipole moment matrix elements related to the
$\varepsilon s_{0}$ and $\varepsilon d_{0}$ waves. Therefore, the phase jumps
occur in the vicinity of $\theta=40^{\circ}$ smaller than the $\varepsilon
d_{0}$ node at $54.5^{\circ}$ where the interference between pathways 3 and 5
has a minimum in strength. Note that this differs from the case in Ref. Busto
_et al._ (2019) where the phase jumps locate around $\theta=75^{\circ}$ larger
than the $\varepsilon d_{0}$ node at $54.5^{\circ}$ because their $\varepsilon
s_{0}$ and $\varepsilon d_{0}$ waves have the same signs when using the
linearly polarized XUV and IR fields. In Figs. 7(d1) and 7(d2), no phase jump
is observed for neither $\rm\phi^{l}(\theta)$ nor $\rm\phi^{h}(\theta)$ in the
vicinity of $\theta=15^{\circ}$ where the interference of pathways 3 and 1 has
a minimum in strength due to its negligible contribution compared to the
interference of pathways 3 and 5. Furthermore, the exact behaviors of the
phase jumps are mainly determined by the relative strength and the phase of
the interferences of pathway 3 with pathways 1 and 5. In particular, the
slightly sharper steepness of the phase jump in $\rm\phi^{h}(\theta)$ compared
to that of $\rm\phi^{l}(\theta)$ emphasizes the fact that there is a $\pi$
difference between the initial phases of the electrons ionized from the
dressed 2p states $\rm|\phi_{2p}^{l}\rangle$ and $\rm|\phi_{2p}^{h}\rangle$
[Eq. (14)].

####  IV.2.2 Polar-angle-integrated RABBIT phases

Figure 8 shows the RABBIT phases as a function of the photoelectron energy,
which are extracted from the 2$\varphi$-oscillations in the photoelectron
spectra integrated along polar emission angle $\theta$ of photoelectrons.
Figures 8(a1) and 8(a2) respectively show the polar-angle-integrated RABBIT
phases $\rm\bar{\phi}^{l}$ and $\rm\bar{\phi}^{h}$ of the peaks the AT doublet
$\rm P^{l}$ and $\rm P^{h}$ in the AT doublet, which are calculated by
including all pathways in Fig. 2. Here the intensity of the left-hand
circularly polarized XUV field is $1\times 10^{13}\rm W/cm^{2}$. The
intensities of the right-hand circularly polarized IR field are $1\times
10^{10}$, $3\times 10^{10}$, and $1\times 10^{11}\rm W/cm^{2}$, respectively.
In Fig. 8(a1), $\rm\bar{\phi}^{l}$ decreases with the photoelectron energy for
the three IR intensities and it is less positive for the higher IR intensity.
In Fig. 8(a2), $\rm\bar{\phi}^{h}$ is less negative for the higher IR
intensity and it shows different dependence on the photoelectron energy for
different IR intensities. At the highest (lowest) IR intensity,
$\rm\bar{\phi}^{h}$ increases (decreases) with the photoelectron energy. At
the moderate IR intensity, $\rm\bar{\phi}^{h}$ keeps nearly constant with the
photoelectron energy. Figure 8(a3) shows the relative RABBIT phase
$\rm\Delta\bar{\phi}=\bar{\phi}^{h}-\bar{\phi}^{l}$ as a function of the
photoelectron energy for the three IR intensities. For all the three IR
intensities here, $\rm\Delta\bar{\phi}$ is close to a Rabi $\pi$ phase for all
the photoelectron energies, which is explained by the dominance of the
interference between pathways 3 and 5.

As a comparison, Figs. 8(b1), 8(b2), and 8(b3) respectively show the RABBIT
phases $\rm\bar{\phi}^{l}$, $\rm\bar{\phi}^{h}$, and $\Delta\bar{\phi}$
calculated by only including pathways of Subset three, which shows an
excellent agreement with the PT results including all pathways in Figs. 8(a1),
8(a2), and 8(a3). This agreement stresses the determinant roles of the
interferences of pathway 3 with pathways 1 and 5 in Case 2. Furthermore, the
relative contribution of the interference of pathways 1 and 3 with respect to
that of pathways 3 and 5 increases with the photoelectron energy and with the
IR intensity. Correspondingly, the relative phase $\Delta\bar{\phi}$ deviates
more from a Rabi $\pi$ phase at higher photoelectron energies and for higher
IR intensities in Figs. 8(a3) and 8(b3). Moreover, the opposite dependence of
the $\rm\bar{\phi}^{l}$ and $\rm\bar{\phi}^{h}$ on the photoelectron energy
shown in Figs. 8(a2) and 8(b2) is attributed to the $\pi$ difference between
the initial phases of the electrons ionized from the dressed 2p states
$\rm|\phi_{2p}^{l}\rangle$ and $\rm|\phi_{2p}^{h}\rangle$.

##  V CONCLUSION

We have carefully examined the RABBIT technique applied to a Rabi-cycling atom
by utilizing circularly polarized XUV and IR fields. In the RABBIT
measurements, the circularly polarized IR field induced Rabi oscillations
between the 2s and 2p states of lithium, which ensured minimal population
leakage to the other bound states. Besides, the use of circularly polarized
laser fields circumvented repetitive RABBIT measurements by scanning the time
delay between the XUV and IR fields, which may make experimental verification
more feasible for the proposed phase-matched Cases 1 and 2. In a single time-
delay RABBIT measurement, the interference phase of the ionized electron wave
packets was retrieved from the photoemission anisotropy along the azimuthal
direction. As expected, both polar-angle-integrated and polar-angle-resolved
photoelectron spectra exhibited a near $\pi$ phase difference between each AT
doublet, as a manifestation of the modulating populations of the two Rabi
states Liao  _et al._ (2022). In addition, adopting the circular polarizations
exhibits a potential to steer the phase of the emitted electron wave packets
by controlling the ionization channels. By using co- or counter-rotating XUV
and IR fields, the polar-angle-integrated and polar-angle-resolved RABBIT
phases exhibited different behaviors, which was traced back to different
competition among partial waves determined by propensity rules Fano (1985);
Busto  _et al._ (2019); Bertolino  _et al._ (2020). In the co-rotating case,
the polar-angle-integrated RABBIT phases showed a bending structure as a
function of the photoelectron energy while the polar-angle-resolved RABBIT
phases showed no phase jumps as a function of polar angle. These complex
phenomena were captured by the competition among the essential ionization
channels of Subset five. In the counter-rotating case, the polar-angle-
integrated RABBIT phases varied monotonously with the photoelectron energy
while the polar-angle-resolved RABBIT phases showed phase jumps near
$\theta=40^{\circ}$, which were explained by the competition among the
essential ionization channels of Subset three.

We discussed how the interference pattern in photoelectron spectra encoded the
information of the modulating populations of the two Rabi states. We believe
the underlying physics can be generalized to the Rabi process studied by other
similar interferometric schemes Mao  _et al._ (2023) and in more complicated
molecular Pan  _et al._ (2023); Sun and Lou (2003) and solid systems Fischer
_et al._ (1998); Silva and Jiménez-Galán (2022). Besides, the pump-probe
interferometry used here can also be adopted to investigate the modulating
population in a multi-level system Wu (1996); Sargent and Horwitz (1976); Song
_et al._ (2016). In addition, the analysis of the interference among multiple
pathways is also applicable in other interferometric schemes Laurent  _et al._
(2012); Bharti  _et al._ (2021); Maroju  _et al._ (2023).

Furthermore, the periodicity of the RABBIT process is determined by both the
physical process of interest and the measurement process. In the RABBIT
measurement on the autoionization process Zhao and Lin (2005), the broken
periodicity leads to a broadening Gruson  _et al._ (2016) and a frequency
modulation Jiménez-Galán  _et al._ (2014, 2016) of the spectral SB peaks. In
the RABBIT measurement of Rabi oscillations, the periodic Rabi process
manifests itself as a splitting of the observed SB peaks in the energy domain,
and the periodicity of the RABBIT process is characterized commonly by the
Rabi frequency and the laser probe frequency. When these two frequencies are
incommensurable, there is no repetitive interference pattern as a function of
the time delay between the two fields in the photoelectron spectra due to the
non-periodic RABBIT process. This breakdown of the periodicity becomes visible
when the two SB peaks overlap in the AT doublet by using the laser fields with
broader spectral widths (e.g., the Rainbow RABBIT method Gruson  _et al._
(2016)). In this case, assuming the IR field is monochromatic, the modulations
are like $|\Omega_{\textrm{R}}^{0}|\tau$ within the overlap of the two SB
peaks Jiang  _et al._ (2021) while the modulations are like $2\omega\tau$
outside of their overlap according to Eqs. (21e). Furthermore, the modulation
frequency of SB peaks should include the finite-pulse effect Jiménez-Galán
_et al._ (2014, 2016) when both laser fields are spectrally broad.

###### Acknowledgements.

Y. L. acknowledges stimulating discussions with Zdeněk Mašín, Serguei
Patchkovskii, Zhen Zhao, Mattias Bertolino, David Busto, Dominik Hoff, and
Anne L’Huillier. This work was supported by National Key Research and
Development Program of China (Grant No. 2019YFA0308300), National Natural
Science Foundation of China (Grants No. 12374264, No. 12021004 and No.
92250303), and Shaanxi Science Foundation (Grant No. 2022JC-DW-05). J. M. D.
acknowledges support from the Knut and Alice Wallenberg Foundation (Grant No.
2019.0154) and the Olle Engkvist’s Foundation (Grant No. 194-0734). The
computing work in this paper is supported by the Public Service Platform of
High Performance Computing provided by Network and Computing Center of HUST.

##  Appendix A Solving the Schrödinger equation within the Rabi subspace

In the two-level Rabi subspace $\mathcal{R}$, the Schrödinger equation is
written as

$i\frac{d}{dt}|\Psi_{\mathcal{R}}(t)\rangle=H_{\mathcal{R}}|\Psi_{\mathcal{R}}(t)\rangle,$
(26)

where the Hamiltonian $H_{\mathcal{R}}$ is defined in Eq. (9). The ansatz of
this equation is a superposition of the 2s and 2p states

$|\Psi_{\mathcal{R}}(t)\rangle=C_{2s}(t;\tau)e^{-i\omega_{2s}t}|\psi_{2s}\rangle+C_{2p}(t;\tau)e^{-i\omega_{2p}t}|\psi_{2p}\rangle,$
(27)

which encodes the unitary transform
$U_{r}=\exp(i\omega_{2s}t)|\psi_{2s}\rangle\langle\psi_{2s}|+\exp(i\omega_{2p}t)|\psi_{2p}\rangle\langle\psi_{2p}|$.
In the condition at resonance $\omega=\omega_{2p}-\omega_{2s}$, the full
transform can be decomposed into two successive unitary transforms as
$U_{r}=S_{0}S_{r}$, where
$S_{0}=\exp[i(\omega_{2s}+\omega_{2p})t/2]\mathbb{I}$ subtracts the global
energy offset $(\omega_{2s}+\omega_{2p})/2$ with
$\mathbb{I}=|\psi_{2s}\rangle\langle\psi_{2s}|+|\psi_{2p}\rangle\langle\psi_{2p}|$
and where $S_{r}=\exp(i\omega\hat{\sigma}_{z}/2)$ encodes the change into a
rotating frame Rabi  _et al._ (1954); Cohen-Tannoudji  _et al._ (1986) with
$\hat{\sigma}_{z}=|\psi_{2s}\rangle\langle\psi_{2s}|-|\psi_{2p}\rangle\langle\psi_{2p}|$
the pauli matrix for $\mathcal{R}$ subspace.

Inserting Eq. (27) into Eq. (26), the coefficient equations are obtained as

$\frac{d}{dt}\bar{Y}(t)=\bar{A}(t)\bar{Y}(t),$ (28a) where the coefficient
matrix is $\bar{Y}(t)\coloneqq\begin{pmatrix}C_{2s}(t;\tau)\\\
C_{2p}(t;\tau)\end{pmatrix},$ (28b) and the interaction matrix is
$\bar{A}(t)\coloneqq-\frac{i}{2}\Omega_{\rm R}(t;\tau)\bar{M}$ (28c) with the
skew hermitian matrix $\bar{M}\coloneqq\begin{pmatrix}0&e^{-i\omega\tau}\\\
e^{i\omega\tau}&0\end{pmatrix}.$ (28d)

Note here that $\Omega_{\rm R}(t;\tau)=\Omega^{\ast}_{\rm R}(t;\tau)$ is used
considering that $\Omega_{\rm R}(t;\tau)$ is real number. According to Magnus
expansion Blanes  _et al._ (2009), the solution of Eq. (28a) is

$\bar{Y}(t)=U^{\prime}_{\mathcal{R}}(t,t_{0})\bar{Y}(t_{0}),$ (29a) with the
propagator $\displaystyle U^{\prime}_{\mathcal{R}}(t,t_{0})=\exp$ (29b)
$\displaystyle\left(\int_{t_{0}}^{t}\bar{A}(t_{1})dt_{1}-\frac{1}{2}\int_{t_{0}}^{t}dt_{1}\int_{t_{0}}^{t_{1}}dt_{2}\left[\bar{A}(t_{2}),\bar{A}(t_{1})\right]+\cdots\right).$

Here $U^{\prime}_{\mathcal{R}}$ is the propagator corresponding to a unitary-
transformed Hamiltonian
$H^{\prime}_{\mathcal{R}}=U_{r}H_{\mathcal{R}}U_{r}^{\dagger}-iU_{r}\partial_{t}U_{r}^{\dagger}$.
Due to the commutation relation
$\left[\bar{A}(t),\bar{A}(t^{\prime})\right]=0$, only the first-order term is
kept in the exponential, i.e.,

$\displaystyle
U^{\prime}_{\mathcal{R}}(t,t_{0})=\exp\left(\int_{t_{0}}^{t}\bar{A}(t^{\prime})dt^{\prime}\right)$
(30)
$\displaystyle=\sum_{n=0}^{\infty}\frac{1}{n!}\left(\int_{t_{0}}^{t}\bar{A}(t^{\prime})dt^{\prime}\right)^{n}$
$\displaystyle=\sum_{n=0}^{\infty}\frac{1}{n!}\left(-\frac{i}{2}\int_{t_{0}}^{t}\Omega_{\rm
R}(t^{\prime};\tau)dt^{\prime}\right)^{n}\bar{M}^{n}.$

Note that $\bar{M}$ satisfies

$\bar{M}^{n}=\begin{cases}\bar{M},&\rm for\,\,odd\,\,n\\\ \mathbb{I},&\rm
for\,\,even\,\,n,\end{cases}$ (31)

and thus

$\displaystyle U^{\prime}_{\mathcal{R}}(t,t_{0})=$
$\displaystyle\cos\left[\frac{1}{2}\int_{t_{0}}^{t}\Omega_{\rm
R}(t^{\prime};\tau)\ dt^{\prime}\right]\mathbb{I}$ (32)
$\displaystyle-i\sin\left[\frac{1}{2}\int_{t_{0}}^{t}\Omega_{\rm
R}(t^{\prime};\tau)\ dt^{\prime}\right]\bar{M}.$

Note that Eq. (32) can also be derived from Eq. (30) by directly using the
relations

$\bar{M}=\vec{\sigma}\cdot\hat{n},$ (33a)
$\displaystyle\exp\left[-i\frac{1}{2}\int_{t_{0}}^{t}\Omega_{\rm
R}(t^{\prime};\tau)\ dt^{\prime}\vec{\sigma}\cdot\hat{n}\right]=$ (33b)
$\displaystyle\cos\left[\frac{1}{2}\int_{t_{0}}^{t}\Omega_{\rm
R}(t^{\prime};\tau)\ dt^{\prime}\right]\mathbb{I}$
$\displaystyle-i\sin\left[\frac{1}{2}\int_{t_{0}}^{t}\Omega_{\rm
R}(t^{\prime};\tau)\ dt^{\prime}\right]\vec{\sigma}\cdot\hat{n},$

where $\vec{\sigma}=[\hat{\sigma}_{x},\hat{\sigma}_{y},\hat{\sigma}_{z}]$ is
related to pauli matrices for $\mathcal{R}$ subspace and where
$\hat{n}=[\cos(\omega\tau),\sin(\omega\tau),0]$ is the unity vector within the
$xOy$ plane in the representation of Bloch sphere Rabi  _et al._ (1954);
Harutyunyan  _et al._ (2023); Torosov and Vitanov (2019).

Then the Rabi amplitudes in accordance with the area theorem Eberly (1998) are
given as

$C_{2s}(t;\tau)=\cos\left[\frac{1}{2}\int_{t_{0}}^{t}\Omega_{\rm
R}(t^{\prime};\tau)dt^{\prime}\right],$ (34a)
$C_{2p}(t;\tau)=-ie^{i\omega\tau}\sin\left[\frac{1}{2}\int_{t_{0}}^{t}\Omega_{\rm
R}(t^{\prime};\tau)dt^{\prime}\right],$ (34b)

which is obtained by using
$\bar{Y}(t)=U^{\prime}_{\mathcal{R}}(t,t_{0})\bar{Y}(t_{0})$ with
$\bar{Y}(t_{0})$ corresponding the boundary condition of $C_{2s}(t_{0})=1$ and
$C_{2p}(t_{0})=0$. Here we take the initial time to infinity
$t_{0}\rightarrow-\infty$. Furthermore, the Fourier transform of the Rabi
coefficients is defined as

$\displaystyle\tilde{C}_{2s}(\omega_{\rm
R};\tau)\coloneqq\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}C_{2s}(t;\tau)e^{-i\omega_{\rm
R}t}d\omega_{\rm R}$ (35a)
$\displaystyle=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}\cos\left[\frac{1}{2}\int_{-\infty}^{t}\Omega_{\rm
R}(t^{\prime}-\tau)dt^{\prime}\right]e^{-i\omega_{\rm R}t}d\omega_{\rm R}$
$\displaystyle=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}\cos\left[\frac{1}{2}\int_{-\infty}^{t-\tau}\Omega_{\rm
R}(t^{\prime})dt^{\prime}\right]e^{-i\omega_{\rm R}(t-\tau)}d\omega_{\rm
R}e^{i\omega_{\rm R}\tau}$ $\displaystyle=\tilde{C}_{2s}(\omega_{\rm
R})e^{i\omega_{\rm R}\tau},$ $\displaystyle\tilde{C}_{2p}(\omega_{\rm
R};\tau)\coloneqq\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}C_{2p}(t;\tau)e^{-i\omega_{\rm
R}t}d\omega_{\rm R}$ (35b)
$\displaystyle=-\frac{ie^{i\omega\tau}}{\sqrt{2\pi}}\int_{-\infty}^{\infty}\sin\left[\frac{1}{2}\int_{-\infty}^{t}\Omega_{\rm
R}(t^{\prime}-\tau)dt^{\prime}\right]e^{-i\omega_{\rm R}t}d\omega_{\rm R}$
$\displaystyle=\frac{-ie^{i(\omega+\omega_{\rm
R})\tau}}{\sqrt{2\pi}}\int_{-\infty}^{\infty}\sin\left[\frac{1}{2}\int_{-\infty}^{t-\tau}\Omega_{\rm
R}(t^{\prime})dt^{\prime}\right]e^{-i\omega_{\rm R}(t-\tau)}d\omega_{\rm R}$
$\displaystyle=\tilde{C}_{2p}(\omega_{\rm R})e^{i\omega_{\rm R}\tau}.$

##  Appendix B Derivations of the Dyson series in the Rabi-RABBIT scheme

In verifying that the wavefunction $|\Psi(t)\rangle$ in Eq. (15) satisfies the
Schrödinger equation

$i\frac{d}{dt}|\Psi(t)\rangle=H(t)|\Psi(t)\rangle,$ (36)

we directly substitute into the Schrödinger equation by the ansatz Milošević
(2013); Patchkovskii  _et al._ (2023)

$|\Psi(t)\rangle=U(t,t_{0})|\Psi(t_{0})\rangle.$ (37a) $\displaystyle
U(t,t_{0})=U_{\mathcal{R}}(t,t_{0})-i\int_{t_{0}}^{t}dt^{\prime}U(t,t^{\prime})H_{\rm
int}^{\perp\mathcal{R}}(t^{\prime})U_{\mathcal{R}}(t^{\prime},t_{0}).$ (37b)

Considering the initial condition
$|\Psi_{\mathcal{R}}(t_{0})\rangle=|\Psi(t_{0})\rangle=e^{-i\omega_{2s}t_{0}}|\psi_{2s}\rangle$,
then we obtain

$\displaystyle i\frac{d}{dt}|\Psi(t)\rangle$
$\displaystyle=i\frac{d}{dt}U(t,t_{0})|\Psi(t_{0})\rangle$ (38)
$\displaystyle=i\frac{d}{dt}U_{\mathcal{R}}(t,t_{0})|\Psi_{\mathcal{R}}(t_{0})\rangle+U(t,t)H_{\rm
int}^{\perp\mathcal{R}}(t)U_{\mathcal{R}}(t,t_{0})|\Psi_{\mathcal{R}}(t_{0})\rangle$
$\displaystyle-i\int_{t_{0}}^{t}dt^{\prime}i\frac{d}{dt}U(t,t^{\prime})H_{\rm
int}^{\perp\mathcal{R}}(t^{\prime})U_{\mathcal{R}}(t^{\prime},t_{0})|\Psi_{\mathcal{R}}(t_{0})\rangle-i\int_{t_{0}}^{t}dt^{\prime}H(t)U(t,t^{\prime})H_{\rm
int}^{\perp\mathcal{R}}(t^{\prime})U_{\mathcal{R}}(t^{\prime},t_{0})|\Psi_{\mathcal{R}}(t_{0})\rangle$
$\displaystyle=[H_{\mathcal{R}}(t)+H_{\rm
int}^{\perp\mathcal{R}}(t)]U_{\mathcal{R}}(t,t_{0})|\Psi_{\mathcal{R}}(t_{0})\rangle-
iH(t)\int_{t_{0}}^{t}dt^{\prime}U(t,t^{\prime})H_{\rm
int}^{\perp\mathcal{R}}(t^{\prime})U_{\mathcal{R}}(t^{\prime},t_{0})|\Psi_{\mathcal{R}}(t_{0})\rangle$
$\displaystyle=H(t)[U_{\mathcal{R}}(t,t_{0})-i\int_{t_{0}}^{t}dt^{\prime}U(t,t^{\prime})H_{\rm
int}^{\perp\mathcal{R}}(t^{\prime})U_{\mathcal{R}}(t^{\prime},t_{0})]|\Psi_{\mathcal{R}}(t_{0})\rangle$
$\displaystyle=H(t)U(t,t_{0})|\Psi(t_{0})\rangle$
$\displaystyle=H(t)|\Psi(t)\rangle.$

Note that in our derivation, we use

$i\frac{d}{dt}|\Psi_{\mathcal{R}}(t)\rangle=H_{\mathcal{R}}(t)|\Psi_{\mathcal{R}}(t)\rangle,$
(39a)
$|\Psi_{\mathcal{R}}(t)\rangle=U_{\mathcal{R}}(t,t_{0})|\Psi_{\mathcal{R}}(t_{0})\rangle,$
(39b)
$i\frac{d}{dt}U_{\mathcal{R}}(t,t_{0})|\Psi_{\mathcal{R}}(t_{0})\rangle=H_{\mathcal{R}}(t)U_{\mathcal{R}}(t,t_{0})|\Psi_{\mathcal{R}}(t_{0})\rangle,$
(39c) $i\frac{d}{dt}U(t,t^{\prime})=H(t)U(t,t^{\prime}),$ (39d)

where the propagator of the full Hamiltonian is
$U(t,t^{\prime})\coloneqq\hat{\mathcal{T}}\exp[-i\int^{t}_{t^{\prime}}H(\tau)d\tau]$
with $\hat{\mathcal{T}}$ the time-ordering operator and $U(t,t)=1$.

##  Appendix C The ionization amplitudes of the other pathways

The ionization amplitudes $\mathcal{A}_{(i)}$ and their reduced amplitudes
$\mathcal{M}_{(i),\mathcal{Q}}$ of pathways 2, 4, 7, 8, 10, and 11 in Fig. 2
are

$\mathcal{A}_{(2)}(\vec{k}_{\pm},\tau)=-\frac{i\pi}{4}e^{i\left[\left(\omega\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}\right)\tau+\phi_{2q-1}\right]}E_{\rm
2q-1}E_{\omega}\sum_{\nu_{1}\neq
2p}\kern-15.00002pt\int\frac{\langle\psi_{f\pm}|\hat{O}_{\Omega}|\psi_{\nu_{1}}\rangle\langle\psi_{\nu_{1}}|\hat{O}_{\omega}|\psi_{2s}\rangle}{\omega_{2s}\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}+\omega-\omega_{\nu_{1}}},$
(40a)
$\mathcal{A}_{(4)}(\vec{k}_{\pm},\tau)=-\frac{i\pi}{4}e^{-i\left[\left(\omega\mp\frac{|\Omega_{\textrm{R}}^{0}|}{2}\right)\tau-\phi_{2q+1}\right]}E_{\rm
2q+1}E_{\omega}\sum_{\nu_{1}}\kern-15.00002pt\int\frac{\langle\psi_{f\pm}|\hat{O}_{\Omega}|\psi_{\nu_{1}}\rangle\langle\psi_{\nu_{1}}|\hat{O}^{\dagger}_{\omega}|\psi_{2s}\rangle}{\omega_{2s}\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}-\omega-\omega_{\nu_{1}}},$
(40b) $\displaystyle\mathcal{A}_{(7)}(\vec{k}_{\pm},\tau)=$
$\displaystyle\mp\frac{i\pi}{8}\frac{\Omega_{\rm R}^{0}}{|\Omega_{\rm
R}^{0}|}e^{i\left[\left(3\omega\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}\right)\tau+\phi_{2q-3}\right]}E_{\omega}^{2}E_{\rm
2q-3}\times$ (40c)
$\displaystyle\sum_{\nu_{1},\nu_{2}}\kern-15.00002pt\int\frac{\langle\psi_{f\pm}|\hat{O}_{\omega}|\psi_{\nu_{2}}\rangle\langle\psi_{\nu_{2}}|\hat{O}_{\Omega}|\psi_{\nu_{1}}\rangle\langle\psi_{\nu_{1}}|\hat{O}_{\omega}|\psi_{2p}\rangle}{(\omega_{2p}\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}+\omega+\Omega_{2q-3}-\omega_{\nu_{2}})(\omega_{2p}\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}+\omega-\omega_{\nu_{1}})},$
$\displaystyle\mathcal{A}_{(8)}(\vec{k}_{\pm},\tau)=$
$\displaystyle\mp\frac{i\pi}{8}\frac{\Omega_{\rm R}^{0}}{|\Omega_{\rm
R}^{0}|}e^{i\left[\left(3\omega\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}\right)\tau+\phi_{2q-3}\right]}E_{\omega}^{2}E_{\rm
2q-3}\times$ (40d)
$\displaystyle\sum_{\nu_{1},\nu_{2}}\kern-15.00002pt\int\frac{\langle\psi_{f\pm}|\hat{O}_{\Omega}|\psi_{\nu_{2}}\rangle\langle\psi_{\nu_{2}}|\hat{O}_{\omega}|\psi_{\nu_{1}}\rangle\langle\psi_{\nu_{1}}|\hat{O}_{\omega}|\psi_{2p}\rangle}{(\omega_{2p}\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}+\omega+\omega-\omega_{\nu_{2}})(\omega_{2p}\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}+\omega-\omega_{\nu_{1}})},$
$\displaystyle\mathcal{A}_{(10)}(\vec{k}_{\pm},\tau)=$
$\displaystyle\mp\frac{i\pi}{8}\frac{\Omega_{\rm R}^{0}}{|\Omega_{\rm
R}^{0}|}e^{-i\left[\left(\omega\mp\frac{|\Omega_{\textrm{R}}^{0}|}{2}\right)\tau-\phi_{2q+1}\right]}E_{\omega}^{2}E_{\rm
2q+1}\times$ (40e) $\displaystyle\sum_{\nu_{1}\neq
2s,\nu_{2}}\kern-22.0pt\int\,\,\,\,\frac{\langle\psi_{f\pm}|\hat{O}^{\dagger}_{\omega}|\psi_{\nu_{2}}\rangle\langle\psi_{\nu_{2}}|\hat{O}_{\Omega}|\psi_{\nu_{1}}\rangle\langle\psi_{\nu_{1}}|\hat{O}^{\dagger}_{\omega}|\psi_{2p}\rangle}{(\omega_{2p}\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}-\omega+\Omega_{2q+1}-\omega_{\nu_{2}})(\omega_{2p}\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}-\omega-\omega_{\nu_{1}})},$
$\displaystyle\mathcal{A}_{(11)}(\vec{k}_{\pm},\tau)=$
$\displaystyle\mp\frac{i\pi}{8}\frac{\Omega_{\rm R}^{0}}{|\Omega_{\rm
R}^{0}|}e^{-i\left[\left(\omega\mp\frac{|\Omega_{\textrm{R}}^{0}|}{2}\right)\tau-\phi_{2q+1}\right]}E_{\omega}^{2}E_{\rm
2q+1}\times$ (40f) $\displaystyle\sum_{\nu_{1}\neq
2s,\nu_{2}}\kern-22.0pt\int\,\,\,\,\frac{\langle\psi_{f\pm}|\hat{O}_{\Omega}|\psi_{\nu_{2}}\rangle\langle\psi_{\nu_{2}}|\hat{O}^{\dagger}_{\omega}|\psi_{\nu_{1}}\rangle\langle\psi_{\nu_{1}}|\hat{O}^{\dagger}_{\omega}|\psi_{2p}\rangle}{(\omega_{2p}\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}-\omega-\omega-\omega_{\nu_{2}})(\omega_{2p}\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}-\omega-\omega_{\nu_{1}})},$
$\displaystyle\mathcal{M}_{(2),\mathcal{Q}}(E_{\pm})$
$\displaystyle=\frac{\pi}{3}E_{\rm
2q-1}E_{\omega}i^{-(L+1)}e^{i(\sigma_{L\pm}+\delta_{L\pm}+\phi_{2q-1})}\langle
Y_{L,M}|Y_{1,1}|Y_{\lambda_{1},\mu_{1}}\rangle\langle
Y_{\lambda_{1},\mu_{1}}|Y_{1,1}|Y_{0,0}\rangle$ (41a)
$\displaystyle\times\sum_{(\nu_{1},\lambda_{1})\neq
2p}\kern-26.00009pt\int\frac{\langle
R_{E_{\pm},L}|r|R_{\nu_{1},\lambda_{1}}\rangle\langle
R_{\nu_{1},\lambda_{1}}|r|R_{2,0}\rangle}{\omega_{2s}\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}+\omega-\omega_{\nu_{1}}},$
$\displaystyle\mathcal{M}_{(4),\mathcal{Q}}(E_{\pm})$
$\displaystyle=-\frac{\pi}{3}E_{\rm
2q+1}E_{\omega}i^{-(L+1)}e^{i(\sigma_{L\pm}+\delta_{L\pm}+\phi_{2q+1})}\langle
Y_{L,M}|Y_{1,1}|Y_{\lambda_{1},\mu_{1}}\rangle\langle
Y_{\lambda_{1},\mu_{1}}|Y_{1,-1}|Y_{0,0}\rangle$ (41b)
$\displaystyle\times\sum_{\nu_{1}}\kern-13.00005pt\int\frac{\langle
R_{E_{\pm},L}|r|R_{\nu_{1},\lambda_{1}}\rangle\langle
R_{\nu_{1},\lambda_{1}}|r|R_{2,0}\rangle}{\omega_{2s}\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}-\omega-\omega_{\nu_{1}}},$
$\displaystyle\mathcal{M}_{(7),\mathcal{Q}}(E_{\pm})$
$\displaystyle=-\left(\frac{\pi}{3}\right)^{\frac{3}{2}}\frac{\Omega_{\rm
R}^{0}}{|\Omega_{\rm R}^{0}|}E_{\omega}^{2}E_{\rm 2q-3}i^{-(L\pm
1)}e^{i(\sigma_{L\pm}+\delta_{L\pm}+\phi_{2q-3})}\langle
Y_{L,M}|Y_{1,1}|Y_{\lambda_{2},\mu_{2}}\rangle\langle
Y_{\lambda_{2},\mu_{2}}|Y_{1,1}|Y_{\lambda_{1},\mu_{1}}\rangle\langle
Y_{\lambda_{1},\mu_{1}}|Y_{1,1}|Y_{1,1}\rangle$ (41c)
$\displaystyle\times\sum_{\nu_{1},\nu_{2}}\kern-15.00002pt\int\frac{\langle
R_{E_{\pm},L}|r|R_{\nu_{2},\lambda_{2}}\rangle\langle
R_{\nu_{2},\lambda_{2}}|r|R_{\nu_{1},\lambda_{1}}\rangle\langle
R_{\nu_{1},\lambda_{1}}|r|R_{2,1}\rangle}{(\omega_{2p}\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}+\omega+\Omega_{2q-3}-\omega_{\nu_{2}})(\omega_{2p}\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}+\omega-\omega_{\nu_{1}})},$
$\displaystyle\mathcal{M}_{(8),\mathcal{Q}}(E_{\pm})$
$\displaystyle=-\left(\frac{\pi}{3}\right)^{\frac{3}{2}}\frac{\Omega_{\rm
R}^{0}}{|\Omega_{\rm R}^{0}|}E_{\omega}^{2}E_{\rm 2q-3}i^{-(L\pm
1)}e^{i(\sigma_{L\pm}+\delta_{L\pm}+\phi_{2q-3})}\langle
Y_{L,M}|Y_{1,1}|Y_{\lambda_{2},\mu_{2}}\rangle\langle
Y_{\lambda_{2},\mu_{2}}|Y_{1,1}|Y_{\lambda_{1},\mu_{1}}\rangle\langle
Y_{\lambda_{1},\mu_{1}}|Y_{1,1}|Y_{1,1}\rangle$ (41d)
$\displaystyle\times\sum_{\nu_{1},\nu_{2}}\kern-15.00002pt\int\frac{\langle
R_{E_{\pm},L}|r|R_{\nu_{2},\lambda_{2}}\rangle\langle
R_{\nu_{2},\lambda_{2}}|r|R_{\nu_{1},\lambda_{1}}\rangle\langle
R_{\nu_{1},\lambda_{1}}|r|R_{2,1}\rangle}{(\omega_{2p}\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}+\omega+\omega-\omega_{\nu_{2}})(\omega_{2p}\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}+\omega-\omega_{\nu_{1}})},$
$\displaystyle\mathcal{M}_{(10),\mathcal{Q}}(E_{\pm})$
$\displaystyle=-\left(\frac{\pi}{3}\right)^{\frac{3}{2}}\frac{\Omega_{\rm
R}^{0}}{|\Omega_{\rm R}^{0}|}E_{\omega}^{2}E_{\rm 2q+1}i^{-(L\pm
1)}e^{i(\sigma_{L\pm}+\delta_{L\pm}+\phi_{2q+1})}\langle
Y_{L,M}|Y_{1,-1}|Y_{\lambda_{2},\mu_{2}}\rangle\langle
Y_{\lambda_{2},\mu_{2}}|Y_{1,1}|Y_{\lambda_{1},\mu_{1}}\rangle\langle
Y_{\lambda_{1},\mu_{1}}|Y_{1,-1}|Y_{1,1}\rangle$ (41e)
$\displaystyle\times\sum_{(\nu_{1},\lambda_{1})\neq
2s,\nu_{2}}\kern-31.0001pt\int\,\,\,\,\,\,\,\,\,\,\frac{\langle
R_{E_{\pm},L}|r|R_{\nu_{2},\lambda_{2}}\rangle\langle
R_{\nu_{2},\lambda_{2}}|r|R_{\nu_{1},\lambda_{1}}\rangle\langle
R_{\nu_{1},\lambda_{1}}|r|R_{2,1}\rangle}{(\omega_{2p}\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}-\omega+\Omega_{2q+1}-\omega_{\nu_{2}})(\omega_{2p}\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}-\omega-\omega_{\nu_{1}})},$
$\displaystyle\mathcal{M}_{(11),\mathcal{Q}}(E_{\pm})$
$\displaystyle=-\left(\frac{\pi}{3}\right)^{\frac{3}{2}}\frac{\Omega_{\rm
R}^{0}}{|\Omega_{\rm R}^{0}|}E_{\omega}^{2}E_{\rm 2q+1}i^{-(L\pm
1)}e^{i(\sigma_{L\pm}+\delta_{L\pm}+\phi_{2q+1})}\langle
Y_{L,M}|Y_{1,1}|Y_{\lambda_{2},\mu_{2}}\rangle\langle
Y_{\lambda_{2},\mu_{2}}|Y_{1,-1}|Y_{\lambda_{1},\mu_{1}}\rangle\langle
Y_{\lambda_{1},\mu_{1}}|Y_{1,-1}|Y_{1,1}\rangle$ (41f)
$\displaystyle\times\sum_{(\nu_{1},\lambda_{1})\neq
2s,\nu_{2}}\kern-31.0001pt\int\,\,\,\,\,\,\,\,\,\,\frac{\langle
R_{E_{\pm},L}|r|R_{\nu_{2},\lambda_{2}}\rangle\langle
R_{\nu_{2},\lambda_{2}}|r|R_{\nu_{1},\lambda_{1}}\rangle\langle
R_{\nu_{1},\lambda_{1}}|r|R_{2,1}\rangle}{(\omega_{2p}\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}-\omega-\omega-\omega_{\nu_{2}})(\omega_{2p}\pm\frac{|\Omega_{\textrm{R}}^{0}|}{2}-\omega-\omega_{\nu_{1}})}.$

##  Appendix D The selected numerical results of the dipole transition matrix
elements

Tab. 3 gives the radial and angular integrals of the dipole transition matrix
elements of pathways 1 to 11 in Fig. 2 for the photoelectron energy of
$2q\omega-I^{2s}_{p}+\frac{\Omega_{\rm R}^{0}}{2}$ with $q=5$ and $\Omega_{\rm
R}^{0}=0.0367$ eV (the higher-energy peak of SB 10) in Case 1.

Pathway |  Ionization channel |  Radial integral |  Angular integral
---|---|---|---
1 | $(l_{i},\lambda_{1},L)=(0,1,2)$ | 41.031722524570355 + 53.198466268618446i | 0.365148371670111
2 | $(l_{i},\lambda_{1},L)=(0,1,2)$ | 3.823833748129400 | 0.365148371670111
3 | $(l_{i},\lambda_{1},L)=(0,1,0)$ | 13.544666666666670 + 44.780000000000015i | 0.333333333333333
$(l_{i},\lambda_{1},L)=(0,1,2)$ | -12.984548602944278 - 11.427350863413425i | -0.149071198499986
4 | $(l_{i},\lambda_{1},L)=(0,1,0)$ | -1.398733333333334 | 0.333333333333333
$(l_{i},\lambda_{1},L)=(0,1,2)$ | 0.852016435026670 | -0.149071198499986
5 | $(l_{i},L)=(1,2)$ | 0.213074268742145 | 0.632455532033676
6 | $(l_{i},\lambda_{1},\lambda_{2},L)=(1,2,3,4)$ | 4656.010341994734 + 3491.179677824839i | 0.276026223736942
7 | $(l_{i},\lambda_{1},\lambda_{2},L)=(1,2,3,4)$ | 4082.151822845632 + 1288.738836005407i | 0.276026223736942
8 | $(l_{i},\lambda_{1},\lambda_{2},L)=(1,2,3,4)$ | -839.4509516287872 | 0.276026223736942
9 | $(l_{i},\lambda_{1},\lambda_{2},L)=(1,2,1,0)$ | -2362.517301523949 - 1025.027667919262i | 0.230940107675850
$(l_{i},\lambda_{1},\lambda_{2},L)=(1,2,1,2)$ | 1004.238761782608 - 123.4397252103228i | -0.103279555898864
$(l_{i},\lambda_{1},\lambda_{2},L)=(1,2,3,2)$ | 277.7039715897753 - 23.54641086493669i | -0.044262666813799
$(l_{i},\lambda_{1},\lambda_{2},L)=(1,2,3,4)$ | 169.0151673747687 + 39.07506621875388i | 0.032991443953693
10 | $(l_{i},\lambda_{1},\lambda_{2},L)=(1,0,1,0)$ | 276.6854940046418 + 869.2200752739551i | 0.192450089729875
$(l_{i},\lambda_{1},\lambda_{2},L)=(1,0,1,2)$ | -256.6496964086782 - 221.8186661817862i | -0.086066296582387
$(l_{i},\lambda_{1},\lambda_{2},L)=(1,2,1,0)$ | 20.244209838865040 + 63.312230519334360i | 0.038490017945975
$(l_{i},\lambda_{1},\lambda_{2},L)=(1,2,1,2)$ | -17.338916109487684 - 16.156881592225183i | -0.017213259316477
$(l_{i},\lambda_{1},\lambda_{2},L)=(1,2,3,2)$ | -219.8836499309096 + 75.62719251805706i | -0.044262666813799
$(l_{i},\lambda_{1},\lambda_{2},L)=(1,2,3,4)$ | -98.482759346168690 - 11.215441372057903i | 0.032991443953693
11 | $(l_{i},\lambda_{1},\lambda_{2},L)=(1,0,1,0)$ | -21.007851794913194 | 0.192450089729875
$(l_{i},\lambda_{1},\lambda_{2},L)=(1,0,1,2)$ | 0.800098112918845 | -0.086066296582387
$(l_{i},\lambda_{1},\lambda_{2},L)=(1,2,1,0)$ | 0.049998533311822 | 0.038490017945975
$(l_{i},\lambda_{1},\lambda_{2},L)=(1,2,1,2)$ | -1.710240592647929 | -0.017213259316477
$(l_{i},\lambda_{1},\lambda_{2},L)=(1,2,3,2)$ | -4.063799702841705 | -0.044262666813799
$(l_{i},\lambda_{1},\lambda_{2},L)=(1,2,3,4)$ | -8.459336144166398 | 0.032991443953693
Table 3: The radial and angular integrals of the dipole transition matrix
elements of pathways 1 to 11 for the higher-energy peak of SB 10 in Case 1.

## References

  * Rabi (1937) I. I. Rabi, Phys. Rev. 51, 652 (1937).
  * Pileio  _et al._ (2009) G. Pileio, M. Carravetta, and M. H. Levitt, Phys. Rev. Lett. 103, 083002 (2009).
  * Dudin  _et al._ (2012) Y. O. Dudin, L. Li, F. Bariani, and A. Kuzmich, Nature Physics 8, 790 (2012).
  * Batalov  _et al._ (2008) A. Batalov, C. Zierl, T. Gaebel, P. Neumann, I.-Y. Chan, G. Balasubramanian, P. R. Hemmer, F. Jelezko, and J. Wrachtrup, Phys. Rev. Lett. 100, 077401 (2008).
  * Stievater  _et al._ (2001) T. H. Stievater, X. Li, D. G. Steel, D. Gammon, D. S. Katzer, D. Park, C. Piermarocchi, and L. J. Sham, Phys. Rev. Lett. 87, 133603 (2001).
  * Merlin (2021) R. Merlin, American Journal of Physics 89, 26 (2021).
  * Knight and Milonni (1980) P. Knight and P. Milonni, Physics Reports 66, 21 (1980).
  * Autler and Townes (1955) S. H. Autler and C. H. Townes, Phys. Rev. 100, 703 (1955).
  * Jiang  _et al._ (2021) W.-C. Jiang, H. Liang, S. Wang, L.-Y. Peng, and J. Burgdörfer, Phys. Rev. Research 3, L032052 (2021).
  * Nandi  _et al._ (2022) S. Nandi, E. Olofsson, M. Bertolino, S. Carlström, F. Zapata, D. Busto, C. Callegari, M. Di Fraia, P. Eng-Johnsson, R. Feifel, G. Gallician, M. Gisselbrecht, S. Maclot, L. Neoričić, J. Peschel, O. Plekan, K. C. Prince, R. J. Squibb, S. Zhong, P. V. Demekhin, M. Meyer, C. Miron, L. Badano, M. B. Danailov, L. Giannessi, M. Manfredda, F. Sottocorona, M. Zangrando, and J. M. Dahlström, Nature 608, 488 (2022).
  * Liao  _et al._ (2022) Y. Liao, Y. Zhou, L.-W. Pi, J. Liang, Q. Ke, Y. Zhao, M. Li, and P. Lu, Phys. Rev. A 105, 063110 (2022).
  * Zhang  _et al._ (2022) X. Zhang, Y. Zhou, Y. Liao, Y. Chen, J. Liang, Q. Ke, M. Li, A. Csehi, and P. Lu, Phys. Rev. A 106, 063114 (2022).
  * Olofsson and Dahlström (2023) E. Olofsson and J. M. Dahlström, Phys. Rev. Res. 5, 043017 (2023).
  * Chu and Lin (2013) W.-C. Chu and C. D. Lin, Phys. Rev. A 87, 013415 (2013).
  * Wu  _et al._ (2013) M. Wu, S. Chen, M. B. Gaarde, and K. J. Schafer, Phys. Rev. A 88, 043416 (2013).
  * Ott  _et al._ (2014) C. Ott, A. Kaldun, L. Argenti, P. Raith, K. Meyer, M. Laux, Y. Zhang, A. Blättermann, S. Hagstotz, T. Ding, R. Heck, J. Madroñero, F. Martín, and T. Pfeifer, Nature 516, 374 (2014).
  * Harkema  _et al._ (2021) N. Harkema, C. Cariker, E. Lindroth, L. Argenti, and A. Sandhu, Phys. Rev. Lett. 127, 023202 (2021).
  * Ramsey and Hellwarth (1956) N. F. Ramsey and R. W. Hellwarth, Physics Today 9, 37 (1956).
  * Xu and Heinzen (1999) G. Xu and D. J. Heinzen, Phys. Rev. A 59, R922 (1999).
  * Floquet (1883) G. Floquet, Annales scientifiques de l’École Normale Supérieure 2e série, 12, 47 (1883).
  * Chu and Telnov (2004) S.-I. Chu and D. A. Telnov, Physics Reports 390, 1 (2004).
  * Jaynes and Cummings (1963) E. Jaynes and F. Cummings, Proceedings of the IEEE 51, 89 (1963).
  * Shore and Knight (1993) B. W. Shore and P. L. Knight, Journal of Modern Optics 40, 1195 (1993).
  * Hentschel  _et al._ (2001) M. Hentschel, R. Kienberger, C. Spielmann, G. A. Reider, N. Milosevic, T. Brabec, P. Corkum, U. Heinzmann, M. Drescher, and F. Krausz, Nature 414, 509 (2001).
  * Paul  _et al._ (2001) P. M. Paul, E. S. Toma, P. Breger, G. Mullot, F. Augé, P. Balcou, H. G. Muller, and P. Agostini, Science 292, 1689 (2001).
  * Pazourek  _et al._ (2015) R. Pazourek, S. Nagele, and J. Burgdörfer, Rev. Mod. Phys. 87, 765 (2015).
  * Isinger  _et al._ (2017) M. Isinger, R. J. Squibb, D. Busto, S. Zhong, A. Harth, D. Kroon, S. Nandi, C. L. Arnold, M. Miranda, J. M. Dahlström, E. Lindroth, R. Feifel, M. Gisselbrecht, and A. L’Huillier, Science 358, 893 (2017).
  * Beaulieu  _et al._ (2017) S. Beaulieu, A. Comby, A. Clergerie, J. Caillat, D. Descamps, N. Dudovich, B. Fabre, R. Géneaux, F. Légaré, S. Petit, B. Pons, G. Porat, T. Ruchon, R. Taïeb, V. Blanchet, and Y. Mairesse, Science 358, 1288 (2017).
  * Lucchini  _et al._ (2015) M. Lucchini, L. Castiglioni, L. Kasmi, P. Kliuiev, A. Ludwig, M. Greif, J. Osterwalder, M. Hengsberger, L. Gallmann, and U. Keller, Phys. Rev. Lett. 115, 137401 (2015).
  * Jordan  _et al._ (2020) I. Jordan, M. Huppert, D. Rattenbacher, M. Peper, D. Jelovina, C. Perry, A. von Conta, A. Schild, and H. J. Wörner, Science 369, 974 (2020).
  * Dahlström  _et al._ (2012) J. M. Dahlström, A. L’Huillier, and A. Maquet, J. Phys. B 45, 183001 (2012).
  * Dahlström  _et al._ (2013) J. Dahlström, D. Guénot, K. Klünder, M. Gisselbrecht, J. Mauritsson, A. L’Huillier, A. Maquet, and R. Taïeb, Chemical Physics 414, 53 (2013), attosecond spectroscopy.
  * Toma and Muller (2002) E. S. Toma and H. G. Muller, J. Phys. B 35, 3435 (2002).
  * Véniard  _et al._ (1996) V. Véniard, R. Taïeb, and A. Maquet, Phys. Rev. A 54, 721 (1996).
  * Géneaux  _et al._ (2016) R. Géneaux, A. Camper, T. Auguste, O. Gobert, J. Caillat, R. Taïeb, and T. Ruchon, Nature Communications 7, 12583 (2016).
  * Gruson  _et al._ (2016) V. Gruson, L. Barreau, Á. Jiménez-Galan, F. Risoud, J. Caillat, A. Maquet, B. Carré, F. Lepetit, J.-F. Hergott, T. Ruchon, L. Argenti, R. Taïeb, F. Martín, and P. Salières, Science 354, 734 (2016).
  * Zhong  _et al._ (2020) S. Zhong, J. Vinbladh, D. Busto, R. J. Squibb, M. Isinger, L. Neoričić, H. Laurell, R. Weissenbilder, C. L. Arnold, R. Feifel, J. M. Dahlström, G. Wendin, M. Gisselbrecht, E. Lindroth, and A. L’Huillier, Nature Communications 11, 5042 (2020).
  * Biswas  _et al._ (2020) S. Biswas, B. Förg, L. Ortmann, J. Schötz, W. Schweinberger, T. Zimmermann, L. Pi, D. Baykusheva, H. A. Masood, I. Liontos, A. M. Kamal, N. G. Kling, A. F. Alharbi, M. Alharbi, A. M. Azzeer, G. Hartmann, H. J. Wörner, A. S. Landsman, and M. F. Kling, Nature Physics 16, 778 (2020).
  * Gong  _et al._ (2022) X. Gong, W. Jiang, J. Tong, J. Qiang, P. Lu, H. Ni, R. Lucchese, K. Ueda, and J. Wu, Phys. Rev. X 12, 011002 (2022).
  * Swoboda  _et al._ (2010) M. Swoboda, T. Fordell, K. Klünder, J. M. Dahlström, M. Miranda, C. Buth, K. J. Schafer, J. Mauritsson, A. L’Huillier, and M. Gisselbrecht, Phys. Rev. Lett. 104, 103003 (2010).
  * Kheifets (2021) A. S. Kheifets, Phys. Rev. A 104, L021103 (2021).
  * Autuori  _et al._ (2022) A. Autuori, D. Platzer, M. Lejman, G. Gallician, L. Maëder, A. Covolo, L. Bosse, M. Dalui, D. Bresteau, J.-F. Hergott, O. Tcherbakoff, H. J. B. Marroux, V. Loriot, F. Lépine, L. Poisson, R. Taïeb, J. Caillat, and P. Salières, Science Advances 8, eabl7594 (2022).