#  Brillouin modes in weakly driven dissipative optical lattices: simple
theoretical model vs pump-probe spectroscopy

David Cubero<EMAIL_ADDRESS>Departamento de Física Aplicada I, Escuela
Politécnica Superior, Universidad de Sevilla, Calle Virgen de África 7, 41011
Sevilla, Spain  Kefeng Jiang, Alexander Staron, Casey Scoggins, Daniel
Wingert, Stone Oliver, and Samir Bali<EMAIL_ADDRESS>Department of
Physics, Miami University, Oxford, Ohio 45056-1866, USA

###### Abstract

Atoms confined in a three-dimensional dissipative optical lattice oscillate
inside potential wells, occasionally hopping to adjacent wells, thereby
diffusing in all directions. Illumination by a weak probe beam modulates the
lattice, yielding propagating atomic density waves, referred to as Brillouin
modes which travel perpendicular to the direction of travel of the probe. We
investigate theoretically and experimentally these modes in the case of a
driving potential perturbation whose spatial period is not a multiple of the
period of the underlying optical potential, allowing for a deeper
understanding of Brillouin mode generation in cold confined atoms. The role of
two distinct mechanisms for directed propagation is elucidated, one arising
from a velocity-matching between the propagating modulation and the average
velocity of the atom oscillating inside a well, and the other arising from a
frequency-matching between the the modulation frequency and the oscillation
frequencies.

Light-induced forces on matter over wavelength and sub-wavelength spatial
scales have wide applicability in quantum sensing and metrology [1], ranging
from the novel design of periodic potential landscapes [2, 3, 4] to the
innovative transport and sorting of particles [5]. In particular, considerable
interest has focused on the noise-induced directed motion of particles in the
absence of a net force [6], with special attention devoted to cold atoms
confined in dissipative optical lattices where environmental noise in the form
of spontaneous emission is significant [6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28].

A dissipative optical lattice consists of counter-propagating laser beams
tuned near atomic resonance that yield AC Stark-shifted ground state potential
wells [29]. Atoms in the lattice undergo the well-known process of Sisyphus
cooling, and settle into these wells, where they oscillate with a vibrational
frequency that is determined by the well-depth [11]. The stochastic optical
pumping processes associated with Sisyphus cooling also cause the atoms to
occasionally transfer between adjacent wells, leading to spatial diffusion of
the cold atom sample [11, 30]. The introduction of a weak probe beam along a
symmetry axis of the lattice, results in a time-periodic driving of the
lattice that breaks the symmetry, causing directed atomic density waves to be
set up. The directed propagation proceeds in the absence of a net force. These
propagating atomic density waves are referred to as Brillouin modes in analogy
to acoustic waves rippling through a fluid [8, 13, 15, 28]. Recently, a noise-
induced resonant enhancement of this directed propagation was observed [28],
and a new theory, based on the Fourier decomposition of the current into its
atomic density wave contributions [31], was developed to explain this
stochastic resonance. This theory was also able to successfully predict the
thresholds for the transition to the regime of infinite density in the cold
atom setup [32, 31].

In this paper, we investigate the driving of the lattice by illuminating the
cold atoms with a probe beam propagating at a slight angle to the lattice
symmetry axis. In this case, the spatial period of the driving perturbation is
now not an integer multiple of the period of the underlying lattice potential.
This allows, at least theoretically, for the possibility of the two spatial
driving frequencies to be in irrational ratio, which permits the exploration
of spatial quasiperiodic driving in these systems, in analogy with the time
quasiperiodic case, wherein the lattice is driven by two incommensurate
frequencies [6, 7, 33, 34, 35]. In the present spatial quasiperiodic case,
though, we will show that the generated directed current is not as sensitive
to the nature of the driving as in the time quasiperiodic case, where true
quasiperiodicity is able to suppress the directed motion. Here the transition
from periodicity to quasiperiodicity is not observed to be sharp. However, the
chosen setup, where the spatial periods of the underlying lattice and the
driving are basically uncoupled, sheds light on how the Brillouin modes are
generated.

The paper is organized as follows. In Sec. I we define the system model
studied. New analytical results, based on a Fourier decomposition of the
current, are discussed in Sec. II. Numerical simulations and experiments are
discussed in Sec. III and IV, respectively. Finally, Sec. V ends with the
conclusions.

##  I System models

We consider atoms confined in a so-called 3D-lin$\perp$lin optical lattice
[11], formed by the superposition of four red-detuned laser beams
$\vec{k}_{1-4}$ and frequency $\omega_{l}$ in a tetrahedral configuration, see
Fig. 1(a). For $F_{g}=1/2\rightarrow F_{e}=3/2$ atoms the lattice is formed by
just two light-shifted ground state $\pm 1/2$-spin potentials. An additional
weak probe laser of frequency $\omega_{p}$, forming an angle $\theta_{p}$ with
the $z$-axis and with its polarization parallel to the $y$-axis, is added to
drive the system out of equilibrium and put the atoms in directed motion [8],
$\theta_{p}=0$ in Fig. 1(a, b) and $\theta_{p}\neq 0$ in Fig. 1(c, d). In the
experiments, this model is already a simplification, since the atoms need to
have a more complex transition than $F_{g}=1/2\rightarrow F_{e}=3/2$, but the
theoretical results are still expected to provide good qualitative insight
[36].

Figure 1:  3D-lin$\perp$lin tetrahedral lattice illuminated by a weak probe.
(a) and (b) depict the case of space periodic driving: The probe propagates
along the $z$-axis, which is the lattice symmetry axis. (c) and (d) show the
case where the space period of the driving is not a multiple of the period of
the underlying optical lattice: The probe propagates along a direction forming
an angle $\theta_{p}$ with the $z$-axis. Here,
$|\Delta{\vec{k}_{1}}|=k_{0}+k_{l}\sin\theta_{p}$ and
$|\Delta{\vec{k}_{2}}|=k_{0}-k_{l}\sin\theta_{p}$, where $k_{l}$ is the laser
wave number, and $k_{0}=k_{l}\sin\theta_{x}$.

Following previous studies, [8, 37, 14, 13, 15, 28, 31], we focus on movement
along one of the directions, taken as the $x$-axis. The optical potential
associated with the above setup is then given by (after taking $y=z=0$)

$\displaystyle
U_{\pm}(x,t)=\frac{U_{0}}{2}\Big{[}-\frac{3}{2}-\frac{1}{2}\cos(2k_{0}x)\pm\cos(k_{0}x)$
$\displaystyle+\varepsilon_{p}\cos(-k_{0}x+k_{l}\sin\theta_{p}x-\delta_{p}t)$
$\displaystyle+\varepsilon_{p}\cos(k_{0}x+k_{l}\sin\theta_{p}x-\delta_{p}t)$
$\displaystyle\pm\varepsilon_{p}\cos(k_{l}\sin\theta_{p}x-\delta_{p}t)\Big{]},$
(1)

where $k_{0}=k_{l}\sin\theta_{x}$, $k_{l}=\omega_{l}/c$ is the laser beam wave
number, $\delta_{p}=\omega_{p}-\omega_{l}$ is a small probe detuning relative
to the lattice ($\delta_{p}/\omega_{l}\ll 1$) 111This definition of the probe
detuning $\delta_{p}$ corrects a sign error in Ref. [28].,
$U_{0}=-16\hbar\Delta^{\prime}_{0}/3$, with $\Delta^{\prime}_{0}$ ($<0$) being
the light-shift per lattice field, and $\varepsilon_{p}=E_{p}/(2E_{0})$.

The optical well-depth $U_{0}$ defines a vibrational frequency associated with
an atom of mass $m_{a}$ oscillating at the bottom of a well,

$\Omega_{x}=k_{x}\sqrt{3U_{0}/2m_{a}}=4\sin\theta_{x}\sqrt{|\Delta_{0}^{\prime}|\omega_{r}},$
(2)

where $\omega_{r}=\hbar k_{l}^{2}/(2m_{a})$ is the recoil frequency.

The probe perturbation in (1) appears in three terms: A modulation propagating
to the right with phase velocity
$v_{1}=\delta_{p}/(k_{0}+k_{l}\sin\theta_{p})$, another to the left with
velocity $v_{2}=-\delta_{p}/(k_{0}-k_{l}\sin\theta_{p})$, and a third moving
with velocity $v_{3}=\delta_{p}/(k_{l}\sin\theta_{p})$. Each of these terms
will produce an excitation of atomic density waves on their own. Thus, for the
sake of simplicity, we consider them separately in the theoretical
considerations that follow.

The first system model to consider, denoted as case (a), is the following
optical potential, which accounts for the third probe term in (1), and is
given by

$\displaystyle
U_{\pm}^{(a)}(x,t)=\frac{U_{0}}{2}\Big{[}-\frac{3}{2}-\frac{1}{2}\cos(2k_{0}x)\pm\cos(k_{0}x)$
$\displaystyle\pm\varepsilon_{p}\cos(k_{p}x-\delta_{p}t+\phi_{p}),\Big{]},$
(3)

where $k_{p}=k_{l}\sin\theta_{p}$, and $\phi_{p}$ is a probe phase which has
been introduced for convenience in the analytical calculations to be presented
below. Here the potential perturbation changes sign with the specific atomic
state. Terms of this kind are also generated with a $x$-polarized probe.

In addition, we consider the case of the following optical potential, denoted
as case (b), which accounts for the first two probe terms in (1), and is given
by:

$\displaystyle
U_{\pm}^{(b)}(x,t)=\frac{U_{0}}{2}\Big{[}-\frac{3}{2}-\frac{1}{2}(2k_{0}x)\pm\cos(k_{0}x)$
$\displaystyle+\varepsilon_{p}\cos(k_{p}x-\delta_{p}t+\phi_{p}),\Big{]},$ (4)

By setting $k_{p}=|\Delta{\vec{k}_{1}}|=k_{0}+k_{l}\sin\theta_{p}$, or
$k_{p}=|\Delta{\vec{k}_{2}}|=-(k_{0}-k_{l}\sin\theta_{p})$, we can study the
effects introduced by the first and second probe terms, respectively. If
$\theta_{p}=0$, $|\Delta{\vec{k}_{1,2}}|=|\Delta{\vec{k}}|=k_{0}$. See Figs. 1
(b) and (d) for an illustration of $\Delta{\vec{k}_{1}}$ and
$\Delta{\vec{k}_{2}}$.

Note the on-axis case $\theta_{p}=0$, recently investigated in Refs. [28, 31],
produces a space periodic perturbation, because the wave number of the driving
field $k_{p}$ is the same as that of the underlying lattice, $k_{0}$. On the
other hand, a probe angle $0<\theta_{p}<\pi/2$ such that $k_{p}/k_{0}$ is an
irrational ratio produces a space quasiperiodic drive.

In the semi-classical approximation [36], the atoms in the ground state
$|\pm\rangle$ are described by the phase space density $P_{\pm}(x,p,t)$ at the
position $x$ with momentum $p$, which satisfies the following coupled Fokker-
Planck equations

$\displaystyle\left[\frac{\partial}{\partial
t}+\frac{p}{m_{a}}\frac{\partial}{\partial
x}-U_{\pm}^{\prime}(x)\frac{\partial}{\partial p}\right]P_{\pm}=$
$\displaystyle-\gamma_{\pm}(x)P_{\pm}+\gamma_{\mp}(x)P_{\mp}+\frac{\partial^{2}}{\partial
p^{2}}\left[D_{0}P_{\pm}\right],$ (5)

where $U_{\pm}^{\prime}=\partial U_{\pm}/\partial x$, and

$\gamma_{\pm}(x)=g_{0}\pm g_{1}\cos(k_{0}x)+g_{2}\cos(2k_{0}x)$ (6)

are the transition rates between the ground state sublevels, defined in terms
of $\Gamma^{\prime}$, the photon scattering rate per lattice beam, as
$g_{0}=2\Gamma^{\prime}/3$, $g_{1}=8\Gamma^{\prime}/9$,
$g_{2}=2\Gamma^{\prime}/9$, and $D_{0}=5\hbar^{2}k_{0}^{2}\Gamma^{\prime}/18$
is a noise strength describing the random momentum jumps that result from the
interaction with the photons. As in [31], we are neglecting the probe
contribution to the transition rates, the radiation forces, and noise terms,
since their effect is observed to be small in the simulations.

Directed motion is characterized by the current, defined as the average atomic
velocity

$\displaystyle\langle
v\rangle=\lim_{t\rightarrow\infty}\frac{\langle[x(t)-x(0)]\rangle}{t}=$
$\displaystyle\lim_{t\rightarrow\infty}\frac{1}{t}\int_{0}^{t}\\!\\!dt^{\prime}\int\\!\\!dx\int\\!\\!dp\,\frac{p}{m_{a}}[P_{+}(x,p,t^{\prime})+P_{-}(x,p,t^{\prime})].$
(7)

##  II Analytical results

Following the method presented in Ref. [31], we Fourier-decompose the current
in (7) into contributions arising from atomic density modes excited by the
probe. Using $l$, $n$, $m$ to denote the mode numbers, so that the mode has a
frequency $\omega=l\delta_{p}$ and wave number $k=nk_{0}+mk_{p}$, we obtain,
in terms of the amplitudes of the excited atomic density waves,

$\displaystyle
P^{\pm}[l,n,m]=\frac{\delta_{p}}{2\pi}\int_{0}^{2\pi/\delta_{p}}\\!\\!\\!dt\,e^{-il\delta_{p}t}\int\\!dx\,e^{i(nk_{0}+mk_{p})x}$
$\displaystyle\int\\!dp\,P_{\pm}(x,p,t),$ (8)

the following expansion for the current for case (b), i.e., the setup defined
by (4) for the $U^{(b)}$ system,

$\displaystyle\langle
v\rangle_{\mathrm{(b)}}=\frac{m_{a}}{m_{a}F_{0}g_{1}-2D_{0}k_{0}}\Bigg{[}$
$\displaystyle-\frac{\text{Im}\left[P_{+}[0,1,0]\right]F_{0}\left(8g_{0}^{2}-4g_{2}^{2}/3+F_{0}k_{0}/(2m_{a})\right)}{k_{0}}$
$\displaystyle+\frac{\text{Im}\left[P_{+}[0,2,0]\right]F_{0}\left(-4g_{0}g_{1}{-8g_{1}g_{2}/3}+F_{0}k_{0}/m_{a}\right)}{k_{0}}$
$\displaystyle+\frac{\text{Im}\left[P_{+}[0,3,0]\right]F_{0}\left(-16g_{0}g_{1}/3-2g_{2}^{2}-3F_{0}k_{0}/(2m_{a})\right)}{k_{0}}$
$\displaystyle+\frac{\text{Im}\left[P_{+}[0,4,0]\right]F_{0}\left(-2g_{1}g_{2}/3+F_{0}k_{0}/(2m_{a})\right)}{k_{0}}$
$\displaystyle-\frac{\text{Im}\left[P_{+}[0,5,0]\right]F_{0}2g_{2}^{2}}{3k_{0}}+\frac{\text{Im}\left[e^{i\phi_{p}}P_{+}[1,-2,1]\right]F_{0}F_{p}(2k_{0}-k_{p})}{2m_{a}k_{p}}$
$\displaystyle+\frac{\text{Im}\left[e^{i\phi_{p}}P_{+}[1,-1,1]\right]F_{0}F_{p}(-k_{0}+k_{p})}{m_{a}k_{p}}$
$\displaystyle+\frac{\text{Im}\left[e^{i\phi_{p}}P_{+}[1,0,1]\right]2F_{p}k_{0}\delta_{p}^{2}}{k_{p}^{2}}$
$\displaystyle+\frac{\text{Im}\left[e^{i\phi_{p}}P_{+}[1,1,1]\right]F_{0}F_{p}(k_{0}+k_{p})}{m_{a}k_{p}}$
$\displaystyle-\frac{\text{Im}\left[e^{i\phi_{p}}P_{+}[1,2,1]\right]F_{0}F_{p}(2k_{0}+k_{p})}{2m_{a}k_{p}}$
$\displaystyle+\frac{\text{Im}\left[e^{i2\phi_{p}}P_{+}[2,0,2]\right]F_{p}^{2}k_{0}}{m_{a}k_{p}},\Bigg{]}.$

where the force amplitudes $F_{0}=k_{0}U_{0}/2$ and
$F_{p}=U_{0}\varepsilon_{p}k_{p}/2$. Equation (II) is valid for an arbitrary
$k_{p}$, including the periodic case $k_{p}=k_{0}$, which was validated
numerically in Ref. [31].

Similar expressions are found for case (a), defined by Eq. (3), which are
reported in the appendix. The main difference between the cases (a) and (b),
defined by (3) and (4), respectively, is that case (a) provides different
expressions for the generic case $k_{p}\neq k_{0}$ and the special case
$k_{p}=k_{0}$. In the expansion for the periodic case $k_{p}=k_{0}$, reported
in (A), there are more terms —thus more atomic wave modes activated— than in
the quasiperiodic case where $k_{p}/k_{0}$ is an irrational ratio, reported in
(A). These extra terms cannot be obtained from (A) after taking the limit
$k_{p}\rightarrow k_{0}$.

Furthermore, (A) shows some apparent singularities in the form of coefficients
with denominators proportional to $k_{0}-k_{p}$ or $k_{0}-2k_{p}$, and thus
apparently problematic in the periodic limits $k_{p}\rightarrow k_{0}$ and
$k_{p}\rightarrow 2k_{0}$. This could suggest a special sensitivity to the
periodic/quasiperiodic transition, similar to that observed in the case of
time quasiperiodicity [6, 7, 39, 33, 34, 35], where the large sensitivity in
the system response to time quasiperiodic forces is known to yield sub-Fourier
resonances.

However, our numerical results, reported in the following section, shows that
this is not the case, and the transition is smooth. This is possible because
the amplitudes of the modes involved decay to zero as (or stronger than)
$k_{0}-k_{p}$ or $k_{0}-2k_{p}$, thus removing any possible singularity in the
periodic cases $k_{0}=k_{p}$ or $k_{0}=2k_{p}$.

The analytical expansions, though, still provide a useful decomposition into
atomic waves of directed transport, and are used in the following sections to
interpret the atomic transport provoked by the probe.

##  III Numerical results

Numerical solutions of the equation (5) are obtained by generating a large
number of individual atomic trajectories $x_{j}(\sigma_{j}(t),t)$, where
$\sigma_{j}(t)=+1$ or $-1$ is the occupied state at time $t$ in that
trajectory, using a stochastic algorithm [40]. Averages are computed using
over $1.5\times 10^{6}$ trajectories. Following Ref. [31], the atomic mode
amplitudes (8) are calculated via the formula,

$\displaystyle
P_{\pm}[l,n,m]=\lim_{l^{\prime}\rightarrow\infty}\frac{\delta_{p}}{2\pi
l^{\prime}N}\times$ $\displaystyle\sum_{j=1}^{N}\int_{0}^{2\pi
l^{\prime}/\delta_{p}}\\!\\!\\!dt\,e^{i[(nk_{0}+mk_{p})x_{j}(\sigma_{j}(t),t)-l\delta_{p}t]}\delta_{\sigma_{j}(t),\pm
1}.$ (10)

In all simulations, units are defined such that $m_{a}=\hbar=k_{l}=1$. In
these units, the optical lattice parameters were fixed to $U_{0}=200$,
$\Gamma^{\prime}=2.85$.

We start by first studying the system sensitivity to the transition between
space periodic and quasiperiodic driving. We chose the $U^{(a)}$ system in (3)
(i.e., case (a)), because it gives different expansions in the periodic
$k_{p}=k_{0}$ and quasiperiodic ($k_{p}$ and $k_{0}$ incommensurate) cases.
Specifically, Eq. (A) indicates a potential singularity when $k_{p}\rightarrow
k_{0}$ in the mode $\mathrm{Re}[P_{+}[1,-1,1]]$ due to the presence of
$(k_{p}-k_{0})$ in the denominator of the coefficient.

Figure 2 shows the current and the mode contributions to the current, given by
(A), as a function of $k_{p}$ around the periodic case $k_{p}=k_{0}$. No mode
or contribution is observed to act abruptly, on the contrary, they all are
seen to behave smoothly. The potential singularity in the mode
$\mathrm{Re}[P_{+}[1,-1,1]]$ does not actually occur, because the mode
amplitude itself tends to zero in the limit $k_{p}\rightarrow k_{0}$ as
$(k_{p}-k_{0})$, such that its current contribution is finite, as seen in Fig.
2.

Similar features are observed near $k_{p}=2k_{0}$, demonstrating that the
transition from space quasiperiodicity to periodicity is smooth.

Figure 2:  Smooth space quasiperiodicity. Current and mode contributions to
the current as a function of the driving wave number $k_{p}$ for the $U^{(a)}$
system in (3) (i.e., case (a)), with $\delta_{p}=7.5$, $\epsilon_{p}=0.8$,
$\theta_{x}=30^{0}$ (thus $k_{0}=0.5$), and $\phi_{p}=0$. Each mode $(l,n,m)$
has a frequency $\omega=l\delta_{p}$ and wave number $k=nk_{0}+mk_{p}$. Mode
amplitudes are measured in the simulation via (10), and their precise
contribution to the current determined using (A). The dashed line is the sum
of all current contributions, and the diamonds are the current calculated from
its definition (7). Units are defined such that $m_{a}=\hbar=k_{l}=1$.

Next, in order to understand how Brillouin modes are generated, we study the
atomic waves when varying the driving frequency $\delta_{p}$. A series of
resonances, that is, local maxima at certain values of the driving frequencies
is observed, allowing for a proper rationalization of the transport mechanisms
and the experimental results.

Figure 3 shows the current as a function of $\delta_{p}$ for several values of
the driving amplitude $\varepsilon_{p}$ in the $U^{(b)}$ system defined by (4)
(i.e., case (b)), with space periodic driving $k_{p}=k_{0}$, for several
values of the driving amplitude $\varepsilon_{p}$.

Figure 3:  The atomic current is plotted for the case of space periodic
driving as a function of the driving frequency $\delta_{p}$ in the $U^{(b)}$
system defined by (4) (i.e., case (b)), with $\theta_{x}=25^{0}$,
$\phi_{p}=\pi$, for several values of the driving amplitude $\varepsilon_{p}$.
The driving wave number is fixed to $k_{p}=k_{0}$. Dotted vertical lines are
placed at $\omega_{0}=5.75$, $\omega_{0}/2$, and $2\omega_{0}$. Resonances at
multiples and sub-multiples of the intrinsic frequency $\omega_{0}$ are
observed.

A peak at $\delta_{p}=\omega_{0}=5.75$ is observed, common to all curves,
except the one corresponding to the largest value of $\varepsilon_{p}$, which
is shifted to higher frequency values, about $\delta_{p}\sim 7.5$. This peak
is not far from the vibrational frequency (2) associated with linear
oscillations in the optical wells, $\Omega_{x}=7.3$ with the parameter values
mentioned above for the simulations. In Fig. 3 the presence of additional
peaks at about $\omega_{0}/2$ and $2\omega_{0}$ can be seen. This behavior is
not uncommon, in rocked ratchets the current is also observed [33] to peak at
multiples and sub-multiples, in general a fractional number, of an intrinsic
frequency.

As the driving wave number $k_{p}$ is varied, so does the phase velocity of
the propagating perturbation $v_{p}=\delta_{p}/k_{p}$. It is expected [8] that
a peak is produced when this phase velocity matches the velocity

$v_{0}=\omega_{0}/k_{0}.$ (11)

The intrinsic velocity $v_{0}$ is associated with an average drift in one
direction due to half oscillations in a well, followed by transitions between
the atomic states [8]. This velocity matching ($vm$) mechanism thus yields a
peak at

$(\delta_{p})_{vm}=\omega_{0}k_{p}/k_{0}.$ (12)

The numerical results, presented in Fig. 4 do not contradict this picture.

Figure 4:  Same as in Fig. 3 (i.e., for the $U^{(b)}$ system in (4), case
(b)), but this time for space quasiperiodic driving with
$k_{p0}=k_{0}/\sqrt{5}$. The current is plotted for several values of the
driving wave number, $k_{p}=0.5k_{p0},k_{p0},1.5k_{p0},2.5k_{p0},3.5k_{p0}$,
and $4.5k_{p0}$. In all cases the driving amplitude $\epsilon_{p}$ was varied
so that $\epsilon_{p}k_{p}=0.8$, thus keeping the force amplitude $F_{p}$
fixed. As in Fig. 3, dotted vertical lines are placed at $\omega_{0}=5.75$,
$\omega_{0}/2$, and $2\omega_{0}$, plus an additional one at $3\omega_{0}$.
The truncated vertical dashed lines correspond to the shifted values
$(\delta_{p})_{vm}=\omega_{0}k_{p}/k_{0}$.

However, it is somehow obscured by the above mentioned frequency peaking at a
fractional ratio of the intrinsic frequency. For example, let’s consider the
case of space-quasiperiodic driving with $k_{p}=0.5k_{p0}$, where
$k_{p0}\equiv k_{0}/\sqrt{5}$. In this case, the current-vs-$\delta_{p}$ curve
displays a peak near ($\delta_{p})_{vm}$, but it is also leaning towards
$\omega_{0}/2$. The curves for $k_{p}=k_{p0}$ and $k_{p}=1.5k_{p0}$ show no
velocity matching shift, just peaking at the value $\omega_{0}$. The curve
$k_{p}=2.5k_{p0}$ does peak at ($\delta_{p})_{vm}$, but this value is very
near $\omega_{0}$ in this case. It also shows a small peak at about
$1.5\omega_{0}$. In the case $k_{p}=3.5k_{p0}$, the frequency
($\delta_{p})_{vm}$ lies between the peaks at $\omega_{0}$ and $2\omega_{0}$.
The corresponding curve peaks at $\omega_{0}$, ($\delta_{p})_{vm}$ (or
$1.5\omega_{0}$, since they are very near here), and $2\omega_{0}$. The curve
$k_{p}=4.5k_{p0}$ offers a clear confirmation of the velocity matching
mechanism, because it shows no clear peak at $\omega_{0}$, and a distinct one
at ($\delta_{p})_{vm}$, which is also close to $2\omega_{0}$.

Overall, the discussed shift due to velocity matching is clearly at play in
the system, but the shift takes place through local maxima at a fractional
ratio of the intrinsic frequency $\omega_{0}$.

We may then wonder how a particular atomic mode is affected by the above
discussed resonances. Figure 5 shows the results for the case
$k_{p}=3.5k_{p0}$, which was chosen because the current shows three peaks in
Figure 4, at $\omega_{0}$, $1.5\omega_{0}\approx$ ($\delta_{p})_{vm}$, and
$2\omega_{0}$. The plot shows peaks at these frequency values at most of the
atomic modes, behaving very similarly all of them. Moreover, an extra peak in
most of them is also visible at $\delta_{p}=3\omega_{0}$, a peak which is
difficult to appreciate in the current plot, Fig. 4.

Figure 5:  Current and mode contributions to the current as a function of the
driving frequency $\delta_{p}$ for the quasiperiodic case $k_{p}=3.5k_{p0}$
shown in Fig. 4 (for the $U^{(b)}$ system in (4), case (b)). As in Fig. 4,
dotted vertical lines are placed at $\omega_{0}=5.75$, $\omega_{0}/2$,
$2\omega_{0}$, and $3\omega_{0}$.

##  IV Experiment

Our experiments are performed on 85Rb atoms in a dilutely occupied dissipative
3D lattice in a tetrahedral lin$\perp$lin configuration, as in Fig. 1.
Directed motion is produced by the weak $\hat{y}$-polarized probe which makes
an angle $\theta_{p}$ with the $z$-axis. The probe frequency $\omega_{p}$ is
scanned around the (fixed) frequency $\omega_{l}$ of the lattice beams, and
probe transmission is measured as a function of probe detuning
$\delta=\omega_{p}-\omega_{l}$. As in Ref. [28], $\delta/\omega_{L}<10^{-9}$.
Here, the intensity ratio of probe to lattice (sum of all four beams) is less
than $3\%$,

In accordance with the notation used in Ref. [11, 28],
$\Delta_{0}^{\prime}\equiv\Delta{s_{0}}/2$, and
$\Gamma^{\prime}\equiv\Gamma{s_{0}/2}$ where
$s_{0}\equiv(I/I_{sat})/(1+4\Delta^{2}/\Gamma^{2})$ is just the saturation
parameter. For the $F_{g}=3\rightarrow F_{e}=4$ transition in 85Rb,
$I_{sat}=1.67$ mW/cm2 for $\sigma$-light, $\Gamma/2\pi$ is the natural
linewidth for 85Rb (6.07 MHz), and $\omega_{r}/2\pi$ is the recoil frequency
(3.86 kHz). In our experiment, $\theta_{x}=\theta_{y}=25^{0}$, and each
lattice beam has intensity $I=6.22$ mW/cm2, a $1/e^{2}$-diameter 5.4 mm (the
probe diameter is 1.4 mm), and red-detuning $\Delta=8.75\Gamma$. In order to
determine the intensity $I$ that actually illuminates the atoms, care is taken
to account for the intensity loss through the windows of the vacuum cell that
houses the lattice, and the background Rb vapor. These values yield
$s_{0}=0.012$, $\Gamma^{\prime}/2\pi=36.4$ kHz, and a well-depth
$U_{0}=445\,\hbar\omega_{r}$. Using the definition of the recoil frequency
just after (2), and setting $\hbar=m_{a}=k_{l}=1$ as in the simulations, we
find this $U_{0}$-value corresponds to 223, and $\Gamma^{\prime}$, in units of
the recoil frequency, corresponds to 4.72 - these values are comparable to
those assumed in the simulations in Sec. III.

Figs. 6(a) and (b) show distinctly different probe transmission spectra for
$\theta_{p}=0$ and $17.5^{0}$, respectively. The peaks in the spectra
correspond to photons absorbed from a lattice beam and emitted via stimulated
emission into the probe, while dips correspond to photons absorbed from the
probe and emitted into a lattice beam.

Figure 6:  Probe transmission spectra measured for (a) $\theta_{p}=0$ (space
periodic case $k_{p}=k_{0}$) and (b) $17.5^{0}$ ($k_{p}$ not a multiple of
$k_{0}$). The resonance peaks have been identifed as $\Omega_{Z}$,
$\Omega_{B}$, $\Omega_{X}$, and $\Omega_{B\pm}$.

We discuss first the periodic case shown in Fig. 6 (a), in which the probe is
aligned with the $z$-axis ($\theta_{p}=0$). This system was analyzed in detail
in Ref. [28], where it was confirmed that the spectral features denoted as
$\Omega_{Z}$ arise from probe-induced Raman transitions between adjacent
vibrational levels corresponding to oscillations along the $z$-axis in each
well, and the $\Omega_{B}$ features arise from Brillouin-like directed
transport along the $\pm x$-directions.

Indeed, the vibrational frequency in the $z$-direction is determined in the
system model by [28]

$\Omega_{z}=2(\mbox{cos}\,{\theta_{x}}+\mbox{cos}\,{\theta_{y}})\sqrt{2|\Delta_{0}^{\prime}|\,\omega_{r}},$
(13)

thus yielding $\Omega_{z}=2\pi(180$ kHz), which is in reasonable agreement
with the observed value of about 200 kHz. Moreover, Eq. (2) yields a
vibrational frequency $\Omega_{x}/2\pi=60$ kHz, which is very close to the
observed value for $\Omega_{B}/2\pi$. This $\Omega_{x}$-value equals 15.5
$\omega_{r}$, which in simulation units corresponds to 7.8, not far from the
7.3 value mentioned in Sec. III. The agreement between the predicted and
observed values for $\Omega_{x}$ and $\Omega_{B}$ is remarkable considering
that the theory assumed a $F_{g}=1/2\rightarrow F_{e}=3/2$ atom.

It is important to note that even though $\Omega_{B}$ coincides with
$\Omega_{x}$, the spectral features at $\Omega_{B}$ in Fig. 6(a) could not
arise from intrawell oscillatory motion in the $x$-direction, owing to the
fact that adjacent vibrational levels are of opposite parity, and hence the
overlap integral of the probe operator (i.e., the lattice-probe interference
term) between these two levels is zero (because the interference term goes as
$\vec{E}_{0}\cdot{\vec{E}^{*}}_{p}$ and is quadratic in $x$ for a probe that
propagates purely along $\hat{z}$; here $\vec{E}_{0}$ and $\vec{E}_{p}$ are
the lattice and probe electric field amplitudes, respectively, [28]).

In this space periodic case $\theta_{p}=0$, the three probe terms of (1)
reduces to a single perturbation propagating with phase velocity
$v_{0}=\delta_{p}/k_{0}$, like in the 1D model (4) with $k_{p}=k_{0}$, whose
numerical results are shown in Fig. 3. In agreement with the experimental
results of Fig. 6 (a), the theoretical model predicts a dominant peak at about
$\omega_{0}$, which is identified with $\Omega_{x}$. The secondary peaks
observed at other multiples of $\omega_{0}$ in Fig. 3 for certain values of
the driving amplitudes are absent in Fig. 6 (a), though.

Let us turn, then, our attention to the case where $k_{p}\neq k_{0}$, shown in
Fig. 6 (b). Equation (1) indicates that the probe driving produces three
perturbations propagating with phase velocities
$v_{1}=\delta_{p}/(k_{0}+k_{l}\sin\theta_{p})$,
$v_{2}=-\delta_{p}/(k_{0}-k_{l}\sin\theta_{p})$, and
$v_{3}=\delta_{p}/(k_{l}\sin\theta_{p})$. The first two perturbations are the
result of the interference between the $\hat{y}$-polarized probe and the
$\hat{y}$-polarized lattice beams $\vec{k}_{1}$, $\vec{k}_{2}$, as depicted in
Fig. 1. However, the third perturbation, being produced by the interference of
the probe with the $\hat{x}$-polarized lattice beams $\vec{k}_{3}$ and
$\vec{k}_{4}$, is not expected to show up in the transmission spectrum
depicted in Fig. 6 (b) because of Doppler broadening [9, 28]. Specifically,
the $z$-component of the motion makes the scattered field vary randomly, being
responsible for a wash-out of the spectral features related to the $x$-motion
[9] due to the perturbation propagating with velocity $v_{3}$. In other words,
the pump-probe measurements shown here are sensitive only to the
$U^{(b)}$-perturbation defined by (4) (i.e., case (b)), and are not sensitive
to the $U^{(a)}$-perturbation defined by (3) under case (a). As indicated in
Sec. III, the simulations in Figs. 3 \- 5 all pertain to case (b).

Case (b) refers to the first two perturbations in (1). Following the velocity
matching argument discussed in Sec. III, we expect a peak for the driving
frequency values where the intrinsic velocity $\pm v_{0}$ (11) matches the
velocity of the propagating perturbation, thus leading to two possible values
for the probe detuning $\delta_{p}$, located on either side of $\Omega_{x}$:

$\displaystyle\delta_{1}=\Omega_{x}\frac{\sin\theta_{x}+\sin\theta_{p}}{\sin\theta_{x}},$
(14)
$\displaystyle\delta_{2}=\Omega_{x}\frac{\sin\theta_{x}-\sin\theta_{p}}{\sin\theta_{x}},$
(15)

where we have identified the vibrational frequency $\omega_{0}=\Omega_{x}$.
The theoretical analysis of Sec. III also predicts a peak at a multiple of the
vibrational frequency, specially around the central value at $\Omega_{x}$, as
seen in Fig. 4.

In order to confirm the angle-dependencies in (14) and (15), we have carried
out several measurements with varying values of the probe angle $\theta_{p}$.
The results are plotted in Fig. 7.

Figure 7:  (a) Pump-probe spectra taken at different probe angles $\theta_{p}$
for the lattice in Fig. 1(a). The vertical lines demarcate the fixed
vibrational frequencies $\Omega_{z}$ and $\Omega_{x}$. At $\theta_{p}=0$, the
gray arrow denotes the Brillouin resonance $\Omega_{B}$ that coincides with
$\Omega_{x}$. At $\theta_{p}\neq 0$, the black and white arrows denote
Brillouin resonances $\Omega_{B+}$ and $\Omega_{B-}$, respectively. (b)
$\Omega_{B+}$ and $\Omega_{B-}$ depart from $\Omega_{x}$ in accordance with
Eqs. (14) and (15), respectively: The solid lines are the theoretical
predictions from the equations. The two black-circled datapoints in (b)
correspond to the spectrum in Fig. 6(b). It takes 10 ms to generate those
spectra, and each datapoint here is an average of at least five scans.

The solid lines in Fig. 7 (b) confirm the analytical predictions of Eqs. (14)
and (15), with $\delta_{1}$ and $\delta_{2}$ explaining the observed peaks
$\Omega_{B+}$ and $\Omega_{B-}$, respectively, that are seen in Fig. 7 (a) and
Fig. 6 (b). The quantitative agreement is remarkable, specially taking into
account that the theory is based on a 1D model with a simplified atomic
transition.

The theory also offers an explanation [9] for the fact that the amplitudes of
the $\Omega_{B-}$ resonance are seen to be smaller than the $\Omega_{B+}$
ones. Since the spatial period $2\pi/(k_{0}-k_{l}\sin\theta_{p})$ associated
for the $\Omega_{B-}$ motion is larger than that of $\Omega_{B+}$,
$2\pi/(k_{0}+k_{l}\sin\theta_{p})$, the sequence of half oscillations and
well-transfers for $\Omega_{B-}$ is more likely to be interrupted by random
photon recoils associated with the Sisyphus process, causing a larger damping
of the directed propagation. For the $\Omega_{B+}$ resonance,
$k_{p}=k_{0}+k_{l}\sin\theta_{p}$, yielding $k_{p}/k_{p0}$-values
corresponding to the four angles from $12.5^{0}$ to $20^{0}$ as 3.4, 3.6, 3.8,
and 4.1, respectively. Note that the first two values are close to the
$k_{p}/k_{p0}$-value of 3.5 in Fig. 5, which predicts peaks in the directed
propagation at $\omega_{0}$, $(\delta_{p})_{vm}\approx 1.5\omega_{0}$, and
$2\omega_{0}$. The peak at $\omega_{0}=\Omega_{x}\approx 2\pi\,(60\mbox{kHz})$
is certainly observed in Fig. 7 (a) and (b), as are single peaks at
approximately 105 kHz and 110 kHz for $\theta_{p}=12.5^{0}$ and $15^{0}$,
respectively. The experiment does not resolve each of these single peaks into
possible components at the numerical peaks, but the observed single-peak
values lie in between these, not too far from them. Thus the experimental
findings in Fig. 7 (a-b) and the numerical simulations in Fig. 5 do not
contradict each other.

Finally, figs. 7(a) and (b) show consistently a peak at
$\delta_{p}=\Omega_{x}$, regardless of the specific value of $\theta_{p}$.
While symmetry consideration for $\theta_{p}=0$ forbids the appearance of a
resonance in the spectrum due to localized vibrations about the well bottoms,
that is not the case for $\theta_{p}\neq 0$. However, as shown in Sec. III for
the 1D theoretical model, the observed peak may well also be due to directed
motion.

##  V Conclusions

We have studied the atomic waves generated in a dissipative optical lattice
under a weak beam that produces a directed current which travels perpendicular
to the direction of travel of the probe.

On the theoretical side, a minimal 1D model of the experimental setup is
studied in detail to elucidate the mechanisms of transport. An analytical
method based on a Fourier decomposition of the current is applied to study the
case when the driving potential perturbation has a spatial period which is the
same as that of the underlying lattice, a case of space periodic driving, and
when both periods are incommensurate, the regime of space quasiperiodic
driving.

It is numerically demonstrated that the transition between both regimes is
smooth, despite the fact that the expansion for one of the probe perturbations
is different in each regime. When the frequency of the probe is varied, the
current, and more specifically the mode amplitudes that contribute to the
current, show several peaks. One set is identified with a multiple of the
intrinsic frequency $\omega_{0}$, and another is associated with a velocity
matching mechanism, in which the velocity of the propagating modulation
matches up with the average velocity of the atom in its intrawell oscillation.
Both mechanisms are seen at play in both the space periodic and space
quasiperiodic regime, and even at play at the same time, since the shift due
to velocity matching is observed to take place through local maxima at a
fractional ratio of the intrinsic frequency $\omega_{0}$.

The pump-probe experiments confirm many of the above predictions. In the case
of space-periodic driving $k_{p}=k_{0}$, with a weak beam that is aligned with
the lattice symmetry axis, the probe transmission spectrum indeed reveals a
dominant peak at $\omega_{0}$, shown in Fig. 6(a), that corresponds to a
propagating Brillouin-like mode in a direction perpendicular to probe
propagation. These observations are borne out by the numerical predictions in
Fig. 3, although secondary peaks predicted at multiples and sub-multiples of
$\omega_{0}$ were not observed in this work. In the case of driving with the
weak beam incident at angle $\theta_{p}$ with the lattice symmetry axis, the
spectrum reveals a peak at $\omega_{0}$, and also additional peaks where the
velocity matching mechanism above is satisfied, as shown in Fig. 6(b). The
angle-dependence of these additional spectral features, which correspond to
two distinct propagating modulations with different spatial periods, is shown
in Fig. 7 to be in accordance with the analytical predictions. These findings
are not contradicted by the numerical predictions in Figs. 4 and 5, although
it was not experimentally possible to tease apart the contributions to the
observed resonances made by $(\delta_{p})_{vm}$ versus those made by multiples
and sub-multiples of $\omega_{0}$. We hope that the new understanding of
Brillouin transport modes gained here would pave the road toward ratcheting of
cold atoms confined in a weakly modulated optical lattice, along a precisely
predictable, arbitrary direction [7].

##  VI Acknowledgments

This work is supported by the Army Research Office under award/contract number
W911NF2110120 and the Ministerio de Ciencia e Innovación of Spain of Spain,
Grant No. PID2019-105316GB-I00 (DC). We thank the Instrumentation Laboratory
at Miami University for electronics and LabView support. We gratefully
acknowledge invaluable assistance in the lab from Ian Dilyard and Jordan
Churi.

##  Appendix A Further analytical results

The setup defined by case (a) in (3) requires different expressions for the
quasiperiodic ($k_{p}/k_{0}$ irrational) and periodic $k_{p}=k_{0}$ cases. The
calculation proceeds along the same lines sketched in [31]. The atomic state
symmetry produced by the probe for setup (3) is given by

$P_{-}[l,n,m]=(-1)^{n+l}P_{+}[l,n,m],$ (16)

which can be also written as

$P_{-}[l,n,m]=(-1)^{n+m}P_{+}[l,n,m],$ (17)

the latter being a useful expression in the special (periodic) case
$k_{p}=k_{0}$.

In the quasiperiodic case we find

$\displaystyle\langle
v\rangle_{\mathrm{(a)}}=\frac{m_{a}}{m_{a}F_{0}g_{1}-2D_{0}k_{0}}\Bigg{[}$
$\displaystyle-\frac{\text{Im}\left[P_{+}[0,1,0]\right]F_{0}\left(8g_{0}^{2}-4g_{2}^{2}/3+F_{0}k_{0}/(2m_{a})\right)}{k_{0}}$
$\displaystyle+\frac{\text{Im}\left[P_{+}[0,2,0]\right]F_{0}\left(-4g_{0}g_{1}{-8g_{1}g_{2}/3}+F_{0}k_{0}/m_{a}\right)}{k_{0}}$
$\displaystyle+\frac{\text{Im}\left[P_{+}[0,3,0]\right]F_{0}\left(-16g_{0}g_{1}/3-2g_{2}^{2}-3F_{0}k_{0}/(2m_{a})\right)}{k_{0}}$
$\displaystyle+\frac{\text{Im}\left[P_{+}[0,4,0]\right]F_{0}\left(-2g_{1}g_{2}/3+F_{0}k_{0}/(2m_{a})\right)}{k_{0}}$
$\displaystyle-\frac{\text{Im}\left[P_{+}[0,5,0]\right]F_{0}2g_{2}^{2}}{3k_{0}}+\frac{\text{Im}\left[e^{i\phi_{p}}P_{+}[1,-4,1]\right](2F_{p}g_{2}^{2}k_{0})}{(2k_{0}-k_{p})k_{p}}$
$\displaystyle+\frac{\text{Im}\left[e^{i\phi_{p}}P_{+}[1,-3,1]\right]2F_{p}g_{1}g_{2}k_{0}}{(2k_{0}-k_{p})k_{p}}$
$\displaystyle+\frac{\text{Im}\left[e^{i\phi_{p}}P_{+}[1,-2,1]\right]}{(2k_{0}-k_{p})k_{p}^{2}}[8F_{p}g_{0}g_{2}(-k_{0}^{2}+k_{0}k_{p})$
$\displaystyle+F_{0}F_{p}(2k_{0}^{2}k_{p}-2k_{0}k_{p}^{2}+k_{p}^{3}/2)/m_{a}]$
$\displaystyle-\frac{\text{Re}\left[e^{i\phi_{p}}P_{+}[1,-2,1]\right]4F_{p}g_{2}k_{0}(k_{0}-k_{p})\delta_{p}}{(2k_{0}-k_{p})k_{p}^{2}}$
$\displaystyle+\frac{\text{Im}\left[e^{i\phi_{p}}P_{+}[1,-1,1]\right]}{(2k_{0}-k_{p})k_{p}^{2}}F_{p}\Bigg{[}4g_{0}g_{1}k_{0}(-2k_{0}+k_{p})$
$\displaystyle+2g_{1}g_{2}k_{0}k_{p}+\frac{F_{0}}{m_{a}}(-2k_{0}^{2}k_{p}+3k_{0}k_{p}^{2}-k_{p}^{3})\Bigg{]}$
$\displaystyle{\color[rgb]{0,0,1}-\frac{\text{Re}\left[e^{i\phi_{p}}P_{+}[1,-1,1]\right]2F_{p}g_{1}k_{0}(k_{0}-2k_{p})\delta_{p}}{(k_{0}-k_{p})k_{p}^{2}}}$
$\displaystyle-\frac{\text{Im}\left[e^{i\phi_{p}}P_{+}[1,0,1]\right]}{(2k_{0}-k_{p})k_{p}^{2}(2k_{0}+k_{p})}F_{p}[g_{0}^{2}(32k_{0}^{3}-8k_{0}k_{p}^{2})$
$\displaystyle-4g_{2}^{2}k_{0}k_{p}^{2}+\delta_{p}^{2}(-8k_{0}^{3}+2k_{0}k_{p}^{2})]$
$\displaystyle-\frac{\text{Re}\left[e^{i\phi_{p}}P_{+}[1,0,1]\right]8F_{p}g_{0}k_{0}\delta_{p}}{k_{p}^{2}}$
$\displaystyle+\frac{\text{Im}\left[e^{i\phi_{p}}P_{+}[1,1,1]\right]}{(2k_{0}+k_{p})k_{p}^{2}}F_{p}[-4g_{0}g_{1}k_{0}(2k_{0}+k_{p})$
$\displaystyle-2g_{1}g_{2}k_{0}k_{p}+F_{0}(2k_{0}^{2}k_{p}+3k_{0}k_{p}^{2}+k_{p}^{3})/m_{a}]$
$\displaystyle-\frac{\text{Re}\left[e^{i\phi_{p}}P_{+}[1,1,1]\right]2F_{p}g_{1}k_{0}(k_{0}+k_{p})\delta_{p}}{(k_{0}+k_{p})k_{p}^{2}}$
$\displaystyle+\frac{\text{Im}\left[e^{i\phi_{p}}P_{+}[1,2,1]\right]}{(2k_{0}+k_{p})k_{p}^{2}}[8F_{p}g_{0}g_{2}(-k_{0}^{2}-k_{0}k_{p})$
$\displaystyle+F_{0}F_{p}(-2k_{0}^{2}k_{p}-2k_{0}k_{p}^{2}-k_{p}^{3}/2)/m_{a}]$
$\displaystyle-\frac{\text{Re}\left[e^{i\phi_{p}}P_{+}[1,2,1]\right]4F_{p}g_{2}k_{0}(k_{0}+k_{p})\delta_{p}}{(2k_{0}+k_{p})k_{p}^{2}}$
$\displaystyle-\frac{\text{Im}\left[e^{i\phi_{p}}P_{+}[1,3,1]\right]2F_{p}g_{1}g_{2}k_{0}}{(2k_{0}+k_{p})k_{p}}$
$\displaystyle-\frac{\text{Im}\left[e^{i\phi_{p}}P_{+}[1,4,1]\right](2F_{p}g_{2}^{2}k_{0})}{(2k_{0}+k_{p})k_{p}}$
$\displaystyle+\frac{\text{Im}\left[e^{i2\phi_{p}}P_{+}[2,0,2]\right]F_{p}^{2}k_{0}}{m_{a}k_{p}}\Bigg{]}.$

Equation (A) is not valid when the ratio $k_{p}/k_{0}$ is a rational number.
In this case there are resonances which require a special derivation [31]. An
indication of this fact is that the coefficient associated with the mode
amplitude $\mathrm{Re}[P_{+}[1,-1,1]]$ goes to infinity in the limit
$k_{p}\rightarrow k_{0}$, whereas the actual coefficient when computed
directly in the case $k_{p}=k_{0}$ remains finite, as shown in Eq. (A).

The full expansion in the periodic case $k_{p}=k_{0}$ is given by

$\displaystyle\langle
v\rangle_{\mathrm{(a)}}=\frac{m_{a}}{m_{a}F_{0}g_{1}-2D_{0}k_{0}}\Bigg{[}$
$\displaystyle{\color[rgb]{0,0,1}-\frac{\text{Re}\left[P_{+}[0,1,0]\right]F_{0}\left(F_{p}^{2}g_{1}\right)}{m_{a}\delta_{p}}}$
$\displaystyle-\frac{\text{Im}\left[P_{+}[0,1,0]\right]F_{0}\left(8g_{0}^{2}-4g_{2}^{2}/3+F_{0}k_{0}/(2m_{a})\right)}{k_{0}}$
$\displaystyle+\frac{\text{Im}\left[P_{+}[0,2,0]\right]F_{0}\left(-4g_{0}g_{1}{-8g_{1}g_{2}/3}+F_{0}k_{0}/m_{a}\right)}{k_{0}}$
$\displaystyle+\frac{\text{Im}\left[P_{+}[0,3,0]\right]F_{0}\left(-16g_{0}g_{1}/3-2g_{2}^{2}-3F_{0}k_{0}/(2m_{a})\right)}{k_{0}}$
$\displaystyle+\frac{\text{Im}\left[P_{+}[0,4,0]\right]F_{0}\left(-2g_{1}g_{2}/3+F_{0}k_{0}/(2m_{a})\right)}{k_{0}}$
$\displaystyle-\frac{\text{Im}\left[P_{+}[0,5,0]\right]F_{0}2g_{2}^{2}}{3k_{0}}+\frac{\text{Im}\left[e^{i\phi_{p}}P_{+}[1,-4,1]\right](2F_{p}g_{2}^{2})}{k_{0}}$
$\displaystyle{\color[rgb]{0,0,1}+\frac{\text{Re}\left[e^{i\phi_{p}}P_{+}[1,-3,1]\right](F_{0}F_{p}g_{1})}{m_{a}\delta_{p}}}$
$\displaystyle+\frac{\text{Im}\left[e^{i\phi_{p}}P_{+}[1,-3,1]\right]2F_{p}g_{1}g_{2}}{k_{0}}$
$\displaystyle+\frac{\text{Im}\left[e^{i\phi_{p}}P_{+}[1,-2,1]\right]F_{0}F_{p}}{2m_{a}}$
$\displaystyle{\color[rgb]{0,0,1}-\frac{\text{Re}\left[e^{i\phi_{p}}P_{+}[1,-2,1]\right]F_{0}F_{p}g_{1}}{m_{a}\delta_{p}}}$
$\displaystyle+\frac{\text{Im}\left[e^{i\phi_{p}}P_{+}[1,-1,1]\right]}{k_{0}}F_{p}(-4g_{0}g_{1}+2g_{1}g_{2})$
$\displaystyle{\color[rgb]{0,0,1}-\frac{\text{Re}\left[e^{i\phi_{p}}P_{+}[1,-1,1]\right]2F_{p}g_{1}\delta_{p}}{k_{0}}}$
$\displaystyle-\frac{\text{Im}\left[e^{i\phi_{p}}P_{+}[1,0,1]\right]}{3k_{0}}2F_{p}(12g_{0}^{2}-2g_{2}^{2}-3\delta_{p}^{2})$
$\displaystyle-\frac{\text{Re}\left[e^{i\phi_{p}}P_{+}[1,0,1]\right]({\color[rgb]{0,0,1}-F_{0}F_{p}g_{1}k_{0}/\delta_{p}}+8F_{p}g_{0}\delta_{p})}{k_{0}}$
$\displaystyle+\frac{\text{Im}\left[e^{i\phi_{p}}P_{+}[1,1,1]\right]}{(3k_{0}}F_{p}[-12g_{0}g_{1}-2g_{1}g_{2}+6F_{0}k_{0}/m_{a}]$
$\displaystyle-\frac{\text{Re}\left[e^{i\phi_{p}}P_{+}[1,1,1]\right]{\color[rgb]{0,0,1}(F_{p}F_{0}g_{1}k_{0}/\delta_{p}}+3F_{p}g_{1}\delta_{p})}{k_{0}}$
$\displaystyle+\frac{\text{Im}\left[e^{i\phi_{p}}P_{+}[1,2,1]\right]}{(6k_{0}}F_{p}(-32g_{0}g_{2}-9F_{0}k_{0}/m_{a})$
$\displaystyle-\frac{\text{Re}\left[e^{i\phi_{p}}P_{+}[1,2,1]\right]8F_{p}g_{2}\delta_{p}}{3k_{0}}$
$\displaystyle-\frac{\text{Im}\left[e^{i\phi_{p}}P_{+}[1,3,1]\right]2F_{p}g_{1}g_{2}}{3k_{0}}$
$\displaystyle-\frac{\text{Im}\left[e^{i\phi_{p}}P_{+}[1,4,1]\right](2F_{p}g_{2}^{2})}{3k_{0}}$
$\displaystyle{\color[rgb]{0,0,1}+\frac{\text{Im}\left[e^{i2\phi_{p}}P_{+}[2,-1,2]\right]F_{p}^{2}g_{1}}{\delta_{p}}}$
$\displaystyle+\frac{\text{Im}\left[e^{i2\phi_{p}}P_{+}[2,0,2]\right]F_{p}^{2}}{m_{a}}\Bigg{]}\quad(k_{p}=k_{0}).$

Terms in (A), which are not directly obtained from (A) after taking the limit
$k_{p}\rightarrow k_{0}$, have been highlighted in blue. Note there are also
terms in the quasiperiodic case (A) which are singular when $k_{p}\rightarrow
2k_{0}$. But again, this singularity is removed in the exact limit
$k_{p}=2k_{0}$.

## References

  * Heinrich  _et al._ [2021] A. Heinrich, W. D. Oliver, L. M. K. Vandersypen, A. Ardavan, R. Sessoli, D. Loss, A. B. Jayich, J. Fernandez-Rossier, A. Laucht, and A. Morello, Quantum-coherent nanoscience, Nat. Nanotechnol. 16, 1318 (2021).
  * Wang  _et al._ [2018] Y. Wang, S. Subhankar, P. Bienias, M. Łącki, T.-C. Tsui, M. A. Baranov, A. V. Gorshkova, P. Zoller, J. V. Porto, and S. Rolston, Dark state optical lattice with sub-wavelength spatial structure, Phys. Rev. Lett. 120, 083601 (2018).
  * Łącki  _et al._ [2016] M. Łącki, M. A. Baranov, H. Pichler, , and P. Zoller, Nanoscale dark state optical potentials for cold atoms, Phys. Rev. Lett. 117, 233001 (2016).
  * Miles  _et al._ [2013] J. A. Miles, Z. J. Simmons, and D. D. Yavuz, Subwavelength localization of atomic excitation using electromagnetically induced transparency, Phys. Rev. X 3, 031014 (2013).
  * Zemanek  _et al._ [2020] P. Zemanek, G. Volpe, A. Jonas, and O. Brzobohaty, Perspective on light-induced transport of particles: from optical forces to phoretic motion, Adv. Opt. Photonics 11, 577 (2020).
  * Cubero and Renzoni [2016] D. Cubero and F. Renzoni,  _Brownian Ratchets: From Statistical Physics to Bio and Nano-motors_ (Cambridge University Press, Cambridge, 2016).
  * Cubero and Renzoni [2012] D. Cubero and F. Renzoni, Control of transport in two-dimensional systems via dynamical decoupling of degrees of freedom with quasiperiodic driving fields, Phys. Rev. E 86, 056201 (2012).
  * Courtois  _et al._ [1996] J. Y. Courtois, S. Guibal, D. R. Meacher, P. Verkerk, and G. Grynberg, Propagating elementary excitation in a dilute optical lattice, Phys. Rev. Lett. 77, 40 (1996).
  * Jurczak  _et al._ [1998] C. Jurczak, J.-Y. Courtois, B. Desruelle, C. Westbrook, , and A. Aspect, Spontaneous light scattering from propagating density fluctuations in an optical lattice, Eur. Phys. J. D 1, 53 (1998).
  * Mennerat-Robilliard  _et al._ [1999] C. Mennerat-Robilliard, D. Lucas, S. Guibal, J. Tabosa, C. Jurczak, J.-Y. Courtois, and G. Grynberg, Ratchet for cold rubidium atoms: The asymmetric optical lattice, Phys. Rev. Lett. 82, 851 (1999).
  * Grynberg and Robilliard [2001] G. Grynberg and C. Robilliard, Cold atoms in dissipative optical lattices, Phys. Rep. 355, 335 (2001).
  * Schiavoni  _et al._ [2002a] M. Schiavoni, L. Sanchez-Palencia, F. Carminati, F. Renzoni, and G. Grynberg, Dark propagation modes in optical lattices, Phys. Rev. A 66, 053821 (2002a).
  * Sanchez-Palencia  _et al._ [2002] L. Sanchez-Palencia, F.-R. Carminati, M. Schiavoni, F. Renzoni, and G. Grynberg, Brillouin propagation modes in optical lattices: Interpretation in terms of nonconventional stochastic resonance, Phys. Rev. Lett. 88, 133903 (2002).
  * Carminati  _et al._ [2003] F.-R. Carminati, M. Schiavoni, Y. Todorov, F. Renzoni, and G. Grynberg, Pump-probe spectroscopy of atoms cooled in a 3d lin-perp-lin optical lattice, Eur. Phys. J. D 22, 311 (2003).
  * Sanchez-Palencia and Grynberg [2003] L. Sanchez-Palencia and G. Grynberg, Synchronization of hamiltonian motion and dissipative effects in optical lattices: Evidence for a stochastic resonance, Phys. Rev. A 68, 023404 (2003).
  * Schiavoni  _et al._ [2003] M. Schiavoni, L. Sanchez-Palencia, F. Renzoni, and G. Grynberg, Phase control of directed diffusion in a symmetric optical lattice, Phys. Rev. Lett. 90, 094101 (2003).
  * Jones  _et al._ [2004] P. H. Jones, M. Goonasekara, and F. Renzoni, Rectifying fluctuations in an optical lattice, Phys. Rev. Lett. 93, 073904 (2004).
  * Gommers  _et al._ [2005] R. Gommers, S. Bergamini, and F. Renzoni, Dissipation-induced symmetry breaking in a driven optical lattice, Phys. Rev. Lett. 95, 073003 (2005).
  * Gommers  _et al._ [2006] R. Gommers, S. Denisov, and F. Renzoni, Quasiperiodically driven ratchets for cold atoms, Phys. Rev. Lett. 96, 240604 (2006).
  * Gommers  _et al._ [2007] R. Gommers, M. Brown, and F. Renzoni, Symmetry and transport in a cold atom ratchet with multifrequency driving,, Phys. Rev. A 75, 053406 (2007).
  * Gommers  _et al._ [2008] R. Gommers, V. L. andM. Brown, and F. Renzoni, Gating ratchet for cold atoms, Phys. Rev. Lett. 100, 040603 (2008).
  * Lebedev and Renzoni [2009] V. Lebedev and F. Renzoni, Two-dimensional rocking ratchet for cold atoms, Phys. Rev. A 80, 023422 (2009).
  * Cubero  _et al._ [2010] D. Cubero, V. Lebedev, and F. Renzoni, Current reversals in a rocking ratchet: dynamical vs symmetry-breaking mechanisms, Phys. Rev. E 82, 041116 (2010).
  * Wickenbrock  _et al._ [2011] A. Wickenbrock, D. Cubero, N. A. A. Wahab, P. Phoonthong, and F. Renzoni, Current reversals in a rocking ratchet: The frequency domain, Phys. Rev. E 84, 021127 (2011).
  * Hagman  _et al._ [2011] H. Hagman, M. Zelan, C. M. Dion, and A. Kastberg, Directed transport with real-time steering and drifts along predesigned paths using a brownian motor,, Phys. Rev. A 83, 020101(R) (2011).
  * Zelan  _et al._ [2011] M. Zelan, H. Hagman, G. Labaigt, S. Jonsell, , and C. M. Dion, Experimental measurement of efficiency and transport coherence of a cold-atom brownian motor in optical lattices, Phys. Rev. A 83, 020102(R) (2011).
  * Wickenbrock  _et al._ [2012] A. Wickenbrock, P. C. Holz, N. A. A. Wahab, P. Phoonthong, D. Cubero, and F. Renzoni, Vibrational mechanics in an optical lattice: Controlling transport via potential renormalization, Phys. Rev. Lett. 108, 020603 (2012).
  * Staron  _et al._ [2022] A. Staron, K. Jiang, C. Scoggins, D. Wingert, D. Cubero, and S. Bali, Observation of stochastic resonance in directed propagation of cold atoms, Phys. Rev. Research 4, 043211 (2022).
  * [29] H. J. Metcalf and P. van der Straten,  _Cooling and Trapping_.
  * Hodapp  _et al._ [1995] T. W. Hodapp, C. Gerz, C. Furtlehner, C. I. Westbrook, W. D. Phillips, and J. Dalibard, Three-dimensional spatial diffusion in optical molasses, Appl. Phys. B 60, 135 (1995).
  * Cubero [2023] D. Cubero, Brillouin propagation modes of cold atoms undergoing sisyphus cooling, Phys. Rev. E 107, 034102 (2023).
  * Lutz and Renzoni [2013] E. Lutz and F. Renzoni, Beyond boltzmann–gibbs statistical mechanics in optical lattices, Nature Physics 9, 615 (2013).
  * Cubero  _et al._ [2014] D. Cubero, J. Casado-Pascual, and F. Renzoni, Irrationality and quasiperiodicity in driven nonlinear systems, Phys. Rev. Lett. 112, 174102 (2014).
  * Cubero and Renzoni [2018] D. Cubero and F. Renzoni, Asymptotic theory of quasiperiodically driven quantum systems, Phys. Rev. E 97, 062139 (2018).
  * Cubero  _et al._ [2018] D. Cubero, G. Robb, and F. Renzoni, Avoided crossing and sub-fourier-sensitivity in driven quantum systems, Phys. Rev. Lett. 121, 213904 (2018).
  * Petsas  _et al._ [1999] K. Petsas, G. Grynberg, and J.-Y. Courtois, Semiclassical monte carlo approaches for realistic atoms in optical lattices, Eur. Phys. J. D 6, 29 (1999).
  * Schiavoni  _et al._ [2002b] M. Schiavoni, F.-R. Carminati, L. Sanchez-Palencia, F. Renzoni, and G. Grynberg, Stochastic resonance in periodic potentials: Realization in a dissipative optical lattice, Europhys. Lett. 59, 493 (2002b).
  * Note [1] This definition of the probe detuning $\delta_{p}$ corrects a sign error in Ref. [28].
  * Casado-Pascual  _et al._ [2013] J. Casado-Pascual, D. Cubero, and F. Renzoni, Universal asymptotic behavior in nonlinear systems driven by a two-frequency forcing, Phys. Rev. E 88, 062919 (2013).
  * Kloeden and Platen [1992] P. Kloeden and E. Platen,  _Numerical Solution of Stochastic Differential Equations_ (Springer, New York, 1992).


2023

1]Department of Electrical and Computer Engineering, University of Pittsburgh,
3700 O’Hara St, Pittsburgh, 15213, PA, USA

2]School of Nursing, University of Pittsburgh, 3500 Victoria Street,
Pittsburgh, 15261, PA, USA

# Prediction of Post-Operative Renal and Pulmonary Complications Using
Transformers

Reza Shirkavand<EMAIL_ADDRESS>Fei Zhang<EMAIL_ADDRESS>Heng Huang
<EMAIL_ADDRESS>[  [

###### Abstract

Postoperative complications pose a significant challenge in the healthcare
industry, resulting in elevated healthcare expenses and prolonged hospital
stays, and in rare instances, patient mortality. To improve patient outcomes
and reduce healthcare costs, healthcare providers rely on various
perioperative risk scores to guide clinical decisions and prioritize care. In
recent years, machine learning techniques have shown promise in predicting
postoperative complications and fatality, with deep learning models achieving
remarkable success in healthcare applications. However, research on the
application of deep learning models to intra-operative anesthesia management
data is limited. In this paper, we evaluate the performance of transformer-
based models in predicting postoperative acute renal failure, postoperative
pulmonary complications, and postoperative in-hospital mortality. We compare
our method’s performance with state-of-the-art tabular data prediction models,
including gradient boosting trees and sequential attention models, on a
clinical dataset. Our results demonstrate that transformer-based models can
achieve superior performance in predicting postoperative complications and
outperform traditional machine learning models. This work highlights the
potential of deep learning techniques, specifically transformer-based models,
in revolutionizing the healthcare industry’s approach to postoperative care.

###### keywords:

Transformers, Postoperative Renal Complications, Postoperative Pulmonary
Complications

##  1 Main

Post-operative complications present a significant challenge in the healthcare
industry, resulting in elevated healthcare expenses, prolonged hospital stays,
and, in rare instances pearse2006mortalityrate1 ; pearse2012mortalityrate2 ,
patient mortality. With approximately 230 million surgeries conducted globally
each year weiser2008surgnum , it is crucial for healthcare providers to
predict potential complications and take preventive measures to improve
patient outcomes and reduce healthcare costs.

Various perioperative risk scores are used to guide clinical decisions and
prioritize care for patients undergoing surgery. Clinical risk scores help
plan clinical management and allocate resources, while administrative risk
scores help compare hospitals. The American Society of Anesthesiologists
Physical Status Classification (ASA) dripps1963asa  and the Surgical Apgar
gawande2007apgar  score are commonly used in the perioperative setting, but
the latter has limitations such as subjectivity and variability. Logistic
regression is used to develop new scoring systems le2016pospom , and there are
other recent perioperative administrative risk scores such as the Risk
Stratification Index (RSI) sessler2010rsi  and the Risk Quantification Index
(RQI) dalton2011rqi .

Healthcare providers can use predictive models to identify high-risk patients
and implement preemptive measures to prevent fatal post-operative
complications. Machine learning techniques can be applied to analyze vast
amounts of patient data, such as medical history, vital signs, laboratory
values, and radiological findings, to identify patients at high risk of
developing complications after surgery.

Deep Learning (DL) is a cutting-edge method for data classification and
prediction, and is seen as the most promising solution to problems of data
complexity. DL has achieved remarkable success in industry, as evidenced by
recent studies nvidia2018 . In healthcare, DL has been used for computational
phenotyping lasko2013computationalpheno ; che2015deepcomppheno , risk
prediction cheng2016risk ; miotto2016deep , medical imaging analysis
esteva2017dermatologist , and natural language processing muneeb2015evaluating
.

Although there have been significant achievements in deep learning models for
other areas of healthcare, research on the application of deep learning models
to intra-operative anesthesia management data is limited. This type of data
has many advantages over typical EMR data(i.e. low-resolution AIMS data and
high-resolution data directly from intra-operative physiological monitors),
including data heterogeneity, longitudinal irregularity, inherent noise, and
incomplete nature. However, the research in this area has mainly targeted
static data such as demographic data, medical history, and post-operative
complications muneeb2015evaluating ; fritz2019deep . In fact, there is a
scarcity of research on the application of deep learning models to dynamic
intra-operative physiological data. In a real-time post-cardiosurgical care
setting, meyer2018rnncomplication  utilized recurrent neural networks to
predict severe complications like mortality, renal failure requiring renal
replacement therapy, and postoperative bleeding leading to operative revision.
In another study xue2021mlpulmonary , logistic regression and decision tree
variants were used to predict postoperative pulmonary complications.
Similarly, lee2018mortalityprediction  employed feed-forward networks to
predict post-operative mortality in critically ill patients.

Traditional machine learning methods, such as Gradient Boosting Decision Trees
(GBDT) and its variants, including XGBoost chen2016xgboost , LightGBM
ke2017lightgbm , and CatBoost prokhorenkova2018catboost , have dominated the
tabular data domain. Nonetheless, recent research has shown that deep learning
(DL) methods can also achieve impressive results. For example, TabNet
arik2021tabnet  introduced sequential attention to reason based on the most
salient features, thereby providing interpretability. Neural Oblivious
Decision Ensembles (NODE) popov2019node  incorporate hierarchical
representation learning into the network. GrowNet badirli2020grownet  utilizes
gradient boosting frameworks to construct a complex deep neural network from
simple shallow components. Self Normalizing Neural Networks (SNN)
klambauer2017snn  apply strong regularization and encourage zero mean and unit
variance activations to facilitate high-level feature learning.

Transformers vaswani2017transformer , initially introduced in the field of
Natural Language Processing (NLP), are a type of neural network architecture
specifically designed to handle sequential data. While they have been highly
effective in NLP tasks, such as machine translation and text summarization,
they have also been shown to be well-suited for tabular data analysis
gorishniy2021rdtl .

In this work, we aim to evaluate the performance of transformer-based models
in predicting postoperative acute renal failure (PO-ARF), postoperative
pulmonary complications (PPC), and postoperative in-hospital mortality from a
tabular dataset of dynamic intra-operative anesthesia management variables. We
compare the performance of our method with state-of-the-art prediction models
suited to the dataset, including gradient boosting trees, and assess their
ability to accurately predict the aforementioned complications using clinical
data.

##  2 Results

In this section, we present the results of our experiments comparing the
performance of our transformer model to gradient boosting decision tree (GBDT)
variants. The goal of our tasks were to predict postoperative pulmonary
complications (PPC), the postoperative acute renal failure (PO-ARF), and
postoperative in-hospital mortality.

Table 1 and Table 2 show the accuracy, precision, recall, and F1 score of the
models across 5-fold cross validation for both post-operative complication and
post-operative in-hospital mortality prediction tasks. The transformer-based
model outperforms all GBDT variants and deep baselines in all evaluation
metrics.

Table 1: Comparison of performance of baselines with the transformer model.
Values correspond to the mean and standard deviation of each classification
metric across 5-fold cross-validation of predicting post-operative pulmonary
and renal complication

Method | Accuracy | Precision | Recall | F1-score
---|---|---|---|---
PPC | PO-ARF | PPC | PO-ARF | PPC | PO-ARF | PPC | PO-ARF
SVM | 0.803 ± 0.020 | 0.914 ± 003 | 0.779 ± 0.014 | 0.755 ± 0.007 | 0.802 ± 0.009 | 0.857 ± 0.017 | 0.785 ± 0.017 | 0.793 ± 0.008
Logistic Regression | 0.796 ± 0.024 | 0.896 ± 0.002 | 0.775 ± 0.015 | 0.730 ± 0.003 | 0.803 ± 0.011 | 0.882 ± 0.009 | 0.780 ± 0.020 | 0.777 ± 0.004
Random Forest | 0.813 ± 0.011 | 0.933 ± 0.007 | 0.807 ± 0.017 | 0.855 ± 0.025 | 0.748 ± 0.017 | 0.708 ± 0.033 | 0.765 ± 0.015 | 0.756 ± 0.034
KNN | 0.787 ± 0.015 | 0.916 ± 0.004 | 0.767 ± 0.022 | 0.824 ± 0.033 | 0.724 ± 0.014 | 0.596 ± 0.018 | 0.737 ± 0.016 | 0.634 ± 0.026
XGboost | 0.817 ± 0.023 | 0.944 ± 0.009 | 0.797 ± 0.035 | 0.848 ± 0.025 | 0.751 ± 0.019 | 0.818 ± 0.030 | 0.766 ± 0.023 | 0.832 ± 0.028
LightGBM | 0.822 ± 0.021 | 0.944 ± 0.009 | 0.794 ± 0.026 | 0.846 ± 0.023 | 0.775 ± 0.013 | 0.819 ± 0.036 | 0.782 ± 0.019 | 0.832 ± 0.031
Catboost | 0.822 ± 0.022 | 0.945 ± 0.009 | 0.797 ± 0.028 | 0.849 ± 0.022 | 0.772 ± 0.014 | 0.819 ± 0.039 | 0.780 ± 0.019 | 0.833 ± 0.031
MLP | 0.831 ± 0.015 | 0.891 ± 0.006 | 0.798 ± 0.016 | 0.727 ± 0.010 | 0.807 ± 0.007 | 0.898 ± 0.004 | 0.802 ± 0.013 | 0.776 ± 010
SNN | 0.837 ± 0.002 | 0.900 ± 0.008 | 0.805 ± 0.003 | 0.737 ± 0.008 | 0.810 ± 0.003 | 0.894 ± 0.006 | 0.807 ± 0.003 | 0.785 ± 0.005
Node | 0.830 ± 0.003 | 0.922 ± 0.014 | 0.796 ± 0.004 | 0.778 ± 0.026 | 0.807 ± 0.005 | 0.879 ± 0.024 | 0.801± 0.004 | 0.814 ± 0.016
Ours | 0.865 ± 0.005 | 0.967 ± 0.007 | 0.830 ± 0.005 | 0.848 ± 0.009 | 0.839 ± 0.004 | 0.917 ± 0.008 | 0.834 ± 0.006 | 0.885 ± 0.007

Table 2: Comparison of performance of baselines with the transformer model.
Values correspond to the mean and standard deviation of each classification
metric across 5-fold cross-validation of predicting postoperative in-hospital
mortality

Method | Accuracy | Precision | Recall | F1-score
---|---|---|---|---
SVM | 0.981 ± 0.001 | 0.748 ± 0.019 | 0.733 ± 0.018 | 0.740 ± 0.013
Logistic Regression | 0.926 ± 0.008 | 0.597 ± 0.008 | 0.906 ± 0.006 | 0.641 ± 0.012
Random Forest | 0.982 ± 0.001 | 0.875 ± 0.056 | 0.559 ± 0.004 | 0.598 ± 0.006
KNN | 0.981 ± 0.001 | 0.804 ± 0.045 | 0.520 ± 0.005 | 0.533 ± 0.009
XGboostchen2016xgboost  | 0.990 ± 0.001 | 0.923 ± 0.010 | 0.781 ± 0.024 | 0.836 ± 0.019
Light GBMke2017lightgbm  | 0.990 ± 0.001 | 0.909 ± 0.006 | 0.789 ± 0.021 | 0.838 ± 0.016
Catboostprokhorenkova2018catboost  | 0.990 ± 0.001 | 0.919 ± 0.012 | 0.792 ± 0.014 | 0.844 ± 0.009
MLP | 0.984 ± 0.002 | 0.771 ± 0.034 | 0.856 ± 0.016 | 0.806 ± 0.017
SNNklambauer2017snn  | 0.983 ± 0.004 | 0.766 ± 0.041 | 0.867 ± 0.029 | 0.804 ± 0.022
Nodepopov2019node  | 0.990 ± 0.001 | 0.881 ± 0.001 | 0.800 ± 0.037 | 0.834± 0.022
Ours | 0.995 ± 0.005 | 0.930 ± 0.006 | 0.902 ± 0.004 | 0.918 ± 0.005

We further analyzed the performance of the models by computing the precision
recall curve. Figure 1 shows the Precision-Recall curve of all methods. The
transformer-based model achieved the highest area under the curve (AUPRC),
indicating that it has a better ability to accurately predict the
complications. The error bars of AUPRC is displayed in Figure 2

(a) Pulmonary complications

(b) Renal complications

(c) In-hospital mortality

Figure 1: Precision Recall Curve of baselines as well as the transformer
architecture for the prediction of pulmonary, renal complications, and in-
hospital mortality

(a) Pulmonary complications

(b) Renal complications

(c) In-hospital mortality

Figure 2: Mean and standard deviation of area under precision recall curve
(AUPRC) for all three task.

We attribute the substantial superiority of our method to several factors:

  1. 1.

Ability to learn complex relationships: Transformers are powerful models that
are capable of learning complex relationships between variables. They can
capture long-range dependencies in the data and model non-linear interactions
between features, which may not be possible with simpler models like decision
trees and MLPs.

  2. 2.

Feature interactions: Tabular data often contains many features, and the
interactions between these features can be crucial for making accurate
predictions. Transformers are designed to handle feature interactions, and can
learn to represent features in a way that captures their relationships with
other features.

  3. 3.

Robustness to missing data: Tabular data may have missing values (as was in
our case where missing lab values where imputed with zeros), which can be a
challenge for some models. Transformers can handle missing data by leveraging
the attention mechanism to focus on relevant features and ignoring irrelevant
ones.

It is worth noting that transformers have been shown gorishniy2021rdtl  to be
a more universal model for structured data like our dataset.

As shown in Figure 3, the top 50 most important features of the dataset are
plotted using a bar chart. The features were determined using permutation
feature importance. Permutation feature importance works by randomly permuting
the values of a single feature and measuring the decrease in model performance
(e.g. F1 score) on a validation set. The larger the decrease in performance,
the more important the feature is considered to be.

To calculate permutation feature importance for a transformer model, we:

  1. 1.

Train the transformer model on the dataset and evaluate the performance of the
trained model on a separate validation set using the F1 score.

  2. 2.

For each feature in the validation set, randomly permute its values and re-
evaluate the F1 score. The difference between the original F1 score and the
permuted F1 score is the feature importance score for that feature.

  3. 3.

Repeat step 3 for all features in the validation set and rank them based on
their importance scores.

Figure 3: Feature Importance Plot: The graph shows the relative importance of
features in predicting the target variable, as determined by permutation
feature importance.

##  3 Discussion

In this study, we conducted a retrospective analysis of the relationship
between perioperative anesthesia management and patient outcomes for major
surgeries under general anesthesia. We found that transformer-based models can
effectively predict in-hospital morbidity and mortality based on automatically
extractable and objective perioperative data.

Major surgery with general anesthesia can cause renal hypoperfusion,
insufficient oxygen delivery, and eventually postoperative renal failure.
Additionally, intraoperative neuromuscular blocking drugs and ventilation can
change respiratory drive and muscle function, increasing the risk of
postoperative pulmonary complications. These complications are associated with
increased in-hospital mortality and are difficult to predict due to their
heterogeneity. Risk factors for postoperative renal and pulmonary
complications are numerous, thus anesthesia providers must recognize
modifiable and non-modifiable risk factors and optimize patient care
accordingly. Many validated risk prediction models are described.

We compared the transformer-based model’s performance to other validated risk
prediction models, demonstrating that the transformer-based model outperforms
other models due to its ability to capture inter-dependencies in the extracted
variables important for predicting postoperative complications. Our findings
have significant implications for improving postoperative care by providing
simple and reliable models for individualized patient outcome monitoring and
prediction, and learning capability.

Our study has limitations, too. The study was conducted at a single
institution, which may limit the generalizability of our findings. Future
research should aim to validate the models’ performance in larger and more
diverse datasets across multiple institutions.

In conclusion, our study highlights the potential of deep learning models,
especially transformer-based models, in predicting postoperative complications
and mortality, which could revolutionize the healthcare industry’s approach to
postoperative care by improving patient outcomes, reducing healthcare costs,
and providing healthcare providers with tools to make informed clinical
decisions.

##  4 Methods

###  4.1 Data Extraction

The Department of Anesthesiology and Perioperative Medicine at the University
of Pittsburgh (Pitt) and University of Pittsburgh Medical Center (UPMC) is
widely considered to be the largest academic anesthesiology program in the
United States. UPMC has more than 300 surgical sites that are equipped with
Anesthesia Information Management System (AIMS) since 2012. UPMC uses Cerner
Powerchart solution as its EMR system and Surginet as its AIMS solution. The
necessary raw clinical data were retrieved from the UPMC clinical database by
the UPMC Health Record Research Request team based on the research team’s data
retrieval strategy. We are interested in the major surgeries that may involve
hemodynamic instability and require invasive hemodynamic monitoring. Totally
perioperative clinical data from 67,295 perioperative cases with155 different
types of surgeries are collected for further data analysis. Given the array of
metrics needed and the complexity of the data structure in the Cerner
database, a 2-stage data warehouse is constructed to reduce the need to join
and optimize multiple tables (e.g., variables and files).The first stage,
termed ”Base Tables”, is designed to serve as a middle layer decreasing the
number of tables and simplifying the joins between them. Conceptually, the
tables in Cerner’s Millennium coalesce into 3 groups: (1) patient-centered
variables (laboratories, allergies, medical history, etc.); (2) encounter
centered files (admission notes, discharge notes, pre-op anesthesia evaluation
note, etc.); and (3) patient-encounter hybrid tables (Anesthesia Record,
etc.). Tables are joined by the following fields: (1) a patient identifier
(e.g., name, medical number, fin number) in Powerchart and SurgiNet; and (2) a
case identifier (e.g., location, name of surgery, case number) in Cerner
SurgiNet.

###  4.2 Target Definition

To identify the presence of a specific postoperative complication, the study
relied on the documentation of a corresponding complication code in the EMR.
The ”present_on_admit” field in the EMR was used to indicate whether the
complication was present on admission, and if it was marked as ”N” or left
empty, it suggested that the complication developed after operation.

To identify in-hospital mortality, the study used a binary event [0, 1]. The
event was considered present if a related note, such as a death summary, death
pronouncement note, or death certificate, was found in the patient’s record
associated with the death. Importantly, the definition of in-hospital
mortality was not dependent on the patient’s length of stay in the hospital.

###  4.3 Extracted Input Features

For each hospital admission, a surgical record was created that contained 151
distinct features, which were computed or extracted at the conclusion of the
surgery. The features consisted of various intra-operative vital signs, such
as minimum and maximum blood pressure values, and a summary of drugs and
fluids interventions, including the total amount of blood infused and
Vasopressin administered. Patient anesthesia descriptions were also included,
such as the presence of an arterial line and the type of anesthesia utilized.
To review the complete list of features, refer to Table

###  4.4 Transformers

Transformers are a type of neural network architecture used for a wide range
of natural language processing (NLP) tasks such as machine translation,
sentiment analysis, and language modeling. They were introduced in the seminal
paper vaswani2017transformer  as a replacement for recurrent neural networks
(RNNs) and have since become the state-of-the-art method for many NLP tasks.

At the core of the transformer architecture is the concept of self-attention,
which allows the model to weigh the importance of different tokens in a
sequence when making predictions. Self-attention is a type of attention
mechanism that can be computed in parallel across all tokens in a sequence.
The transformer architecture uses a multi-head self-attention mechanism that
is capable of attending to multiple positions at once.

The self-attention mechanism is composed of several layers, each of which has
three main components: the query, the key, and the value. These components are
derived from the input sequence, and are used to calculate the attention
weights.

The query, key, and value components are obtained by applying learned linear
transformations to the input sequence. Specifically, given an input sequence
$\mathbf{X}=[\mathbf{x}_{1},\mathbf{x}_{2},\dots,\mathbf{x}_{n}]$, the query,
key, and value components are calculated as follows:

$\mathbf{Q}=\mathbf{X}\mathbf{W}_{Q},\mathbf{K}=\mathbf{X}\mathbf{W}_{K},\mathbf{V}=\mathbf{X}\mathbf{W}_{V}$
(1)

where $\mathbf{W}_{Q}$, $\mathbf{W}_{K}$, and $\mathbf{W}_{V}$ are learned
weight matrices.

The attention weights are then calculated using the dot product between the
query and key vectors, scaled by the square root of the dimension of the key
vectors:

$\mathrm{Attention}(\mathbf{Q},\mathbf{K},\mathbf{V})=\mathrm{softmax}\left(\frac{\mathbf{Q}\mathbf{K}^{\top}}{\sqrt{d_{k}}}\right)\mathbf{V}$
(2)

where $d_{k}$ is the dimension of the key vectors. The softmax function is
applied row-wise to ensure that the attention weights for each position in the
sequence sum to one.

Multi-head self-attention is an extension of self-attention that allows the
model to attend to information from multiple representation subspaces. In
other words, instead of computing a single set of query, key, and value
vectors for each position in the sequence, multi-head attention computes
multiple sets of vectors, each of which attends to a different subspace of the
input representation. These multiple attention heads allow the model to
capture more nuanced relationships between different parts of the input
sequence and to learn more complex interactions between them. The multi-head
self-attention result is obtained by concatenating the results of multiple
attention heads:

$\mathrm{MultiHead}(\mathbf{Q},\mathbf{K},\mathbf{V})=\mathrm{Concat}(\mathrm{head}_{1},\mathrm{head}_{2},\dots,\mathrm{head}_{k})\mathbf{W}_{O}$
(3)

$\mathrm{head}_{i}=\mathrm{Attention}(\mathbf{Q}\mathbf{W}_{i}^{Q},\mathbf{K}\mathbf{W}_{i}^{K},\mathbf{V}\mathbf{W}_{i}^{V})$
(4)

and $W_{i}^{Q}$, $W_{i}^{K}$, $W_{i}^{V}$, and $W_{O}$ are trainable weight
matrices. The output of the multi-head self-attention layer is then passed
through a feed-forward network with a non-linear activation function.

Our architecture is based on the Transformer design proposed in
gorishniy2021rdtl . The model takes the input $x$ and maps it to an embedding
space $h$ using an affine transformation given by:

$h_{j}=x_{j}.W_{j}+b_{j}$ (5)

Here, $x_{j}$ is the $j$-th feature of the data point $x$, and $W_{j}$ and
$b_{j}\in R^{d}$ are the corresponding weight and bias vectors. The final
embedding $h$ is obtained by concatenating all the embeddings of the input
data:

$h=Concat[h_{1},\cdots,h_{d}],.$ (6)

The embeddings $h$ are then processed by $L$ Transformer blocks
$f_{\theta}^{k}$ to generate the final representation $z$, which is given by:

$z=f_{\theta}^{L}(\cdots f_{\theta}^{1}(h)),.$ (7)

The final prediction $\hat{y}$ is obtained by passing a classification head
$g$, as follows:

$\hat{y}=g(z)$ (8)

Figure 4 depicts our transformer architecture.

Figure 4: Overview of the transformer-based framework. First, 1-D features are
mapped to an embedding space by an affine transformation. Then several
transformer blocks are applied to produce the final output.

###  4.5 Baselines

We conduct a comparative analysis of our framework and previous research on
tabular data prediction. Our comparison is primarily focused on variants of
Gradient Boosting Decision Trees (GBDT)friedman2001gbdt , such as XGBoost
chen2016xgboost , LightGBM ke2017lightgbm , and CatBoost
prokhorenkova2018catboost , which are widely used. However, we also examine
recently proposed deep learning methods, including Neural Oblivious Decision
Ensembles (NODE) popov2019node , GrowNet badirli2020grownet , Self-Normalizing
Networks (SNN) klambauer2017snn , as well as simple Multi-Layer Perceptron
(MLP).

###  4.6 Experimental Setup

During the training process, we used the stratified 5-fold cross-validation
method, which is an important technique for ensuring that the model is not
biased towards any particular class. In stratified cross-validation, the
dataset is split into k folds, where each fold contains roughly the same
proportion of each class as the original dataset. This helps to ensure that
the model is trained on a representative sample of the data and can generalize
well to unseen data. The use of this technique helps to prevent overfitting
and ensure that the model can accurately classify all classes, regardless of
their frequency in the dataset.

We used the Binary Cross Entropy Loss function as the loss function during
training. To ensure that the model learns to classify all classes accurately,
we used the class-balanced version of the BCE function. This assigns higher
weights to underrepresented classes and lower weights to overrepresented
classes, effectively balancing the contribution of each class to the overall
loss. The use of this function helps to prevent the model from ignoring the
minority classes and ensures that all classes are learned equally. We used the
AdamW optimizer with learning rate (lr) of 0.0003, betas of (0.9,0.999), and
weight decay of 0.001.

We used a batch size of 256 and trained the model for 200 epochs. A large
batch size helps to reduce the noise in the gradient estimation and can lead
to faster convergence. We also used early stopping with a patience of 10,
which means that the training was stopped if the validation loss did not
improve for 10 consecutive epochs. This helped to prevent overfitting and
ensure that the model was not trained for too long, which can lead to poor
generalization performance. The value of all hyperparameters were chosen based
on previous experiments and empirical observations

We trained the baselines using a set of carefully chosen hyperparameters
selected from a pool of possible values. Our aim was to ensure that the models
could achieve the highest possible performance on the given tasks. The
training settings employed to train the baselines can be found in the publicly
available GitHub repository of our paper. This repository comprises all the
essential code, scripts, and instructions required for replicating our
experiments and the results reported in the paper. Our goal in providing this
information is to facilitate other researchers’ ability to expand on our work
and explore new avenues of research.

## References

  *   * (1) Pearse, R. M. _et al._ Identification and characterisation of the high-risk surgical population in the united kingdom.  _Critical care_ 10 (3), 1–6 (2006) .
  * (2) Pearse, R. M. _et al._ Mortality after surgery in europe: a 7 day cohort study.  _The Lancet_ 380 (9847), 1059–1065 (2012) .
  * (3) Weiser, T. G. _et al._ An estimation of the global volume of surgery: a modelling strategy based on available data.  _The Lancet_ 372 (9633), 139–144 (2008) .
  * (4) Dripps, R.  New classification of physical status.  _Anesthesiology_ 24, 111 (1963) .
  * (5) Gawande, A. A., Kwaan, M. R., Regenbogen, S. E., Lipsitz, S. A. & Zinner, M. J.  An apgar score for surgery.  _Journal of the American College of Surgeons_ 204 (2), 201–208 (2007) .
  * (6) Le Manach, Y. _et al._ Preoperative score to predict postoperative mortality (pospom) derivation and validation.  _Anesthesiology_ 124 (3), 570–579 (2016) .
  * (7) Sessler, D. I. _et al._ Broadly applicable risk stratification system for predicting duration of hospitalization and mortality.  _The Journal of the American Society of Anesthesiologists_ 113 (5), 1026–1037 (2010) .
  * (8) Dalton, J. E. _et al._ Development and validation of a risk quantification index for 30-day postoperative mortality and morbidity in noncardiac surgical patients.  _The Journal of the American Society of Anesthesiologists_ 114 (6), 1336–1344 (2011) .
  * (9) NVIDIA.  Deep learning and ai solutions for industries (2018).  URL https://www.nvidia.com/en-us/deep-learning-ai/industries/.
  * (10) Lasko, T. A., Denny, J. C. & Levy, M. A.  Computational phenotype discovery using unsupervised feature learning over noisy, sparse, and irregular clinical data.  _PloS one_ 8 (6), e66341 (2013) .
  * (11) Che, Z., Kale, D., Li, W., Bahadori, M. T. & Liu, Y.  Deep computational phenotyping.  _Proceedings of the 21th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ 507–516 (2015) .
  * (12) Cheng, Y., Wang, F., Zhang, P. & Hu, J.  Risk prediction with electronic health records: A deep learning approach.  _Proceedings of the 2016 SIAM international conference on data mining_ 432–440 (2016) .
  * (13) Miotto, R., Li, L., Kidd, B. A. & Dudley, J. T.  Deep patient: an unsupervised representation to predict the future of patients from the electronic health records.  _Scientific reports_ 6 (1), 1–10 (2016) .
  * (14) Esteva, A. _et al._ Dermatologist-level classification of skin cancer with deep neural networks.  _nature_ 542 (7639), 115–118 (2017) .
  * (15) Muneeb, T., Sahu, S. & Anand, A.  Evaluating distributed word representations for capturing semantics of biomedical concepts (2015).
  * (16) Fritz, B. A. _et al._ Deep-learning model for predicting 30-day postoperative mortality.  _British journal of anaesthesia_ 123 (5), 688–695 (2019) .
  * (17) Meyer, A. _et al._ Machine learning for real-time prediction of complications in critical care: a retrospective study.  _The Lancet Respiratory Medicine_ 6 (12), 905–914 (2018) .
  * (18) Xue, Q. _et al._ Developing machine learning algorithms to predict pulmonary complications after emergency gastrointestinal surgery.  _Frontiers in Medicine_ 8, 655686 (2021) .
  * (19) Lee, C. K., Hofer, I., Gabel, E., Baldi, P. & Cannesson, M.  Development and validation of a deep neural network model for prediction of postoperative in-hospital mortality.  _Anesthesiology_ 129 (4), 649–662 (2018) .
  * (20) Chen, T. & Guestrin, C.  Xgboost: A scalable tree boosting system.  _Proc. of the 22nd acm SIGKDD_ 785–794 (2016) .
  * (21) Ke, G. _et al._ Lightgbm: A highly efficient gradient boosting decision tree.  _Advances in NeurIPS_ 30 (2017) .
  * (22) Prokhorenkova, L., Gusev, G., Vorobev, A., Dorogush, A. V. & Gulin, A.  Catboost: unbiased boosting with categorical features.  _Advances in NeurIPS_ 31 (2018) .
  * (23) Arik, S. Ö. & Pfister, T.  Tabnet: Attentive interpretable tabular learning.  _Proc. of the AAAI Conf. on AI_ 35 (8), 6679–6687 (2021) .
  * (24) Popov, S., Morozov, S. & Babenko, A.  Neural oblivious decision ensembles for deep learning on tabular data.  _arXiv preprint arXiv:1909.06312_ (2019) .
  * (25) Badirli, S. _et al._ Gradient boosting neural networks: Grownet.  _arXiv preprint arXiv:2002.07971_ (2020) .
  * (26) Klambauer, G., Unterthiner, T., Mayr, A. & Hochreiter, S.  Self-normalizing neural networks.  _Advances in NeurIPS_ 30 (2017) .
  * (27) Vaswani, A. _et al._ Attention is all you need.  _Advances in neural information processing systems_ 30 (2017) .
  * (28) Gorishniy, Y., Rubachev, I., Khrulkov, V. & Babenko, A.  Revisiting deep learning models for tabular data.  _Advances in NeurIPS_ 34, 18932–18943 (2021) .
  * (29) Friedman, J. H.  Greedy function approximation: a gradient boosting machine.  _Annals of statistics_ 1189–1232 (2001) .