0.8 in0.8 in1in1.4in

††thanks: E-mail<EMAIL_ADDRESS>E-mail:
<EMAIL_ADDRESS>
# Flying-Qubit Control via a Three-level Atom with Tunable Waveguide Couplings

Wenlong Li  Center for Intelligent and Networked Systems, Department of
Automation, Tsinghua University, Beijing 100084, China  Xue Dong  Center for
Intelligent and Networked Systems, Department of Automation, Tsinghua
University, Beijing 100084, China  Guofeng Zhang  Department of Applied
Mathematics, The Hong Kong Polytechnic University, Hong Kong, China  Re-Bing
Wu  Center for Intelligent and Networked Systems, Department of Automation,
Tsinghua University, Beijing 100084, China

###### Abstract

The control of flying qubits is at the core of quantum networks. As often
carried by single-photon fields, the flying-qubit control involves not only
their logical states but also their shapes. In this paper, we explore a
variety of flying-qubit control problems using a three-level atom with time-
varying tunable couplings to two input-output channels. It is shown that one
can tune the couplings of a $\Lambda$-type atom to distribute a single photon
into the two channels with arbitrary shapes, or use a $V$-type atom to catch
an arbitrary-shape distributed single photon. The $\Lambda$-type atom can also
be designed to transfer a flying qubit from one channel to the other, with
both the central frequency and the photon shape being converted. With a
$\Xi$-type atom, one can use the tunable coupling to shape a pair of
correlated photons via cascaded emission. In all cases, analytical formulas
are derived for the coupling functions to fulfil these control tasks, and
their physical limitations are discussed as well. These results provide useful
control protocols for high-fidelity quantum information transmission over
complex quantum networks.

##  I Introduction

The engineering of flying qubits is fundamentally important in coherent
information transmission over quantum networks Divincenzo (1995); Cirac et al.
(1997); Kimble (2008); Kuhn et al. (2002); Duan et al. (2001); Gheri et al.
(1998). Usually, the logical states of flying qubits are encoded by the number
or polarization of photons contained in a pulsed electromagnetic field.
Moreover, since the photon field often involves a continuous band of modes,
its spectrum (or equivalently its temporal shape) is also to be controlled so
as to physically match the receiver quantum nodes.

The control of flying qubits must be actuated by some manipulatable standing
quantum system (i.e., an atom). From the viewpoint of quantum input-output
theory Gardiner and Collett (1985); Combes et al. (2017), the incoming and
outgoing flying qubits form the quantum input and output of the standing
quantum system, leading to the following three classes of control tasks: (1)
the generation of flying qubits with vacuum quantum inputs Kurtsiefer et al.
(2002); Houck et al. (2007); Pierre et al. (2014); Gough and Zhang (2015);
Trivedi et al. (2018); Fischer et al. (2018), (2) the catching of flying
qubits with vacuum quantum outputs Stobinska et al. (2009); Wang et al.
(2011); Yin et al. (2013); Aljunid et al. (2013); Pierre et al. (2014); Nurdin
et al. (2016); Sounas (2020), and (3) the transformation of flying qubits with
quantum inputs and outputs both non-vacuum Zhang (2014, 2017); Leong et al.
(2016); Trivedi et al. (2018); Hurst and Kok (2018); Kiilerich and Mølmer
(2019).

Depending on the purposes of flying-qubit control, the actuating standing
quantum system may have various structures. For the generation or catching of
arbitrary-shape single flying qubits, it is sufficient to apply a simple two-
level atom Nurdin et al. (2016); Li et al. (2020, 2021). However, when
processing multiple flying qubits or converting flying qubits between
different channels, multi-level atoms are required Pan and Zhang (2017);
Schöll et al. (2020); Sbresny et al. (2022). As the simplest extension, a
three-level atom can be simultaneously coupled to two channels, which enables
the conversion of flying qubits from one channel to the other or the
generation of an entangled pair of flying qubits.

In superconducting quantum circuits, all kinds of three-level artificial atoms
can be realized with tunable frequency and coupling strength Liu et al.
(2005); Peng et al. (2016); Gu et al. (2017), including the $\Lambda-$type
Manucharyan et al. (2009), the $\Xi-$type Abdumalikov et al. (2010), the
$V-$type Srinivasan et al. (2011) and the cyclic $\Delta-$type Liu et al.
(2005) that does not exist in natural atoms. However, little attention has
been paid to the flying-qubit control with such systems. To our knowledge,
$\Lambda$-type atoms were proposed for the generation of persistent single
photons You et al. (2007); Astafiev et al. (2007), the nonreciprocal routing
of single photons among two coupled chiral waveguides Gonzalez-Ballestero et
al. (2016) and the few-photon scattering processes Hurst and Kok (2018), while
$\Xi$-type atoms can be taken as transistors for switching single photons
scattered by the atom Kyriienko and Sørensen (2016). One can also apply the
$\Delta$-type atom to the generation of microwave photon pairs Marquardt
(2007); Peng et al. (2015) and the routing of single photons Zhou et al.
(2013).

There was also not much attention paid to the shape control of continuous-mode
flying qubits Forn-Diaz et al. (2017). Since the underlying control processes
often involve time-varying parameters, the broadly used frequency-domain
scattering analysis are not applicable as it can only deal with static
systems. In the literature, various approaches have been proposed to model
such dynamical processes, including the master equations Scully and Zubairy
(1997), quantum Langevin equations Gardiner and Collett (1985), quantum
trajectories Carmichael (1993); Gough et al. (2014), quantum scattering theory
Zhang (2014); Shi et al. (2015); Zhang (2017); Trivedi et al. (2018); Fischer
et al. (2018), or quantum stochastic differential equations (QSDEs) Li et al.
(2020).

In this paper, we will adopt the QSDE-based approach to explore a variety of
flying-qubit control problems with a three-level atom Liu et al. (2014) whose
couplings to two quantum input-output channels are tunable. Analytic
conditions will be derived for the coupling functions, as well as the time-
varying detuning functions, to generate, catch or convert flying qubits with
different types of three-level atoms, and their physical limitations will be
discussed as well.

The remainder of this paper is organized as follows. Section II describes the
flying-qubit control model based on our previous work Li et al. (2020, 2021).
In Sec. III, we design control protocols for the generation of entangled
flying qubits and cascaded generation of flying qubit pairs, following which
the catching and conversion of flying qubits are discussed in Sec. IV and Sec.
V, respectively. Finally, conclusion is drawn in Sec. VI.

##  II The Model of flying-qubit control systems

In this section, we will introduce the representation of flying-qubit states
and their joint state with the standing three-level quantum system. Then, we
will provide the formulation for calculating outgoing flying-qubit states
under time-dependent controls.

Fig 1: Schematics of a standing quantum system that is coupled to two quantum
input-output channels, where the channels may physically correspond to (a) two
semi-infinite waveguides, (b) one bidirectional waveguide, or (c) two
unidirectional chiral waveguides.

###  II.1 The representation of flying-qubit states

Flying qubits are usually carried by traveling photon fields, and their
logical states are encoded by the number of photons contained in the field,
i.e., the logical states $\lvert 0\rangle$ and $\lvert 1\rangle$ are
represented by the field’s vacuum state $\lvert\rm vac\rangle$ and single-
photon state $\lvert 1_{\xi}\rangle$. Here, the shape function $\xi(t)$ is an
intrinsic property of the flying qubit in addition to the superposition of
logical states, because only properly shaped photon fields can be perfectly
received by remote quantum systems Wang et al. (2011); Aljunid et al. (2013);
Nurdin et al. (2016).

To define the single-photon state $\lvert 1_{\xi}\rangle$ with shape function
$\xi(t)$, let us start with the temporal annihilation operator
$b(t)=\int_{-\infty}^{+\infty}e^{-i\omega t}b(\omega){\rm d}\omega$ of the
traveling field in a waveguide, where $b(\omega)$ is the annihilation operator
associated with mode $\omega$. The state of a flying qubit that contains
exactly one photon in a single channel can be described as Zhang (2017, 2021)

$|1_{\xi}\rangle=\int_{-\infty}^{+\infty}\xi(\tau)b^{\dagger}(\tau)|{\rm
vac}\rangle{\rm d}\tau,$ (1)

where the shape function $\xi(t)$ represents the probability density amplitude
of finding a photon at moment $t\in(-\infty,\infty)$, and the preservation of
the overall probability requires that $\xi(t)$ be normalized, i.e.,
$\int_{-\infty}^{+\infty}|\xi(\tau)|^{2}{\rm d}\tau=1$.

The single-photon field can also be distributed in multiple channels, which
can be represented by the following superposition of $m$ single-photon
components:

$|1_{\xi}\rangle=\sum_{j=1}^{m}\int_{-\infty}^{+\infty}\xi_{j}(\tau)b^{\dagger}_{j}(\tau)|{\rm
vac}\rangle{\rm d}\tau,$ (2)

where $b_{j}^{\dagger}(\tau)$ and $\xi_{j}(\tau)$ are the field operator and
photon shape function in the $j$th channel, respectively.

More generally, we can define $n$-photon states in $m$ channels as follows

$|n_{\xi}\rangle=\sum_{j_{1},\cdots,j_{n}=1}^{m}\int_{-\infty}^{+\infty}{\rm
d}\tau_{n}\int_{-\infty}^{\tau_{n}}{\rm
d}\tau_{n-1}\cdots\int_{-\infty}^{\tau_{2}}{\rm d}\tau_{1}\
\xi_{j_{1},\cdots,j_{n}}(\tau_{1},\cdots,\tau_{n})b^{\dagger}_{j_{n}}(\tau_{n})\cdots
b^{\dagger}_{j_{1}}(\tau_{1})|{\rm vac}\rangle,$ (3)

where the shape function $\xi_{j_{1},\cdots,j_{n}}(\tau_{1},\cdots,\tau_{n})$
indicates the probability density amplitude of generating the $k$th photon in
the $j_{k}$th channel at time $\tau_{k}$. For notational consistence, we
denote the shape function with zero photon as $\xi$ that does not have any
superscripts or subscripts.

###  II.2 The computation of outgoing flying-qubit states

Consider a three-level atom that is coupled to two quantum input-output
channels. Let $|g\rangle$, $|e\rangle$ and $|f\rangle$ be the atom’s
eigenstates, where $|f\rangle$ is the highest excited state, and $\mathcal{H}$
is the Hilbert space they span. As is shown in Fig. 1, the two channels can be
physically realized by two semi-infinite waveguides, a bidirectional waveguide
or two unidirectional waveguides. Throughout this paper, we take the third
scenario for the description of problem setup, but the obtained results can be
directly extended to the other two scenarios.

To analyze the underlying quantum input-output processes, we generalize the
above flying-qubit state representation to the joint state of the interacting
atom-waveguide system, as follows:

$\begin{split}|\Psi(t)\rangle&=\sum_{n=0}^{\infty}\quad\sum_{j_{1},\cdots,j_{n}=1}^{2}\int_{-\infty}^{t}{\rm
d}\tau_{n}\int_{-\infty}^{\tau_{n}}{\rm
d}\tau_{n-1}\cdots\int_{-\infty}^{\tau_{2}}{\rm d}\tau_{1}\
|\psi_{j_{1},\cdots,j_{n}}(\tau_{1},\cdots,\tau_{n}|t)\rangle\otimes
b^{\dagger}_{j_{n}}(\tau_{n})\cdots b^{\dagger}_{j_{1}}(\tau_{1})|{\rm
vac}\rangle,\end{split}$ (4)

where
$|\psi_{j_{1},\cdots,j_{n}}(\tau_{1},\cdots,\tau_{n}|t)\rangle\in\mathcal{H}$
is the atom’s associated state at time $t$ when $n$ photons are observed at
moments $\tau_{1}\leq\cdots\leq\tau_{n}$ in the $j_{1}$th, $\cdots$, $j_{n}$th
channels, respectively. Note that each
$|\psi_{j_{1},\cdots,j_{n}}(\tau_{1},\cdots,\tau_{n}|t)\rangle$ is generally
unnormalized, but the total state $|\Psi(t)\rangle$ is always normalized.

Assume that the quantum input field (i.e., incoming flying qubits) is in the
vacuum state $|{\rm vac}\rangle$, and the atom is initially prepared at state
$|\psi(-\infty)\rangle=|\psi_{0}\rangle$. In Refs. Li et al. (2020, 2021), it
is proven that the time evolution of the joint state $\lvert\Psi(t)\rangle$ is
determined by a non-unitary evolution operator $V(t)$ on $\mathcal{H}$, which
is governed by the following equation Gardiner and Zoller (2004)

$\dot{V}(t)=\left[-iH(t)-\frac{1}{2}\sum_{j=1}^{2}L^{\dagger}_{j}(t)L_{j}(t)\right]V(t),$
(5)

where $V(-\infty)=\mathbb{I}$, $H(t)$ is the atom’s Hamiltonian including its
internal part and its interaction part with external coherent control fields,
and $L_{i}(t)$ ($i=1,2$) is its coupling operator to the $j$th channel. Our
previous studies showed that each
$|\psi_{j_{1},\cdots,j_{n}}(\tau_{1},\cdots,\tau_{n}|t)\rangle$ can be
calculated as follows:

$\begin{split}&\quad|\psi_{j_{1},\cdots,j_{n}}(\tau_{1},\cdots,\tau_{n}|t)\rangle\\\
&=V(t)\widetilde{L}_{j_{n}}(\tau_{n})\cdots\widetilde{L}_{j_{1}}(\tau_{1})|\psi_{0})\rangle,\end{split}$
(6)

where
$\widetilde{L}_{j_{i}}(\tau_{i})=V^{-1}(\tau_{i})L_{j_{i}}(\tau_{i})V(\tau_{i})$.

The outgoing flying-qubit state can be extracted from the asymptotic limit of
$|\Psi(t)\rangle$ when $t\rightarrow\infty$. Owing to the decaying nature of
$V(t)$, each component
$\quad|\psi_{j_{1},\cdots,j_{n}}(\tau_{1},\cdots,\tau_{n}|t)\rangle$ must
decay to either $|g\rangle$ or $|e\rangle$, and can thus be decomposed as:

$\displaystyle|\psi_{j_{1},\cdots,j_{n}}(\tau_{1},\cdots,\tau_{n}|\infty)\rangle$
$\displaystyle=$
$\displaystyle~{}\xi^{g}_{j_{1},\cdots,j_{n}}(\tau_{1},\cdots,\tau_{n})|g\rangle+\xi^{e}_{j_{1},\cdots,j_{n}}(\tau_{1},\cdots,\tau_{n})|e\rangle.$
(7)

Taking the asymptotic limit, we can reexpress the joint state as

$\displaystyle|\Psi(\infty)\rangle$ $\displaystyle=$
$\displaystyle|g\rangle\otimes\left[\sum_{n=0}^{\infty}\sum_{j_{1},\cdots,j_{n}}\int_{-\infty}^{t}{\rm
d}\tau_{n}\cdots\int_{-\infty}^{\tau_{2}}{\rm
d}\tau_{1}\xi^{g}_{j_{1},\cdots,j_{n}}(\tau_{1},\cdots,\tau_{n})b^{\dagger}_{j_{n}}(\tau_{n})\cdots
b^{\dagger}_{j_{1}}(\tau_{1})|{\rm vac}\rangle\right]$ (8)
$\displaystyle+|e\rangle\otimes\left[\sum_{n=0}^{\infty}\sum_{j_{1},\cdots,j_{n}}\int_{-\infty}^{t}{\rm
d}\tau_{n}\cdots\int_{-\infty}^{\tau_{2}}{\rm
d}\tau_{1}\xi^{e}_{j_{1},\cdots,j_{n}}(\tau_{1},\cdots,\tau_{n})b^{\dagger}_{j_{n}}(\tau_{n})\cdots
b^{\dagger}_{j_{1}}(\tau_{1})|{\rm vac}\rangle\right],$

in which the summations in the brackets represent the state of the emitted
photon field when the atom decays $|g\rangle$ or $|e\rangle$. When the number
of excitation numbers is preserved, the summations usually contain a finite
number of terms and hence can be effectively calculated to analyze the control
processes.

When the quantum input of some channel is non-empty (i.e., the field is
initially at a non-vacuum state), we can generalize the above procedure by
cascading an ancillary system that generates the quantum input. As is shown in
Fig. 2, the cascaded system receives vacuum quantum inputs, to which the above
procedure can be directly applied.

Suppose that the non-vacuum quantum input is fed into the first channel. Let
$H_{A}(t)$ and $L_{A}(t)$ be the Hamiltonian and the coupling operator of the
ancillary system $A$, which can be artificially chosen according to the
desired input state $|1_{\xi}\rangle$. According to the $(S,L,H)$ formula
Gough and James (2009), the equivalent Hamiltonian of the joint system is

$\displaystyle\bar{H}(t)$ $\displaystyle=$ $\displaystyle
H_{A}(t)\otimes\mathbb{I}+\mathbb{I}_{A}\otimes H(t)$ (9)
$\displaystyle+\frac{1}{2i}\left[L_{1}^{\dagger}(t)L_{A}(t)-L_{A}^{\dagger}(t)L_{1}(t)\right].$

and the equivalent coupling operators are

$\displaystyle\bar{L}_{1}(t)=$ $\displaystyle~{}\mathbb{I}_{A}\otimes
L_{1}(t)+L_{A}(t)\otimes\mathbb{I},$ (10a) $\displaystyle\bar{L}_{2}(t)=$
$\displaystyle~{}\mathbb{I}_{A}\otimes L_{2}(t).$ (10b)

The equivalent $\bar{H}(t)$ and $\bar{L}_{j}(t)$ can then be applied to
construct the nonunitary evolution equation (5) for the calculation of output
fields. The same operations can be done when the second channel or both
channels have non-vacuum inputs.

Fig 2: Schematics for transforming a non-vacuum input system (the box on the
left) to an equivalent system with vacuum input (the box on the right). The
single-photon input state $\lvert 1_{\xi}\rangle$ is generated by a properly
designed ancillary system $A$.

###  II.3 Existing results on flying-qubit control with two-level atoms

Before studying flying-qubit control actuated by a three-level atom, we
briefly review the known results with a two-level atom.

Suppose that the internal Hamiltonian of the two-level system is
$H(t)=\epsilon(t)\sigma_{z}$, where $\epsilon(t)$ is the detuning between the
atomic transition frequency and the central frequency of the field in the
channel. The coupling operator is $L(t)=\sqrt{\gamma(t)}\sigma_{-}$, where
$\sigma_{-}$ is the standard Pauli lowering operator. In previous works Li et
al. (2020, 2021), it is derived that, to generate a single-photon pulse
$\xi(t)=|\xi(t)|e^{-i\phi(t)}$ from an excited two-level atom, the coupling
function should be set to

$\bar{\bar{\gamma}}(t)=\frac{|\xi(t)|^{2}}{\int_{t}^{\infty}|\xi(s)|^{2}{\rm
d}s},$ (11)

and the detuning frequency needs to match the phase condition
$\epsilon(t)=\dot{\phi}(t)$, where $\phi(t)$ must start from
$\phi(-\infty)=0$. Hereafter, the double bars over $\gamma(t)$ are used to
denote coupling functions associated with two-level atoms.

One can also derive that, to catch a single photon $|1_{\xi}\rangle$ with
$\xi(t)=|\xi(t)|e^{-i\phi(t)}$, the coupling function $\gamma(t)$ should be
set to

$\bar{\bar{\gamma}}(t)=\frac{|\xi(t)|^{2}}{\int_{-\infty}^{t}|\xi(s)|^{2}{\rm
d}s},$ (12)

while the same phase condition needs to be matched.

The above conditions for controlling flying qubits exhibit an elegant symmetry
in that the integral in the denominator points to the future when generating a
single photon and to the past when catching a single photon. Owing the
preservation of energy, it can be verified that these two denominators both
equals the population of the excited state at the present time $t$. Later we
will see that this rule also holds for control protocols using a three-level
atom.

##  III The Generation of Flying Qubits with a Three-level Atom

In this section, we will use the above flying-qubit control model to design
control functions for generating flying qubits with a $\Lambda$-type or a
$\Xi$-type atom.

###  III.1 Entangled flying qubits generated by a $\Lambda$-type atom

Fig 3: Schematics for the generation of a distributed single photon with a
$\Lambda$-type three-level atom, which forms a pair of flying qubits that is
entangled with the atom at state $\alpha_{1}|g\rangle|1_{\xi_{1}}\rangle|{\rm
vac}\rangle+\alpha_{2}|e\rangle|{\rm vac}\rangle|1_{\xi_{2}}\rangle$. Here,
$\epsilon_{1}(t)$ and $\epsilon_{2}(t)$ are the detunings between the central
frequency of the incident field and the corresponding atomic transition
frequencies, respectively, and $\gamma_{1}(t)$ and $\gamma_{2}(t)$ are the
coupling functions.

Consider a $\Lambda$-type atom that is coupled to two quantum channels (see
Fig. 3). Being initially excited to the state $\lvert f\rangle$, the atom will
emit a single photon distributed into the two channels, forming a pair of
flying qubits that are entangled with the atom after they leave. Their joint
state can be written as

$|\Psi(\infty)\rangle=\alpha_{1}|g\rangle\otimes|1_{\xi_{1}}\rangle|{\rm
vac}\rangle+\alpha_{2}|e\rangle\otimes|{\rm vac}\rangle|1_{\xi_{2}}\rangle,$
(13)

where the coefficients $\alpha_{1}$ and $\alpha_{2}$ are complex numbers and
the shape functions $\xi_{1}(t)$ and $\xi_{2}(t)$ are normalized.

To generate such distributed single photons, we can apply the following
coupling operators:

$L_{1}(t)=\sqrt{\gamma_{1}(t)}|g\rangle\langle f|,\quad
L_{2}(t)=\sqrt{\gamma_{2}(t)}|e\rangle\langle f|,$ (14)

where the tunable coupling functions $\gamma_{1}(t)$ and $\gamma_{2}(t)$ alter
the instantaneous rate of field emission into the corresponding channels. In
the rotating-wave frame, the system’s Hamiltonian reads

$H(t)=\left[\epsilon_{1}(t)+\epsilon_{2}(t)\right]|f\rangle\langle f|,$ (15)

where $\epsilon_{1}(t)$ and $\epsilon_{2}(t)$ are the detunings between the
transition frequencies and central frequencies of the channels.

Fig 4: The generation of a pair of flying qubits at the joint state as
$\frac{1}{\sqrt{2}}(|g\rangle|1_{\xi_{1}}\rangle|{\rm
vac}\rangle+|e\rangle|{\rm vac}\rangle|1_{\xi_{2}}\rangle)$ using a
$\Lambda$-type atom. The output single photons have (a) exponentially decaying
shapes or (b) Gaussian shapes, whose corresponding coupling strengths are
shown in (c) and (d), respectively.

According to Eqs. (6)-(8), we can derive that (see Appendix A.1 for details),
for arbitrary given coefficients $\alpha_{1}$ and $\alpha_{2}$ and shape
functions $\xi_{1}(t)$ and $\xi_{2}(t)$, one can generate the entangled state
(13) using the following coupling functions

$\displaystyle\gamma_{1}(t)=\frac{|\alpha_{1}\xi_{1}(t)|^{2}}{\int_{t}^{\infty}|\alpha_{1}\xi_{1}(s)|^{2}{\rm
d}s+\int_{t}^{\infty}|\alpha_{2}\xi_{2}(s)|^{2}{\rm d}s},$ (16a)
$\displaystyle\gamma_{2}(t)=\frac{|\alpha_{2}\xi_{2}(t)|^{2}}{\int_{t}^{\infty}|\alpha_{1}\xi_{1}(s)|^{2}{\rm
d}s+\int_{t}^{\infty}|\alpha_{2}\xi_{2}(s)|^{2}{\rm d}s},$ (16b)

which have a common denominator because the coupling operators of the two
channels share the excited state $|f\rangle$.

Let $\xi_{j}(t)=|\xi_{j}(t)|e^{-i\phi_{j}(t)}$ ($j=1,2$) with $\phi_{j}(t)$
being the phase of the shape function. To match the phases, the detuning
frequencies should satisfy

$\epsilon_{1}(t)+\epsilon_{2}(t)=\dot{\phi}_{1}(t)=\dot{\phi}_{2}(t),$ (17)

where $\phi_{1}(-\infty)=\phi_{2}(-\infty)=0$. This actually requires that
$\xi_{1}(t)$ and $\xi_{2}(t)$ must have identical phase functions
$\phi_{1}(t)=\phi_{2}(t)$.

The obtained solution (16) is closely related to the single-photon generation
condition (11) with a two-level atom. They are similar in that both
denominators represent the total number of photons to be released in the
future. Because the number of excitation numbers (or energy) is preserved, the
denominators are also equal to the remaining population of the initially
excited state ($|e\rangle$ for the two-level atom and $|f\rangle$ for the
three-level atom) at present time $t$.

For illustration, we simulate two typical classes of single-photon shape
functions. In the first case, we assume that both $\xi_{1}(t)$ and
$\xi_{2}(t)$ are in the following exponentially decaying form

$|\xi_{k}(t)|=\sqrt{\gamma_{ck}}e^{-\gamma_{ck}t/2},\quad k=1,2,$ (18)

where $\gamma_{c1}/2\pi=15$MHz and $\gamma_{c2}/2\pi=5$MHz. According to Eq.
(16), it is easy to calculate that, to generate a pair of arbitrary-shape
flying qubits with the above waveforms at the joint state as
$\frac{1}{\sqrt{2}}(|g\rangle|1_{\xi_{1}}\rangle|{\rm
vac}\rangle+|e\rangle|{\rm vac}\rangle|1_{\xi_{2}}\rangle)$ , the coupling
functions need to be

$\displaystyle\gamma_{1}(t)=\frac{\gamma_{c1}}{1+e^{(\gamma_{c1}-\gamma_{c2k})t}},$
(19a)
$\displaystyle\gamma_{2}(t)=\frac{\gamma_{c2}}{1+e^{(\gamma_{c2}-\gamma_{c1})t}}.$
(19b)

Recall that, according to Eq. (11), the coupling function for generating the
two exponentially decaying waveforms with two-level atoms are
${\bar{\bar{\gamma}}}_{1}(t)\equiv\gamma_{c1}$ and
${\bar{\bar{\gamma}}}_{2}(t)\equiv\gamma_{c2}$, respectively. There is a
competition between the photon emission in the two channels in that
$\gamma_{2}(t)$ asymptotically approaches to ${\bar{\bar{\gamma}}}_{2}(t)$ but
the more rapidly decaying $\gamma_{1}(t)$ gradually vanishes.

This competition also exists when the two shapes are both Gaussian functions:

$|\xi_{k}(t)|=\left(\frac{\Omega_{k}^{2}}{2\pi}\right)^{\frac{1}{4}}e^{-(\frac{\Omega_{k}t}{2})^{2}},\quad
k=1,2,$ (20)

where $\Omega_{1}=2$MHz and $\Omega_{2}=4$MHz are the widths of the shapes. We
can see in Fig 4(d) that, as predicted, the required coupling function
$\gamma_{1}(t)$ approaches to ${\bar{\bar{\gamma}}}_{1}(t)$ because
$\xi_{1}(t)$ decays more slowly, and $\gamma_{2}(t)$ gradually vanishes.

###  III.2 Cascaded generation of correlated flying qubits by a $\Xi$-type
atom

Fig 5: Schematics for the generation of a pair of correlated single photons by
a $\Xi$-type atom, including $|1_{\xi_{1}}\rangle$ in the 1st channel and
$|1_{\xi_{2}}\rangle$ in the 2nd channel. Here, $\epsilon_{1}(t)$ and
$\epsilon_{2}(t)$ are the detuning frequencies between the incident field and
the $|f\rangle\leftrightarrow|e\rangle$ and
$|e\rangle\leftrightarrow|g\rangle$ transistion frequencies, respectively.
$\gamma_{1}(t)$ and $\gamma_{2}(t)$ are the coupling coefficients between the
$\Xi$-type atom and the channels.

Consider a $\Xi$-type atom shown in Fig. 5. When the atom is initially
prepared at the excited state $|f\rangle$, a pair of single photons will be
sequentially emitted into the first and the second channels. Since the two
photons are statistically correlated, their state cannot be decomposed as the
direct product of two single-photon states like
$|1_{\xi_{1}}\rangle\otimes|1_{\xi_{2}}\rangle$. Nevertheless, one can define
the single-photon shape observed in each channel as the corresponding marginal
probability distribution function:

$\displaystyle|\xi_{1}({\tau_{1}})|^{2}$ $\displaystyle=$
$\displaystyle\int_{\tau_{1}}^{\infty}|\xi_{1,2}^{g}(\tau_{1},\tau_{2})|^{2}{\rm
d}\tau_{2},$ (21) $\displaystyle|\xi_{2}({\tau_{2}})|^{2}$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{\tau_{2}}|\xi_{1,2}^{g}(\tau_{1},\tau_{2})|^{2}{\rm
d}\tau_{1},$ (22)

where $\xi_{1,2}^{g}(t_{1},t_{2})$ is the shape of generated two-photon field
that can be calculated from Eq. (6). In this case, $\xi_{1}(t)$ and
$\xi_{2}(t)$ do not possess definite phases $\phi_{1}(t)$ and $\phi_{2}(t)$ as
the two photons are not separable.

We seek proper coupling functions $\gamma_{1}(t)$ and $\gamma_{2}(t)$ to shape
$|\xi_{1}(t)|^{2}$ and $|\xi_{2}(t)|^{2}$ as demanded. The underlying internal
Hamiltonian and coupling operators of a $\Xi$-type atom are as follows:

$\displaystyle H(t)$ $\displaystyle=$
$\displaystyle\epsilon_{1}(t)|f\rangle\langle
f|+\epsilon_{2}(t)|e\rangle\langle e|,$ (23) $\displaystyle L_{1}(t)$
$\displaystyle=$ $\displaystyle\sqrt{\gamma_{1}(t)}|e\rangle\langle f|,$ (24)
$\displaystyle L_{2}(t)$ $\displaystyle=$
$\displaystyle\sqrt{\gamma_{2}(t)}|g\rangle\langle e|,$ (25)

from which we can derive that (see Appendix A.2 for details)

$\displaystyle\gamma_{1}(t)=\frac{|\xi_{1}(t)|^{2}}{\int_{t}^{\infty}|\xi_{1}(s)|^{2}{\rm
d}s},$ (26a)
$\displaystyle\gamma_{2}(t)=\frac{|\xi_{2}(t)|^{2}}{\int_{t}^{\infty}|\xi_{2}(s)|^{2}{\rm
d}s-\int_{t}^{\infty}|\xi_{1}(s)|^{2}{\rm d}s}.$ (26b)

The solution (26) indicates that $\gamma_{1}(t)$ is only dependent on the
shape function $|\xi_{1}(t)|^{2}$ and thus its design can be treated as that
with a two-level atom. However, the design of $\gamma_{2}(t)$ is not only
dependent on the shape function $|\xi_{2}(t)|^{2}$ but also on
$|\xi_{1}(t)|^{2}$ of the previously emitted photon.

The observed population rule in the above discussion can be verified here as
the denominator of Eq. (26a) is equal to the current population of
$|f\rangle$, whose decay leads to the emission of the first photon. Also,
owing to the conservation of energy, the denominator of Eq. (26b) represents
the accumulated population of the state $|e\rangle$ that is associated with
the emission of the second photon, as it is equal to the difference between
the energies dumped from the upper level $|f\rangle$ and to the lower level
$g\rangle$ indicated by the denominator.

Due to the passivity of the control system, the energy dumped to $|g\rangle$
can only come from the energy dumped from $|f\rangle$, and hence the following
physical realizability condition must be held:

$\int_{t}^{\infty}|\xi_{2}(s)|^{2}{\rm
d}s>\int_{t}^{\infty}|\xi_{1}(s)|^{2}{\rm d}s,\quad\forall t>0,$ (27)

which is also demanded by positivity of $\gamma_{2}(t)$ according to Eq.
(26b). This condition implies that the tail area of $|\xi_{2}(t)|^{2}$ must be
always greater than that of $|\xi_{1}(t)|^{2}$, and hence $|\xi_{1}(t)|^{2}$
should decay more quickly.

We illustrate the above result with two representative examples of Gaussian
photon pulses. In the first case, the two Gaussian pulses have different
widths and identical peak times. As shown in Fig. 6(a) and (c), the tail areas
of $|\xi_{2}(t)|^{2}$ is initially greater than that of $|\xi_{1}(t)|^{2}$,
but it becomes smaller later. This leads to the singularity of $\gamma_{2}(t)$
at the peak time and non-physical negativity after that, which indicates that
such photon shapes are not physically realizable by tunable couplings.

In the second example, the two photons have identical Gaussian shapes, but the
second photon is behind the first with a time delay, as follows:

$\begin{split}&|\xi_{1}(t)|=\left(\frac{\Omega^{2}}{2\pi}\right)^{\frac{1}{4}}e^{-\left[\frac{\Omega
t}{2}\right]^{2}},\\\
&|\xi_{2}(t)|=\left(\frac{\Omega^{2}}{2\pi}\right)^{\frac{1}{4}}e^{-\left[\frac{\Omega(t-t_{p})}{2}\right]^{2}},\end{split}$
(28)

in which $t_{p}=0.2{\rm\mu s}$. As are shown in Figs. 6(b) and (d),
$\gamma_{2}(t)$ is physically realizable because the condition (27) holds.
Because the tail area of $|\xi_{1}(t)|^{2}$ decays more quickly,
$\gamma_{2}(t)$ asymptotically approaches to ${\bar{\bar{\gamma}}}_{2}(t)$
that is associated to a two-level aotom.

Fig 6: The generation of a pair of correlated single photons
$(|1_{\xi_{1}}\rangle|1_{\xi_{2}}\rangle$ using a $\Xi$-type atom. The
incident photon shapes $\xi_{1}(t)$ and $\xi_{2}(t)$ are with (a) different
Gaussian shapes or (b) same Gaussian shapes. The corresponding coupling
strengths are shown in (c) and (d) respectively.

##  IV The catching of flying qubits by a three-level atom

Fig 7: Schematics for catching a flying qubit at state $|1_{\xi}\rangle|{\rm
vac}\rangle$ by a $\Lambda$-type atom at state $|g\rangle$. Here,
$\epsilon_{1}(t)$ and $\epsilon_{2}(t)$ are the detuning frequencies between
the incident field and the $|f\rangle\leftrightarrow|g\rangle$ and
$|f\rangle\leftrightarrow|e\rangle$ transistion frequencies, respectively.
$\gamma_{1}(t)$ and $\gamma_{2}(t)$ are the coupling coefficients between the
$\Lambda$-type atom and the channels.

Let us first consider the $\Lambda$-type atom shown in Fig. 7. When the atom
is initially prepared at the state $|g\rangle$, and we wish to catch a single
photon $|1_{\xi}\rangle$ from the first channel, i.e., finding proper control
functions such that the ouputs of both channels are empty, and the atom’s
state transits from $|g\rangle$ to $|f\rangle$. Assuming that
$\xi(t)=|\xi(t)|e^{-i\phi(t)}$, we can introduce an ancillary system $A$ to
generate $|1_{\xi}\rangle$, and apply Eqs. (6)-(8) to the joint standing
system. It can be derived that the coupling function and detuning of the first
channel must be

$\gamma_{1}(t)=\frac{|\xi(t)|^{2}}{\int_{-\infty}^{t}|\xi(s)|^{2}{\rm
d}s},\quad\epsilon_{1}(t)=\dot{\phi}(t).$ (29)

Meanwhile, the coupling to the second channel must be turned off, i.e.,
$\gamma_{2}(t)\equiv 0$, so as to prevent the leakage of the photon into the
second channel (The details of its derivation can see in Appendix A.3). Thus,
the design of such controls is exactly the same as that with a two-level atom.

Fig 8: Schematics for catching flying qubits at state
$\alpha_{1}|1_{\xi_{1}}\rangle|{\rm vac}\rangle+\alpha_{2}|{\rm
vac}\rangle|1_{\xi_{2}}\rangle$ by a $V$-type atom at state $|g\rangle$. Here,
$\epsilon_{1}(t)$ and $\epsilon_{2}(t)$ are the detuning frequencies between
the incident field and the $|f\rangle\leftrightarrow|g\rangle$ and
$|f\rangle\leftrightarrow|e\rangle$ transistion frequencies, respectively.
$\gamma_{1}(t)$ and $\gamma_{2}(t)$ are the coupling functions between the
atom and the channels.

We can also use a $V$-type atom to catch a distributed single photon. As is
shown in Fig. 8, the photon fields coming from both channels, forming an
entangled state $\alpha_{1}|1_{\xi_{1}}\rangle|{\rm
vac}\rangle+\alpha_{2}|{\rm vac}\rangle|1_{\xi_{2}}\rangle$, where
$\xi_{k}(t)=|\xi_{k}(t)|e^{-i\phi_{k}(t)}$ ($k=1,2$). When the atom is
initially prepared at the ground state $|g\rangle$, we can prove that the
coupling functions and detuning frequencies for the entangled qubits to be
caught by the three-level atom is

$\gamma_{k}(t)=\frac{|\xi_{k}(t)|^{2}}{\int_{-\infty}^{t}|\xi_{k}(s)|^{2}{\rm
d}s},\quad\epsilon_{k}(t)=\dot{\phi}_{k}(t),$ (30)

where $k=1,2$ (the derivation is straightforward and will not be provided
here). Consequently, the entangled state is transferred to the atom’s
superposition state $\alpha_{1}|f\rangle+\alpha_{2}|e\rangle$.

The obtained solutions can still be interpreted by the rule we observed in the
above analysis, that is, the denominators are equal to the population of the
associated excited state accumulated from the absorption of the corresponding
single photon.

##  V The conversion of flying qubits by a three-level atom

Fig 9: Schematics for the conversion of flying qubits from state
$|1_{\xi_{1}}\rangle|{\rm vac}\rangle$ to $|{\rm
vac}\rangle|1_{\xi_{2}}\rangle$ by a $\Lambda$-type atom. Here,
$\epsilon_{1}(t)$ and $\epsilon_{2}(t)$ are the detuning frequencies between
the incident field and the $|f\rangle\leftrightarrow|g\rangle$ and
$|f\rangle\leftrightarrow|e\rangle$ transition frequencies, respectively.
$\gamma_{1}(t)$ and $\gamma_{2}(t)$ are the coupling coefficients between the
$\Lambda$-type atom and the channels.

It is natural to use a $\Lambda$-type atom to convert a single photon between
the two channels that it couples to. The conversion can not only reshape the
single photon but also change its carrier frequency (note that $\xi(t)$ is the
envelope of the single-photon shape function in the rotating-wave frame).

As is shown in Fig. 9, we want to design control protocols such that the atom
converts a single-photon (at state $|1_{\xi_{1}}\rangle$) in the first channel
to a single photon (at state $|1_{\xi_{2}}\rangle$) in the second channel.
This consists of the catching of the first photon via transition between
$|g\rangle$ and $|f\rangle$ and then the release of the second photon via
transition between $|e\rangle$ and $|f\rangle$.

Similarly, we can introduce an ancillary system to generate
$|1_{\xi_{1}}\rangle$ and obtain the conditions for the coupling functions and
detuning frequencies to satisfy. In Appendix A.3, we derive that

$\begin{split}\gamma_{1}(t)&=\frac{|\xi_{1}(t)|^{2}}{\int_{-\infty}^{t}|\xi_{1}(\tau)|^{2}{\rm
d}\tau-\int_{-\infty}^{t}|\xi_{2}(\tau)|^{2}{\rm d}\tau},\\\
\gamma_{2}(t)&=\frac{|\xi_{2}(t)|^{2}}{\int_{t}^{\infty}|\xi_{2}(\tau)|^{2}{\rm
d}\tau-\int_{t}^{\infty}|\xi_{1}(\tau)|^{2}{\rm d}\tau},\end{split}$ (31)

and the detuning frequencies are

$\epsilon(t)=\dot{\phi}_{1}(t)=\dot{\phi}_{2}(t).$ (32)

Moreoever, for physically realizable conversions, the phase functions must
satisfy

$\phi_{1}(t)-\phi_{2}(t)=(2n+1)\pi,\quad n\in Z.$

This is different from the phase condition for the generation of single
photons.

The obtained solution for flying-qubit conversion is also physical realizable
under the condition (27) on the tail areas of $|\xi_{1}(t)|^{2}$ and
$|\xi_{2}(t)|^{2}$. This is because the outgoing photon must be released after
the incoming photon is caught, or equivalently saying, the radiated energy
into the outgoing photon must be less than the energy absorbed from the
incoming photon. Their difference, as described by the common denominator in
Eq. (31), is nothing but the gained population of the state $|f\rangle$ shared
by the two channels. This again conforms to the population rule we observed in
the above three control scenarios.

Hence, the shape of the incoming photon $\xi_{1}(t)$ must decay more quickly
than that of the outgoing photon, and consequently, $\gamma_{1}(t)$ will
asymptotically approach to ${\bar{\bar{\gamma}}}_{1}(t)$ (the solution for
catching $|1_{\xi_{2}}\rangle$ with a two-level atom) when
$t\rightarrow-\infty$, and $\gamma_{2}(t)$ will approach to
${\bar{\bar{\gamma}}}_{2}(t)$ (the solution for generating
$|1_{\xi_{2}}\rangle$ with a two-level atom) when $t\rightarrow\infty$.

Fig 10: The conversion of a flying qubit $|1_{\xi_{1}}\rangle$ in the first
channel into a flying qubit $|1_{\xi_{2}}\rangle$ using a $\Lambda$-type atom.
The incident photon shapes $\xi_{1}(t)$ and $\xi_{2}(t)$ are with (a)
different Gaussian shapes or (b) same Gaussian shapes. The corresponding
coupling strengths are shown in (c) and (d) respectively.

For verification, we test the same group of Gaussian shapes simulated in Sec.
III. As is shown in Figs. 10(a) and (c), the violation of the tail area
condition leads to non-physical negative coupling functions $\gamma_{1}(t)$
and $\gamma_{2}(t)$ after the peak time. In Figs 10(b) and (d), the resulting
coupling functions are physically realizable, and one can clearly observe the
predicted asymptotic limits of $\gamma_{1}(t)$ and $\gamma_{2}(t)$ when
$t\rightarrow\pm\infty$.

##  VI Conclusion and outlook

To conclude, we have explored the control of flying qubits using various types
of three-level atoms through tunable couplings to two input-output channels,
including the generation of entangled qubits and correlated flying qubits, as
well as the catching and conversion of single flying qubits. Analytic formulas
are presented for the coupling and detuning functions to fulfil these control
tasks, which lay a basis for systematic design of control protocols for high-
fidelity flying-qubit transmission over complex quantum networks.

The method applied here can be naturally generalized to more complicated cases
that include more atomic levels and input-output channels, as long as the
number of excitation number (or energy) is preserved. We conjecture that,
owing to the conservation law, the population rule concluded in this paper
will still hold, i.e., the flying-qubit control problems can be decomposed
into several single-photon generation or catching problems, and the
corresponding coupling function is equal to the ratio between the square norm
of the desired photon shape and the current population of the excited state
associated with this photon. Whether this holds for most general cases is to
be examined in our future studies.

In addition to the three types of three-level atoms we studied here, it is
more intriguing to apply a $\Delta$-type atom, which is only realizable by
aritificial atoms, for the control of flying qubits. For example, one can use
it to convert a single photon into a pair of correlated photon pair, or vice
versa, or routing a single photon between two different channels. However, the
underlying scattering process is much more complicated, and we have not
obtained any useful analytic solutions.

It should noted that the obtained solutions may confront physical constraints
on the realizable shape functions $\xi_{1}(t)$ and $\xi_{2}(t)$, such as the
phase conditions and the restriction (27) on the pulse tail areas. Towards
more flexible flying-qubit shaping controls, one may have to introduce
coherent controls that breaks the conservation of the number of excitations.
Under such circumstances, analytic solutions are usually unavailable, but
numerical optimization methods can be introduced to assist the design. All
these problems will be explored in our future studies.

##  VII Acknowledgments

This work is supported by the National Natural Science Foundation of China
under Grant Nos. 61833010, 62173201 and 6217023269. GZ acknowledges supports
from the Hong Kong Research Grant council (RGC) grants (Nos. 15203619 and
15208418), the Shenzhen Fundamental Research Program under the Grant No.
JCYJ20190813165207290, and the CAS AMSS-PolyU Joint Laboratory of Applied
Mathematics.

##  Appendix A Proofs of main results

For conciseness., we introduce several notations as follows:

$\displaystyle\Gamma_{j}(t)$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{t}\gamma_{j}(\tau){\rm d}\tau,$ (33)
$\displaystyle\Theta_{j}(t)$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{t}\epsilon_{j}(\tau){\rm d}\tau,$ (34)

where $\gamma_{j}(\tau)$ and $\epsilon_{j}(\tau)$ ($j=1,2$) are the coupling
and detuning functions of the $j$th channel. Based on the above notations, we
further denote $\gamma(t)=\gamma_{1}(t)+\gamma_{2}(t)$ and
$\epsilon(t)=\epsilon_{1}(t)+\epsilon_{2}(t)$.
$\Gamma(t)=\Gamma_{1}(t)+\Gamma_{2}(t)$ and
$\Theta(t)=\Theta_{1}(t)+\Theta_{2}(t)$.

In the following, we will prove the main results presented in the main text.

###  A.1 The derivation of Eq. (16) and Eq. (17)

Applying Eqs. (5)-(8), it is straightforward to find that, with the
Hamiltonian and coupling operator given in Sec. III.1, only the single-photon
emissions can occur when the three-level atom is initially prepared at
$|f\rangle$ and ends up in $|g\rangle$ or $|e\rangle$. The calculation of Eq.
(II.2) on the single-photon emission gives

$\displaystyle\xi_{1}^{g}(\tau)$ $\displaystyle=$
$\displaystyle\sqrt{\gamma_{1}({\tau})}e^{-\frac{1}{2}\Gamma(s)-i\Theta(s)},$
(35) $\displaystyle\xi_{2}^{e}(\tau)$ $\displaystyle=$
$\displaystyle\sqrt{\gamma_{2}({\tau})}e^{-\frac{1}{2}\Gamma(s)-i\Theta(s)}.$
(36)

To generate the target distributed single photon that are entangled with the
atom at state $\alpha_{1}|g\rangle|1_{\xi_{1}}\rangle|{\rm
vac}\rangle+\alpha_{2}|e\rangle|{\rm vac}\rangle|1_{\xi_{2}}\rangle$, we
should have $\xi_{1}^{g}(\tau)=\alpha_{1}\xi_{1}(\tau)$ and
$\xi_{2}^{e}(\tau)=\alpha_{2}\xi_{2}(\tau)$, and hence

$\displaystyle\alpha_{1}\xi_{1}(\tau)$ $\displaystyle=$
$\displaystyle\sqrt{\gamma_{1}({\tau})}e^{-\frac{1}{2}\Gamma(s)-i\Theta(s)},$
(37) $\displaystyle\alpha_{2}\xi_{2}(\tau)$ $\displaystyle=$
$\displaystyle\sqrt{\gamma_{2}({\tau})}e^{-\frac{1}{2}\Gamma(s)-i\Theta(s)}.$
(38)

Let $\xi_{k}(\tau)=\left|\xi_{k}({\tau})\right|e^{-i\phi_{k}(t)}$, $k=1,2$.
The phase condition Eq. (17) can be immediately obtained by comparing the
phases of the left and right hand sides of Eqs. (37)-(38). These two equations
also imply

$\frac{|\alpha_{1}\xi_{1}(\tau)|^{2}}{|\alpha_{2}\xi_{2}(\tau)|^{2}}=\frac{\gamma_{1}(\tau)}{\gamma_{2}(\tau)},$
(39)

and

$\left|\alpha_{1}\xi_{1}({\tau})\right|^{2}+\left|\alpha_{2}\xi_{2}({\tau})\right|^{2}=\gamma(\tau)e^{-\Gamma(\tau)}=-\frac{\rm
d}{{\rm d}\tau}e^{-\Gamma(\tau)}.$ (40)

The integration of both sides of Eq. (40) from $-\infty$ to $t$ gives

$\int_{-\infty}^{t}\left[\left|\alpha_{1}\xi_{1}({\tau})\right|^{2}{\rm
d}\tau+\left|\alpha_{2}\xi_{2}({\tau})\right|^{2}\right]{\rm
d}\tau=-e^{-\Gamma(\tau)},$ (41)

from which we solve

$\gamma(t)=\frac{\left|\alpha_{1}\xi_{1}(t)\right|^{2}+\left|\alpha_{2}\xi_{2}(t)\right|^{2}}{\int_{-\infty}^{t}\left[\left|\alpha_{1}\xi_{1}({\tau})\right|^{2}+\left|\alpha_{2}\xi_{2}({\tau})\right|^{2}\right]{\rm
d}\tau}.$ (42)

The solution (16) of $\gamma_{1}(t)$ and $\gamma_{2}(t)$ can thus be obtained
by (39) and (42).

###  A.2 The derivation of Eq. (26)

We apply Eqs. (5)-(8) with the model of $\Xi$-type atom described in Sec.
III.2. It is found that, when the three-level atom is initially prepared at
$|f\rangle$ and ends up in $|g\rangle$, only the two-photon parts are left in
Eq. (II.2), as follows

$\xi_{1,2}^{g}(\tau_{1},\tau_{2})=\sqrt{\gamma_{1}({\tau_{1}})\gamma_{2}({\tau_{2}})}e^{-\frac{1}{2}\int_{-\infty}^{\tau_{1}}\gamma_{1}(s){\rm
d}s-\frac{1}{2}\int_{\tau_{1}}^{\tau_{2}}\gamma_{2}(s){\rm
d}s}e^{-i\int_{-\infty}^{\tau_{1}}\epsilon_{1}(s){\rm
d}s-i\int_{\tau_{1}}^{\tau_{2}}\epsilon_{2}(s){\rm d}s}.$ (43)

According to the control target, the desired shape functions
$|\xi_{1}({\tau_{1}})|^{2}$ and $|\xi_{2}({\tau_{2}})|^{2}$ should be the
marginal distributions of $|\xi_{1,2}^{g}(\tau_{1},\tau_{2})|^{2}$,
respectively, from which we derive that

$|\xi_{1}({\tau_{1}})|^{2}=\int_{\tau_{1}}^{\infty}|\xi_{1,2}^{g}(\tau_{1},\tau_{2})|^{2}{\rm
d}\tau_{2}=\gamma_{1}({\tau_{1}})e^{-\Gamma_{1}(\tau_{1})},$ (44)

and

$\begin{split}|\xi_{2}({\tau_{2}})|^{2}&=\int_{-\infty}^{\tau_{2}}|\xi_{1,2}^{g}(\tau_{1},\tau_{2})|^{2}{\rm
d}\tau_{1}\\\
&=\gamma_{2}(\tau_{2})e^{-\Gamma_{2}(\tau_{2})}\int_{-\infty}^{\tau_{2}}|\xi_{1}({\tau_{1}})|^{2}e^{\Gamma_{2}(\tau_{1})}{\rm
d}\tau_{1}.\end{split}$ (45)

Equation (44) can be simplified to obtain Eq. (26a) for the solution of
$\gamma_{1}(t)$. As for the solution of $\gamma_{2}(t)$, we can rewrite Eq.
(45) as

$\int_{-\infty}^{\tau_{2}}|\xi_{1}({\tau_{1}})|^{2}e^{\Gamma_{2}(\tau_{1})}{\rm
d}\tau_{1}=\frac{|\xi_{2}({\tau_{2}})|^{2}}{\gamma_{2}(\tau_{2})}e^{\Gamma_{2}(\tau_{2})},$
(46)

and differentiate both sides of Eq. (46) to obtain the solution (26b).

###  A.3 The derivation of Eq. (31) and Eq. (32)

We first construct an auxiliary system $A$ to virtually generate the incident
single photon $|1_{\xi_{1}}\rangle$, for which one can set

$H_{A}(t)=0,\quad L_{A}(t)=\sqrt{\gamma_{A}(t)}\sigma_{-},$ (47)

with the coupling function being

$\gamma_{A}(t)=\frac{|\xi_{1}(t)|^{2}}{\int_{t}^{\infty}|\xi_{1}(s)|^{2}{\rm
d}s},$ (48)

where, $\gamma_{A}(t)$ is the coupling function between the auxiliary system
and first channel. According to the $(S,L,H)$ formula Gough and James (2009),
the equivalent Hamiltonian of the joint system as Eq. (9) is

$\begin{split}\bar{H}(t)&=\mathbb{I}_{A}\otimes\left[\epsilon_{1}(t)+\epsilon_{2}(t)\right]|f\rangle\langle
f|\\\ &+\frac{1}{2i}\sqrt{\gamma_{A}(t)\gamma_{1}(t)}\left[|g,f\rangle\langle
e,g|-H.c\right],\end{split}$ (49)

where the notation $|\alpha,\beta\rangle=|\alpha\rangle\otimes|\beta\rangle$
is adopted. The corresponding equivalent coupling operators as Eq. (10) are

$\begin{split}\bar{L}_{1}(t)&=\mathbb{I}_{A}\otimes\sqrt{\gamma_{1}(t)}|g\rangle\langle
f|+\sqrt{\gamma_{A}(t)}|g\rangle\langle e|\otimes\mathbb{I},\\\
\bar{L}_{2}(t)&=\mathbb{I}_{A}\otimes\sqrt{\gamma_{2}(t)}|e\rangle\langle
f|.\end{split}$ (50)

Applying Eq. (6) with the above equivalent $\bar{H}(t)$ and $\bar{L}_{j}(t)$,
we can see that only the following three components of correlated system
states are non-vanishing

$\begin{split}&|\psi(t)\rangle=e^{-\frac{1}{2}\Gamma_{A}(t)}|e,g\rangle-\Xi(t)e^{-\frac{1}{2}\Gamma(t)}|g,f\rangle,\\\
&|\psi_{1}(\tau|t)\rangle=\left[\xi_{1}(\tau)-\sqrt{\gamma_{1}({\tau})}e^{-\frac{1}{2}\Gamma(\tau)-i\Theta(\tau)}\Xi({\tau})\right]|g,g\rangle,\\\
&|\psi_{2}(\tau|t)\rangle=-\sqrt{\gamma_{2}({\tau})}e^{-\frac{1}{2}\Gamma(\tau)-i\Theta(\tau)}\Xi({\tau})|g,e\rangle,\end{split}$
(51)

where $\Gamma_{A}(t)=\int_{-\infty}^{t}\gamma_{A}(s){\rm d}s$ and

$\Xi(t)=\int_{-\infty}^{t}\sqrt{\gamma_{1}(\tau)}e^{\frac{1}{2}\Gamma(\tau)+i\Theta(\tau)}\xi_{1}(\tau){\rm
d}\tau.$ (52)

The three non-vanishing components correspond to the cases when both channels
have vacuum output, the first channel has a single-photon output, and the
second channel has a single-photon output, respectively.

Note that $|\psi_{1}(\tau|t)\rangle$ and $|\psi_{2}(\tau|t)\rangle$ are both
independent of $t$, and hence we have

$\displaystyle\xi_{1}^{g}(\tau)\\!$ $\displaystyle=$
$\displaystyle\\!\xi_{1}(\tau)-\sqrt{\gamma_{1}({\tau})}e^{-\frac{1}{2}\Gamma(\tau)-i\Theta(\tau)}\Xi({\tau}),$
(53) $\displaystyle\xi_{2}^{e}(\tau)\\!$ $\displaystyle=$
$\displaystyle\\!-\sqrt{\gamma_{2}({\tau})}e^{-\frac{1}{2}\Gamma(\tau)-i\Theta(\tau)}\Xi({\tau}),$
(54)

To convert the incident flying qubit $|1_{\xi_{1}}\rangle$ into the outgoing
flying qubit $|1_{\xi_{2}}\rangle$, the single photon must completely go out
through the second channel, and hence $\xi_{1}^{g}(\tau)=0$ and
$\xi_{1}^{e}(\tau)=\xi_{2}(\tau)$, which implies that the two shape functions
satisfy

$\displaystyle\xi_{1}(\tau)$ $\displaystyle=$
$\displaystyle\sqrt{\gamma_{1}({\tau})}e^{-\frac{1}{2}\Gamma(\tau)-i\Theta(\tau)}\Xi({\tau}),$
(55) $\displaystyle\xi_{2}(\tau)$ $\displaystyle=$
$\displaystyle-\sqrt{\gamma_{2}({\tau})}e^{-\frac{1}{2}\Gamma(\tau)-i\Theta(\tau)}\Xi({\tau}).$
(56)

Now we apply the normalization condition of the atom-field joint state (4),
which requires that

$\begin{split}1&\equiv\langle\Psi(t)|\Psi(t)\rangle=\langle\psi(t)|\psi(t)\rangle\\\
&\quad+\int_{-\infty}^{t}\left[\langle\psi_{1}(\tau|t)|\psi_{1}(\tau|t)\rangle+\langle\psi_{2}(\tau|t)|\psi_{2}(\tau|t)\rangle\right]{\rm
d}\tau.\end{split}$ (57)

Combining with Eqs. (51), (55) and (56), we can transform Eq. (57) as follows:

$1\equiv\int_{t}^{\infty}|\xi_{1}(\tau)|^{2}{\rm
d}\tau+|\Xi(t)|^{2}e^{-\Gamma(t)}+\int_{-\infty}^{t}|\xi_{2}(\tau)|^{2}{\rm
d}\tau,$ (58)

where the relation $e^{-\Gamma_{A}(t)}=\int_{t}^{\infty}|\xi_{1}(s)|^{2}{\rm
d}s$ is used.

This equation precisely describes the conservation of excitation number (or
energy in the system). The system is initially prepared at the ground state
with a single-photon input, and hence its total number of excitations is 1, as
indicated by the left hand side of Eq. (58). At any time $t>0$, the energy
contained in the photon may stay in the input field (the first term on the
right hand side) or be transferred to the population of $|f\rangle$ (the
second term on the right hand side ) and the output field (the third term on
the right hand side). When $t\rightarrow\infty$, all the energy flows form the
incoming photon to the outgoing photon, while the first and second terms decay
to zero.

Further, Eq. (58) implies that

$\begin{split}|\Xi(t)|^{2}e^{-\Gamma(t)}&=\int^{t}_{-\infty}|\xi_{1}(\tau)|^{2}{\rm
d}\tau-\int_{-\infty}^{t}|\xi_{2}(\tau)|^{2}{\rm d}\tau,\\\
&=\int_{t}^{\infty}|\xi_{2}(\tau)|^{2}{\rm
d}\tau-\int^{\infty}_{t}|\xi_{1}(\tau)|^{2}{\rm d}\tau,\end{split}$ (59)

as both $\xi_{1}(\tau)$ and $\xi_{2}(\tau)$ are normalized functions.
Replacing this equation back into Eqs. (55) and (56), we obtain the solutions
for the coupling functions:

$\begin{split}\gamma_{1}(t)&=\frac{|\xi_{1}(t)|^{2}}{\int_{-\infty}^{t}|\xi_{1}(\tau)|^{2}{\rm
d}\tau-\int_{-\infty}^{t}|\xi_{2}(\tau)|^{2}{\rm d}\tau},\\\
\gamma_{2}(t)&=\frac{|\xi_{2}(t)|^{2}}{\int_{t}^{\infty}|\xi_{2}(\tau)|^{2}{\rm
d}\tau-\int_{t}^{\infty}|\xi_{1}(\tau)|^{2}{\rm d}\tau}.\end{split}$ (60)

As for the solutions of the detuning functions, one can directly observe that
Eqs. (55) and (56) are satisfied when

$\epsilon(t)=\dot{\phi}_{1}(t),\quad\epsilon(t)=\dot{\phi}_{2}(t),$ (61)

where $\phi_{1}(t)$ and $\phi_{2}(t)$ are the phase functions of $\xi_{1}(t)$
and $\xi_{2}(t)$, respectively. Moreover, Eqs. (55) and (56) indicates that
$\xi_{1}(t)$ and $\xi_{2}(t)$ must have inverted phases, i.e.,

$\phi_{1}(t)-\phi_{2}(t)=(2n+1)\pi,\quad n\in Z.$

The obtained equations also provide the solution for catching an incident
single photon from the first channel, which requires that
$\xi_{1}^{g}(\tau)=\xi_{2}^{e}(\tau)=0$. According to Eq. (56),
$\gamma_{2}(t)$ must be turned off so that the resulting design is equivalent
to that with a two-level atom.

## References

  * Divincenzo (1995) D. Divincenzo, Science 270, 255 (1995).
  * Cirac et al. (1997) J. I. Cirac, P. Zoller, H. J. Kimble, and H. Mabuchi, Physical Review Letters 78, 3221 (1997).
  * Kimble (2008) H. J. Kimble, Nature 453, 1023 (2008).
  * Kuhn et al. (2002) A. Kuhn, M. Hennrich, and G. Rempe, Physical Review Letters 89, 067901 (2002).
  * Duan et al. (2001) L. Duan, M. D. Lukin, J. I. Cirac, and P. Zoller, Nature 414, 413 (2001), ISSN 1476-4687.
  * Gheri et al. (1998) K. Gheri, K. Ellinger, T. Pellizzari, and P. Zoller, Fortschritte der Physik 46, 401 (1998).
  * Gardiner and Collett (1985) C. W. Gardiner and M. J. Collett, Physical Review A 31, 3761 (1985).
  * Combes et al. (2017) J. Combes, J. Kerckhoff, and M. Sarovar, Advances in Physics: X 2, 784 (2017).
  * Kurtsiefer et al. (2002) C. Kurtsiefer, S. Mayer, P. Zarda, and H. Weinfurter, Physical Review Letters 85, 290 (2002).
  * Houck et al. (2007) A. Houck, D. Schuster, J. Gambetta, J. Schreier, B. Johnson, J. Chow, L. Frunzio, J. Majer, M. Devoret, S. Girvin, et al., Nature 449, 328 (2007).
  * Pierre et al. (2014) M. Pierre, I. Svensson, S. Sathyamoorthy, G. Johansson, and P. Delsing, Applied Physics Letters 104, 232604 (2014).
  * Gough and Zhang (2015) J. Gough and G. Zhang, EPJ Quantum Technology 2 (2015).
  * Trivedi et al. (2018) R. Trivedi, K. Fischer, S. Xu, S. Fan, and J. Vuckovic, Physical Review B 98, 144112 (2018).
  * Fischer et al. (2018) K. Fischer, R. Trivedi, V. Ramasesh, I. Siddiqi, and J. Vučković, Quantum 2, 69 (2018).
  * Stobinska et al. (2009) M. Stobinska, G. Alber, and G. Leuchs, Europhysics Letters 86, 14007 (2009).
  * Wang et al. (2011) Y. Wang, J. Minar, L. Sheridan, and V. Scarani, Physical Review A 83, 063842 (2011).
  * Yin et al. (2013) Y. Yin, Y. Chen, D. Sank, P. O’Malley, T. C. White, R. Barends, J. Kelly, E. Lucero, M. Mariantoni, A. Megrant, et al., Physical Review Letters 110, 107001 (2013).
  * Aljunid et al. (2013) S. Aljunid, G. Maslennikov, Y. Wang, H. Dao, V. Scarani, and C. Kurtsiefer, Physical Review Letters 111, 103001 (2013).
  * Nurdin et al. (2016) H. I. Nurdin, M. R. James, and N. Yamamoto, in _55th IEEE Conference on Decision and Control (CDC)_ (2016), pp. 2513–2518.
  * Sounas (2020) D. Sounas, Physical Review B 101, 104303 (2020).
  * Zhang (2014) G. Zhang, Automatica 50, 442 (2014).
  * Zhang (2017) G. Zhang, Automatica 83, 186 (2017).
  * Leong et al. (2016) V. Leong, M. Seidler, M. Steiner, A. Cere, and C. Kurtsiefer, Nature Communications 7, 13716 (2016).
  * Hurst and Kok (2018) D. Hurst and P. Kok, Physical Review A 97, 043850 (2018).
  * Kiilerich and Mølmer (2019) A. Kiilerich and K. Mølmer, Physical Review Letters 123, 123604 (2019).
  * Li et al. (2020) W. Li, G. Zhang, and R. Wu, in _20th IFAC World Congress_ (2020), vol. 50, pp. 11755 – 11759, 20th IFAC World Congress.
  * Li et al. (2021) W. Li, G. Zhang, and R. Wu, _On the control of flying qubits_ (2021).
  * Pan and Zhang (2017) Y. Pan and G. Zhang, Journal of Physics A: Mathematical and Theoretical 50, 345301 (2017).
  * Schöll et al. (2020) E. Schöll, L. Schweickert, L. Hanschke, K. D. Zeuner, F. Sbresny, T. Lettner, R. Trivedi, M. Reindl, S. F. C. da Silva, R. Trotta, et al., Physical Review Letters 125 (2020).
  * Sbresny et al. (2022) F. Sbresny, L. Hanschke, E. Schöll, W. Rauhaus, B. Scaparra, K. Boos, E. Zubizarreta Casalengua, H. Riedl, E. del Valle, J. J. Finley, et al., Phys. Rev. Lett. 128, 093603 (2022).
  * Liu et al. (2005) Y.-x. Liu, J. Q. You, L. F. Wei, C. P. Sun, and F. Nori, Phys. Rev. Lett. 95, 087001 (2005).
  * Peng et al. (2016) Z. H. Peng, S. E. D. Graaf, J. S. Tsai, and O. V. Astafiev, Nature Communications 7 (2016).
  * Gu et al. (2017) X. Gu, A. Kockum, A. Miranowicz, Y.-x. Liu, and F. Nori, Physics Reports 718, 1 (2017).
  * Manucharyan et al. (2009) V. E. Manucharyan, J. Koch, L. I. Glazman, and M. H. Devoret, Science 326, 113 (2009).
  * Abdumalikov et al. (2010) A. A. Abdumalikov, O. Astafiev, A. M. Zagoskin, Y. A. Pashkin, Y. Nakamura, and J. S. Tsai, Phys. Rev. Lett. 104, 193601 (2010).
  * Srinivasan et al. (2011) S. J. Srinivasan, A. J. Hoffman, J. M. Gambetta, and A. A. Houck, Phys. Rev. Lett. 106, 083601 (2011).
  * You et al. (2007) J. Q. You, Y.-x. Liu, C. P. Sun, and F. Nori, Phys. Rev. B 75, 104516 (2007).
  * Astafiev et al. (2007) O. Astafiev, K. Inomata, A. O. Niskanen, T. Yamamoto, Y. A. Pashkin, Y. Nakamura, and J. S. Tsai, Nature 449, 1476 (2007).
  * Gonzalez-Ballestero et al. (2016) C. Gonzalez-Ballestero, E. Moreno, F. J. Garcia-Vidal, and A. Gonzalez-Tudela, Phys. Rev. A 94, 063817 (2016).
  * Kyriienko and Sørensen (2016) O. Kyriienko and A. S. Sørensen, Phys. Rev. Lett. 117, 140503 (2016).
  * Marquardt (2007) F. Marquardt, Phys. Rev. B 76, 205416 (2007).
  * Peng et al. (2015) Z. H. Peng, Y.-x. Liu, J. T. Peltonen, T. Yamamoto, J. S. Tsai, and O. Astafiev, Phys. Rev. Lett. 115, 223603 (2015).
  * Zhou et al. (2013) L. Zhou, L.-P. Yang, Y. Li, and C. P. Sun, Phys. Rev. Lett. 111, 103604 (2013).
  * Forn-Diaz et al. (2017) P. Forn-Diaz, C. W. Warren, C. W. S. Chang, A. M. Vadiraj, and C. M. Wilson, Physical Review Applied 8, 054015 (2017).
  * Scully and Zubairy (1997) M. Scully and M. Zubairy, _Quantum optics_ (Cambridge University Press, New York, 1997).
  * Carmichael (1993) H. Carmichael, _An Open Systems Approach to Quantum Optics_ (Springer Berlin Heidelberg, Berlin, 1993).
  * Gough et al. (2014) J. Gough, M. James, and H. Nurdin, New Journal of Physics 16, 075008 (2014).
  * Shi et al. (2015) T. Shi, D. Chang, and J. Cirac, Physical Review A 92, 053834 (2015).
  * Liu et al. (2014) Y.-X. Liu, H.-C. Sun, Z. Peng, A. Miranowicz, J. Tsai, and F. Nori, Scientific Reports 4 (2014).
  * Zhang (2021) G. Zhang, Control Theory and Technology 19, 544 (2021).
  * Gardiner and Zoller (2004) C. Gardiner and P. Zoller, _Quantum Noise_ (Springer-Verlag, Berlin, 2004).
  * Gough and James (2009) J. Gough and M. R. James, IEEE Transactions on Automatic Control 54, 2530 (2009).


|  Bouncing-to-wetting transition for impact of water droplets on soft solids
---|---
|  Surjyasish Mitra a, Quoc Vo b,∗, and Tuan Tran b,∗∗
|  Soft surfaces impacted by liquid droplets trap more air underneath than
their rigid counterpart. The extended lifetime of the air film not only
facilitates bouncing behaviours of the impacting droplets but also increases
the possibility of an interaction between the air film itself and the air
cavity formed inside the droplets by capillary waves. Such interaction may
cause rupture of the trapped air film by a so-called dimple inversion
phenomenon and suppress bouncing. In this work, we systematically investigate
the relation between air cavity collapse and air film rupture for water
droplets impacting on soft, hydrophobic surfaces. By constructing a bouncing-
to-wetting phase diagram based on the rupturing dynamics of the trapped air
film, we observe that the regime in which air film rupture is induced by
dimple inversion consistently separates the bouncing regime and the one in
which wetting is caused by random rupture. We also find that air film rupture
by dimple inversion, in-turn, affects both the collapsing dynamics of the air
cavity and the resulting high-speed jet. We then provide a detailed
characterisation of the collapsing dynamics of the air cavity and subsequent
jetting.

##

††footnotetext: a School of Physical & Mathematical Sciences, Nanyang
Technological University, 50 Nanyang Avenue, 639798 Singapore††footnotetext: b
School of Mechanical & Aerospace Engineering, Nanyang Technological
University, 50 Nanyang Avenue, 639798 Singapore. E-mail:∗<EMAIL_ADDRESS>∗∗
<EMAIL_ADDRESS>
##  1 Introduction

The air film separating an impacting liquid droplet with an impacted surface
plays a crucial role in dictating impact outcomes. For high velocity impacts,
the lubrication pressure built-up in the air film causes liquid splashes 1, 2.
For low velocity impacts, the presence of a sustained air film leads to
bouncing of the droplets regardless of surface’s wettability 3, 4, 5, 6, 7.
Studies of the bouncing-to-wetting transition of low velocity water droplets
impacting on solid surfaces are motivated by both fundamental and practical
interests. The former motivation comes from numerous intangible physical
phenomena, e.g., capillary waves 8, air entrapment dynamics 5, involved in
dictating such transition, while the latter one is from applications requiring
design of dynamic surfaces such as anti-bacterial, self-cleaning surfaces 9 or
improvement of industrial processes including droplet deposition 10 and ink-
jet printing 11. Although numerous investigations have been focusing on the
bouncing-to-wetting transition, in particular on the dynamics of the
intervening air film at the moment wetting occurs 3, 4, 5, 6, 7, 12, the
mechanisms causing such transition remains elusive.

Typically, the prelude to the final touchdown between an approaching droplet
and a solid surface is the formation of a thin air film in which lubrication
pressure is built up. The lubrication pressure subsequently becomes
sufficiently large that it deforms the droplet’s bottom surface, creating a
central dimple surrounded by an outer edge with one or two kinks, the regions
where the air film thickness is minimum 13, 6, 5. For impact on hydrophilic
surfaces, bouncing is ensured at low impact velocity. For higher impact
velocity, the air film typically ruptures either at inner or outer kink
leading to wetting initiation 6, 5. In other words, bouncing-to-wetting
transition for droplets impacting smooth hydrophilic surfaces, e.g., glass or
mica, is mainly determined by air film’s dynamics and the surface properties
that cause random wetting initiation 6, 5.

Fig. 1: (a) Schematic of the experimental setup. (b) Impact dynamics of a 1.16
mm water droplet on soft PDMS 60:1 substrates with Young’s modulus
$E=4.8\,{\rm kPa}$ and two impact velocities: 0.33 ${\rm m\,s^{-1}}$
($\rm{We}=1.7$) (upper panel), 0.39 ${\rm m\,s^{-1}}$ ($\rm{We}=2.4$) (lower
panel). The snapshots corresponding to the last column in both panels show two
different impact outcomes: bouncing and deposition with entrapped air bubble.
Scale bars represent 0.5 mm.

Impact outcomes on soft surfaces are markedly different 14, 15, 16. On the one
hand, an elastomer or gel coated surface is usually hydrophobic 17. Thus,
impacts of water droplets on such soft surfaces appear to have characteristics
similar to those on smooth, rigid hydrophobic surfaces 18, 19. An impacting
droplet above a critical impact velocity may develop a pyramidal structure due
to capillary waves induced upon impact. This leads to formation of a
cylindrical air cavity that penetrates deeply into the droplet and eventually
collapses, shooting out a liquid jet 14, 15, 18, 19. On the other hand, due to
the deformation of the soft solids, air film trapped beneath the impacting
droplets have a longer lifetime compared to those impacting on a rigid one 16.
The sustained air film not only inhibits wetting initiation and facilitates
bouncing of the impacting droplet 14, 15, but also increases the probability
for interaction between the air film and the air cavity formed in the liquid
bulk. We note that the latter condition is rarely met for impacts on smooth,
rigid hydrophilic or hydrophobic surfaces due to random air film rupture and
wetting initiation 6, 5, 18, 19. Whether the interactions between the trapped
air film and the air cavity dictates the wetting initiating mechanism of the
bouncing-to-wetting transition for impact of water droplets on soft surfaces
remains an open question.

In this work, we experimentally study the bouncing-to-wetting transition of
water droplets impacting on soft surfaces. We hypothesise that such transition
is determined by interactions between the air cavity and the air film trapped
beneath the impacting droplets. To test this hypothesis, we first use both
high-speed laser interferometry and high-speed optical imaging to characterize
the rupture mechanism of the air film and the dynamics of the air cavity in
the bulk. We reveal how interaction between the air film and the cavity’s
collapsing dynamics results in different types of air film rupture. We then
construct a phase diagram describing the dependence of the bouncing-to-wetting
transition on the stiffness of the substrate and the impact characteristics
based on the characterized air film rupture mechanisms. Finally, we provide an
explanation for the collapsing dynamics of the air cavity and the resulting
high-speed jets and bubble entrapment.

##  2 Experiments

Our test substrates were glass slides coated with a layer of
Polydimethylsiloxane (PDMS) of thickness $70\,\upmu$m. The elasticity $E$ of
the substrate was varied from $4.8\,$kPa to $263.6\,$kPa by changing the ratio
of the monomer to crosslinker from 60:1 to 30:1 (Tab. 1). The equilibrium
contact angle $\theta_{\rm Y}$ of a small water droplet on the test substrates
varied between 103.6∘ and 113.0∘ (Tab. 1). Impact experiments were conducted
using water droplets with radius $r_{0}$ varying from $0.63\,$mm to $1.4\,$mm;
the velocity $v$ of the impacting droplets was varied between $0.30$ ${\rm
m\,s^{-1}}$ and $0.51$ ${\rm m\,s^{-1}}$. The Weber number, defined as ${\rm
We}=\rho r_{0}v^{2}/\gamma$, thus varies from $2.1$ to $5.1$. Here,
$\rho=1000\,{\rm kg\,m^{-3}}$ and $\gamma=72\,{\rm Nm^{-1}}$ respectively are
the density and surface tension of water.

Table 1: Substrates used in our experiments. The numbers show the weight ratio of the monomer to crosslinker. The value of stiffness $E$ is measured experimentally using a rheometer. Equilibrium contact angles are measured using the sessile drop technique. Substrates: | P60:1 | P50:1 | P40:1 | P35:1 | P30:1
---|---|---|---|---|---
$E$ (kPa) | 4.8 | 29.5 | 70.9 | 142.7 | 263.6
$\theta_{\rm Y}$ | 109.5∘ $\pm$ 3.5∘ | 108.2∘ $\pm$ 1.7∘ | 107.7∘ $\pm$ 0.8∘ | 107.0∘ $\pm$ 1.3∘ | 106.8∘ $\pm$ 3.2∘

Impacting droplets were recorded synchronously from the bottom and the side
using two high-speed cameras (SA-X2 and SA-5, Photron) operating at imaging
rates from 30,000 to 200,000 frames-per-second and shutter time 1/800,000.
This setup has a temporal resolution at 3.75 $\mu$s. A cold LED light source
was used for side view illumination, while a red diode laser (wavelength
$\lambda=640$ nm) was used to illuminate the impacted surface from below (Fig.
1a). The laser illumination coupled with a 5$\times$ optical zoom lens enabled
us to observe the dimple profile with a height resolution of $\lambda/4=160$
nm 6, 20, and a horizontal resolution of 3.5 $\mu$m/pixel.

##  3 Results and Discussions

###  3.1 Bouncing-to-wetting transition

For impacts of low viscosity liquids, the intrinsic length scale of the
capillary waves induced upon impact is the wavelength
${\rm{\lambda}_{c}}\sim\gamma{\rho}^{-1}{v}^{-2}$ 18. The condition for
capillary wave formation is $\lambda_{\rm c}>r_{0}$, or equivalently ${\rm
We}>1$. As the Weber number in our experiment varies between 2 and 5, both
inertial and capillary forces are significant. In this Weber number range,
impacting droplets deform into pyramidal shapes, as shown in Fig. 1b
($t=3.4\,$ms).

From side-view images, we observe two distinct macroscopic behaviours when the
velocity $v$ of the impacting droplet increases: the droplet either bounces
off from the substrate at low impact velocity (Fig. 1b, upper panel) or is
deposited onto the substrate as the impact velocity in increased (Fig. 1b,
lower panel). In both cases, we observe that the air cavity forms roughly at
the moment the droplet reaches the maximum deformation and starts retracting
(Fig. 1b, $t=6.7\,$ms). Subsequently, the air cavity collapses resulting in
liquid jets from the droplet. Typically, we observe that the jet velocity is
higher in the case of depositing droplets compared to that of bouncing
droplets from the substrate.

Fig. 2:  Air film rupture dynamics of 1.4 mm radius water drops impacting on
soft PDMS surfaces at different impact velocity $v$. (a) At $v=0.36\,{\rm
m\,s^{-1}}$, the droplet bounces of after impact without wetting the solid
surface. (b) At $v=0.38\,{\rm m\,s^{-1}}$, a perfect inversion of the air
cavity causes rupture of the air film at the dimple’s center. Liquid jet and
air bubble entrapment caused by the collapse of air cavity are observed. (c)
At $v=0.47\,{\rm m\,s^{-1}}$, rupture of air film occurs at the rim of the air
cavity. (d) At $v=0.51\,{\rm m\,s^{-1}}$, rupture of air film occurs randomly
at the rim location where the air film is thinnest. The scale bars in both
side-view and bottom-view images represent 1 mm.

To reveal mechanisms causing the different impact outcomes, i.e., bouncing or
deposition, in Fig. 2, we show bottom-view interferometric recordings and the
corresponding extracted profiles of the air film trapped between the liquid
and the solid upon impact. In all cases, the liquid droplet initially skates
on a thin film of air with the bottom liquid-air interface expanding at a
lateral speed $\sim$ 0.2 ${\rm m\,s^{-1}}$, in close proximity to those
reported in a recent study for impacts in a similar range of Weber number 21.
While expanding, the bottom surface of the droplet deforms into the familiar
dimple profile 4, 6 with two distinct kinks due to excess pressure in the
trapped air film, evident from the recorded interference signatures (see Fig.
2 a-d, $t=0.34$ ms).

The evolution of the trapped air film is sensitive to the impact velocity. At
low impact velocity, $v\leq 0.36\,{\rm m\,s^{-1}}$ (Fig. 2a), the air film
remains intact during the entire duration of impact. Therefore, the droplet
bounces off from the substrate similar to previous studies involving droplet
bouncing on glass 5 and PDMS 15 surfaces. For $0.38\,{\rm m\,s^{-1}}\leq v\leq
0.47\,{\rm m\,s^{-1}}$, the air cavity in the bulk forces the dimple downward
causing shape inversion. The air film separating the inverted dimple and the
soft solid eventually ruptures, initiating wetting. We observe two distinct
types of wetting initiation caused by the inversion of the dimple: wetting
either first occurs at the center of the dimple (Fig. 2b, $t=6.76$ ms), or at
the dimple’s inner rim (Fig. 2c, $t=6.84$ ms). While wetting initiation at
dimple’s center had been previously observed in a similar study of droplet
impact on PDMS surfaces 14, wetting initiation at dimple’s inner rim due to
dimple inversion has not been reported. In our experiment, the wetting
initiation at the dimple’s inner rim consistently happens at impact velocity
slightly higher than that at the dimple’s center. For a higher range of impact
velocity, $v\geq 0.51\,{\rm m\,s^{-1}}$, rupture of the air film initiates at
a random position near the inner or outer kink (Fig. 2d, $t=0.34$ ms)
consistent with previous drop impact studies on rigid and soft surfaces 4, 6,
15.

Fig. 3:  Phase diagram showing the general behaviours, i.e., bouncing, dimple
inversion and wetting, of droplets impacting on soft solid. The behaviours are
obtained by varying Weber number ${\rm We}$ and Young modulus $E$ of the soft
substrates. The droplet is water having fixed radii $r_{0}=1.4\,$mm. Fig. 4:
Snapshots showing the formation and collapse of air cavity during impacts of
$r_{0}=1.4\,$mm water droplets onto PDMS 40:1 substrates with $E=70.9\,{\rm
kPa}$ at different impact velocities $v$. From (a) to (f), the impact velocity
$v$ increases from $0.33\,{\rm m\,s^{-1}}$ to $0.46\,{\rm m\,s^{-1}}$
resulting in different collapse dynamics of the air cavity. The dashed boxes
show the snapshots at the moment the air film ruptures. The schematics show
the air cavity shapes right before collapse, i.e., $t=t_{\rm c}$ (see the
dashed-dotted box). The arrows in each schematic show the position where the
collapse occurs. In all cases, the scale bars represent 1 mm.

In Fig. 3, we show a phase diagram presenting how the Weber number ${\rm We}$
and surface elasticity $E$ affect the behaviour of the trapped air film and
the resulting impact outcomes. Generally, we observe that dimple inversion
causing air film rupture at the center of the dimple is sandwiched between two
other major behaviours: bouncing and random rupture. Wetting initiated by
rupture at the dimple’s rim (open triangles) appears as a transitional
behaviour separating the center rupture and the random rupture behaviours. In
other words, for a fixed soft substrate, increasing the ${\rm We}$ number
causes the impact dynamics to change from bouncing to center rupture, rim
rupture and finally random rupture. For surfaces with $E>142$ kPa, we do not
observe the rim rupture behaviour. We also observe that increasing the
substrate elasticity $E$, while causing insignificant effect to the transition
between bouncing and wetting considerably reduces the transitional ${\rm We}$
number between center rupture and rim rupture. This is qualitatively
consistent with the fact that increasing the surface stiffness results in less
air entrapment and thinner air film, which eventually leads to random rupture
of the air film and wetting initiation5, 6.

###  3.2 Collapse dynamics of air cavity

The dimple and air film dynamics not only affect the general impact outcomes
but also alter the formation and collapse of the air cavity formed upon
impact. In Fig. 4, we show several series of snapshots highlighting the
formation and collapse of the air cavity upon impact of $r_{0}=1.4\,$mm water
droplets on soft PDMS surfaces having $E=70.9\,\rm{kPa}$. At low impact
velocity, i.e., $v=0.33\,{\rm m\,s^{-1}}$, we observe that a U-shaped cavity
is partially formed (Fig. 4a). Subsequently, the droplet’s surface restores to
its minimum surface area pushing the cavity out. No liquid jet or air bubble
entrapment inside the liquid bulk is observed. At $v\geq 0.36\,{\rm
m\,s^{-1}}$, the capillary waves on the droplet’s surface are sufficiently
strong to form a cylindrical cavity through the droplet width (Fig. 4b-f, the
first snapshots). Subsequently the cavity radius retracts and eventually
collapses resulting in liquid jet and bubble entrapment. Depending on the
impact velocity $v$, the cavity either collapses at the bottom of the
cylindrical cavity (Fig. 4b, $t=7.36\,{\rm ms}$), or at the cavity’s waist
(Fig. 4c, $t=6.82\,{\rm ms}$), cavity’s top (Fig. 4d, $t=7.16\,{\rm ms}$) or
at both top and bottom (Fig. 4e, $t=7.60\,{\rm ms}$). We attribute the
different collapsing dynamics of the cavity to the capillary waves on the
cavity surface which is clearly observed in (Fig. 4b-e): the collapse position
occurs at the wave’s peak. At high impact velocity, where the air film
ruptures randomly, (Fig. 4f), the wetting initiation at the rupture point
disturbs the retracting cavity in the bulk.

Fig. 5: Snapshots showing changes in temporal order of air cavity collapse and
wetting initiation by varying droplet radius: (a) $r_{0}=1.4\,$mm, wetting
occurs before cavity collapse; (b) $r_{0}=1.16\,$mm, wetting initiates right
at the moment cavity collapses; (c) $r_{0}=0.63\,$mm, wetting initiates after
cavity completely collapsed. Parameters: (a) $E=29.5\,$kPa, ${\rm We}=2.8$;
(b) $E=4.8\,$kPa, ${\rm We}=2.4$; (c) $E=4.8\,$kPa, ${\rm We}=2.3$. Scale bars
in both the side-view and bottom-view snapshots represent 1 mm.

Interestingly, we observe that the relative temporal order between wetting
initiation and complete collapse of the air cavity depends on both impact
velocity and droplet radius. As shown in Fig. 4c-e, with a fixed droplet
radius ($r_{0}=1.4\,$mm), the duration $t_{\rm c}-t_{\rm r}$ decreases from
$0.36\,$ms to $0$ ms when the impact velocity $v$ increases from $0.38\,{\rm
m\,s^{-1}}$ (Fig. 4c) to $0.45\,{\rm m\,s^{-1}}$ (Fig. 4e). Here, $t_{\rm r}$
is the moment at which rupture by dimple inversion happens (dashed-blue boxes
in Fig. 4). In Fig. 5, we show that wetting occurs _before_ the final collapse
of the air cavity for larger droplets (see Fig. 5a, $r_{0}=1.4$ mm) and
_after_ the collapse of the air cavity for smaller droplets (see Fig. 5c,
$r_{0}=0.63$ mm). The radius at which the temporal order between wetting
initiation and cavity collapse switches is $r_{0}=1.16$ mm (Fig. 5b). This is
consistent with a previously reported experimental observation by Chen et al.
14 in which wetting initiation at the dimple center occurs after the air
cavity completely collapses for water droplets of radius $r_{0}=1\,$mm and a
similar Weber number impacting on PDMS coated substrates.

Fig. 6: (Color online) Plots showing the dependence of air cavity radius
$r_{\rm c}$ with $t_{\rm c}-t$ for different impact outcomes, i.e., bounce or
wet, and different experimental parameters: (a) droplets bounce,
$r_{0}=1.4\,$mm, $v=0.36\,{\rm m\,s^{-1}}$, varying $E$; (b) droplets bounce,
$E=4.8\,$kPa, ${\rm We}\approx 2$, varying $r_{0}$; (c) droplets wet the solid
by dimple inversion, $r_{0}=1.4\,$mm, $v=0.39\,{\rm m\,s^{-1}}$, varying $E$;
(d) droplet wet the solid by dimple inversion, $E=4.8\,$kPa, varying both
$r_{0}$ and impact velocity $v$.

To quantitatively characterize the collapsing dynamics of the air cavity, we
measure the evolution of the retracting cavity radius $r_{\rm c}$. In Fig. 6,
we show the dependence of $r_{\rm c}$ on $t_{\rm c}-t$ for different impact
outcomes, i.e., bounce or wet, while varying the surface elasticity $E$,
droplet radius $r_{0}$, and impact velocity $v$. Here, $r_{\rm c}$ is the
cavity radius measured along the location where the air cavity eventually
collapses, $t_{\rm c}$ the time instant when the retracting cavity completely
collapses. It should be noted here that we only show values of $r_{\rm c}$
measured in the case that the air cavity symmetrically collapses, which case
occurs when the droplets bounces or wet the substrate by dimple inversion
(Fig. 4b-e). In the case that random rupture occurs, it is not possible to
accurately measure $r_{\rm c}$ as the cavity collapses asymmetrically (Fig.
4f). We observe two distinct collapsing behaviours depending on the impact
outcomes: $r_{\rm c}\sim(t_{\rm c}-t)^{1/2}$ for bouncing droplets (Fig.
6a,b), and $r_{\rm c}\sim(t_{\rm c}-t)^{2/3}$ for wetting induced by dimple
inversion (Fig. 6c,d). We observe that elasticity of the substrate $E$ does
not affect the dependence of $r_{\rm c}$ on $t_{\rm c}-t$ in both bouncing and
dimple inversion regimes. In other words, although surface elasticity
significantly affects the transition between bouncing and dimple inversion,
the formation and collapse of the air cavity is mainly determined by the
hydrodynamical properties of the droplets.

For impacts in the bouncing regime (Figs. 2a and 4a,b), it was shown in
previous studies 18, 19, 22 that the air cavity’s dynamics is dominated by
inertia, leading to the relation $r_{\rm c}\sim(\gamma
r_{0})^{1/4}{\rho}^{-1/4}(t_{\rm c}-t)^{1/2}$. This is indeed consistent with
the power law $r_{c}\sim A(t_{\rm c}-t)^{1/2}$ observed in our experiment
(Fig. 6a,b). Fitting this power law to our experimental data collected for
$r_{0}$ in the range from $0.63\,\rm{mm}$ to $1.4\,\rm{mm}$ respectively
yields $A$ from $0.010$ to $0.022$ consistent with the calculated values of
$(\gamma r_{0})^{1/4}{\rho}^{-1/4}$ from $0.014$ to $0.020$.

For impacts in which the air film ruptures due to dimple inversion, we observe
a different power-law behaviour for the collapsing cavity radius. We argue
that the fast contact line motion at the rupture point generates capillary
waves on the cavity’s surface (Fig. 4); these waves act as the driving factor
for the subsequent collapse of the cavity. As the collapsing dynamics is in
the inertial regime, the capillary waves are self-similar with phase velocity
$u\sim\gamma^{1/3}[\rho(t_{\rm c}-t)]^{-1/3}$ 23, 24, 25, 26 in the present
context. Consequently, balancing the dynamical pressure $\rho u^{2}$ with the
Laplace pressure $\gamma r_{\rm c}^{-1}$ at the collapsing position results in
an expression for retracting cavity radius:

$r_{\rm c}\sim{{\left(\frac{\gamma}{\rho}\right)}^{1/3}}{(t_{\rm
c}-t)^{2/3}}.$ (1)

Eq. 1 indicates that the dependence of $r_{\rm c}$ on $t_{\rm c}-t$ in the
case that dimple inversion happens does not depend on either droplet radius
$r_{0}$ or the impact velocity $v$, consistent with our experimental data
shown in Fig. 6d.

###  3.3 Liquid jet dynamics

Fig. 7: (a) Snapshots showing the jetting dynamics by air cavity collapsing at
its waist. Parameters: for the left snapshot (bounce), $E=4.8\,$kPa,
$v=0.34\,{\rm m\,s^{-1}}$, $r_{0}=1.16\,$mm; for the right snapshot (wet),
$E=70.9\,$kPa, $v=0.38\,{\rm m\,s^{-1}}$, $r_{0}=1.4\,$mm. (b) Snapshots
showing the jetting dynamics by air cavity collapsing at its ends, i.e.,
bottom end, top end or both. Parameters: for the left snapshot (bounce),
$E=70.9\,$kPa, $v=0.36\,{\rm m\,s^{-1}}$, $r_{0}=1.4\,$mm; for the right
snapshot (wet) $E=29.5\,$kPa, $v=0.45\,{\rm m\,s^{-1}}$, $r_{0}=1.4\,$mm; (c)
Snapshots showing the rupture moment by dimple inversion (upper panel) and the
subsequent high-speed jet (lower panel) of a special case in which the jet
originates from the bottom of the air cavity. This suggests that the high-
speed jet originating from the bottom of the air cavity directly relates to
dimple inversion phenomenon. At the time the jet appears in the air cavity,
i.e., at $t=7.8\,$ms, the liquid already wets all the impact area (snapshots
not shown here). Parameters: $E=4.8\,$kPa, $v=0.47\,{\rm m\,s^{-1}}$,
$r_{0}=1.4\,$mm. The scale bars in all the side-view snapshots represent 1 mm.

When the air cavity collapses due to necking, it is divided into two air
cavities which subsequently evolve differently under surface tension effects.
On one hand, the enclosed geometry of the lower air cavity prevents the air
from escaping and eventually results in an air bubble entrapped inside the
liquid bulk. On the other hand, the upper air cavity, which is connected to
the ambient air at its top continues to collapse and finally generates a
liquid jet as illustrated in Fig. 1b. As a result, the jetting dynamics, i.e.,
the jet’s radius and velocity, strongly depends on the collapsing dynamics of
the air cavity.

In Fig. 7a and b, we show two different jetting dynamics as a result of air
cavity collapse dynamics, i.e., cylindrical cavity (Fig. 4c) or conical cavity
(Fig. 4b, d). When jetting is induced by cylindrical air cavity collapse, we
observe that the jets have small radius and high velocity (Fig. 7a). The jet’s
shape does not depend on whether the droplets bounce or wet the solid.
However, we observe that the jetting velocity is much higher when the droplets
wet the solid by dimple inversion compared to that when it bounces. For
wetting by dimple inversion, we observe jet velocity as high as $37\,{\rm
m\,s^{-1}}$. Whereas, for bouncing scenarios, the jet velocity observed is
typically in the range $2\,{\rm m\,s^{-1}}-10\,{\rm m\,s^{-1}}$. When the air
cavity collapses at its ends, the so-called conical collapse (Fig. 7b), the
jets have larger radius and lower velocity compared to those of the
cylindrical collapse. Typical jet velocity for conical collapse is
approximately 2 ${\rm m\,s^{-1}}$. Whether the air film ruptures or not does
not significantly affect either the radius or velocity of the jets. The
maximum jet velocity observed in our experiments is approximately 100 times
the impact velocity. This is higher than the maximum jet velocity of droplets
impact on hydrophilic, hydrophobic, or superhydrophobic surfaces reported in
previous studies 18, 19, 27, 28, 29, i.e., $16-55$ times the impact velocity.
We note that the lower bound of jet velocity in our case, i.e., approximately
1 ${\rm m.s^{-1}}$ is similar to those observed in the mentioned studies 18,
19, 27, 28, 29. The jet velocity appears to be more dependent on the nature of
air cavity at collapse, i.e., conical or cylindrical, rather than the
elasticity of the soft surfaces or their hydrophobic character.

For the lower air cavity which fails to escape, we consistently observe that
it eventually gets trapped inside the liquid bulk as an air bubble. From our
experiments, we find the radius of the entrapped air bubble to vary between
100 $\mu$m - 400 $\mu$m. The entrapped air bubble freely floats around in the
liquid bulk.

In Fig. 7c, we show snapshots of a special case in which a thin liquid jet
emanates from the bottom of the retracting air cavity and shoots upwards
($t\geq 7.8\,$ms) after wetting initiates at the center due to dimple
inversion ($t=7.48\,$ms). This jet is similar to the one reported by Chen et
al. 19 for impacting of droplets on structured PDMS surfaces. We note here
that the onset of the jet occurs approximately $0.30-0.35\,$ms after air film
rupture and wetting initiation suggesting that the jet directly links to
dimple inversion causing air film rupture at the impact center. Typically, the
jet radius varies between 15 $\upmu$m and 45 $\upmu$m. It appears that the
liquid jet expands radially briefly, for a duration of 60 $\upmu$s - 80
$\upmu$s, before terminating into tiny droplets. The liquid jetting is
immediately followed by air cavity collapse and its accompanying jet. Hence,
for such impact scenarios, we observe jets from both thin liquid film collapse
and air cavity collapse.

##  4 Concluding remarks

In summary, we have experimentally investigated the impact characteristics of
water droplets on soft PDMS surfaces at low velocity. We observed that soft
surfaces promotes bouncing for a higher range of impact Weber number as
compared to smooth hydrophilic or hydrophobic surfaces 6, 5. Further, using
interference, we probed the interaction between the collapsing air cavity in
the bulk of an impacting droplet and the air film separating the droplet with
the impacted soft surface. Above a threshold impact velocity, the interaction
induces inversion of the dimple in the air film and subsequently ruptures the
air film. This behaviour is in a stark contrast to wetting initiation caused
by random air film rupture commonly witnessed 6, 16. We then constructed a
bouncing-to-wetting phase diagram based on the observed air film rupture
dynamics. We experimentally confirmed that for impact of water droplets on
soft surfaces, bouncing-to-wetting transition is determined by interactions
between the air cavity in the bulk and the air film trapped beneath the
droplets, the so-called dimple inversion phenomenon. The stiffness of the
substrates, on the other hand, does not significantly alter the bouncing-to-
wetting transition. We have also found that, rupture of the air film by dimple
inversion, in-turn, pins the base of the retracting bulk cavity to the surface
and influences its collapsing dynamics. In such scenarios, the collapsing
radius exhibits self-similarity and obeys a $2/3$ power-law in time, similar
to the inertia-capillary break-away of a tapered fluid sheet 30 or other self-
similar phenomena observed in free surface flows 23, 31. For such collapsing
mechanism, we further observed that the lower half of the collapsed air cavity
causes entrapment and formation of an air bubble, while the upper half induces
a liquid jet with velocity which can be as high as 100 folds of the impact
velocity.

We highlight that within the tested velocity range, dimple inversion and its
subsequent rupture causes wetting initiation at the center, as opposed to
random rupture locations 6, 16. This mode of liquid-solid contact is a more
controlled event where the liquid-solid footprint radius expands in a radially
symmetric manner, similar to onset of quasi-static spreading of liquid drops
on surfaces 32, 33. Such controlled wetting characteristics may contribute to
minimizing air bubble entrapment upon droplet impact and thus would be helpful
in applications such as droplet deposition 10 or ink-jet printing 11.

## Conflicts of interest

There are no conflicts to declare.

## Acknowledgments

The authors thank Maurice H. W. Hendrix for valuable discussions regarding the
procedure to extract dimple profile shape from interference fringes. The
authors also thank Marcus Lin for his help with rheometer measurements of the
soft surfaces. This study is supported by Nanyang Technological University
(NTU) and A*STAR (SERC Grant No. 1523700102). S.M. is supported by NTU
Research Scholarship.

## Notes and references

  * Xu _et al._ 2005 L. Xu, W. W. Zhang and S. R. Nagel, _Phys. Rev. Lett._ , 2005, 94, 184505\.
  * Driscoll and Nagel 2011 M. M. Driscoll and S. R. Nagel, _Phys. Rev. Lett._ , 2011, 107, 154502\.
  * Kolinski _et al._ 2012 J. M. Kolinski, S. M. Rubinstein, S. Mandre, M. P. Brenner, D. A. Weitz and L. Mahadevan, _Phys. Rev. Lett._ , 2012, 108, 074503.
  * de Ruiter _et al._ 2012 J. de Ruiter, J. M. Oh, D. van den Ende and F. Mugele, _Phys. Rev. Lett._ , 2012, 108, 074505.
  * de Ruiter _et al._ 2015 J. de Ruiter, R. Lagraauw, D. van den Ende and F. Mugele, _Nat. Phys._ , 2015, 11, 48.
  * de Ruiter _et al._ 2015 J. de Ruiter, D. van den Ende and F. Mugele, _Phys. Fluids_ , 2015, 27, 012105.
  * van der Veen _et al._ 2012 R. C. van der Veen, T. Tran, D. Lohse and C. Sun, _Phys. Rev. E_ , 2012, 85, 026315.
  * Yarin 2006 A. L. Yarin, _Annu. Rev. Fluid Mech._ , 2006, 38, 159–192.
  * Bhushan and Jung 2011 B. Bhushan and Y. C. Jung, _Prog. Mater. Sci._ , 2011, 56, 1–108.
  * Breitenbach _et al._ 2018 J. Breitenbach, I. V. Roisman and C. Tropea, _Exp. Fluids_ , 2018, 59, 55.
  * Ashgriz 2011 N. Ashgriz, _Handbook of Atomization and Sprays: Theory and Applications_ , Springer Science & Business Media, 2011.
  * Li and Thoroddsen 2015 E. Li and S. T. Thoroddsen, _J. Fluid Mech._ , 2015, 780, 636–648.
  * Klaseboer _et al._ 2014 E. Klaseboer, R. Manica and D. Y. Chan, _Phys. Rev. Lett._ , 2014, 113, 194501.
  * Chen _et al._ 2011 L. Chen, J. Wu, Z. Li and S. Yao, _Colloids Surf. A_ , 2011, 384, 726–732.
  * Chen _et al._ 2016 L. Chen, E. Bonaccurso, P. Deng and H. Zhang, _Phys. Rev. E_ , 2016, 94, 063117.
  * Langley _et al._ 2020 K. Langley, A. A. Castrejon-Pita and S. T. Thoroddsen, _Soft Matter_ , 2020\.
  * Andreotti and Snoeijer 2020 B. Andreotti and J. H. Snoeijer, _Annu. Rev. Fluid Mech._ , 2020, 52, 285–308.
  * Bartolo _et al._ 2006 D. Bartolo, C. Josserand and D. Bonn, _Phys. Rev. Lett._ , 2006, 96, 124501.
  * Chen _et al._ 2017 L. Chen, L. Li, Z. Li and K. Zhang, _Langmuir_ , 2017, 33, 7225–7230.
  * Hendrix _et al._ 2012 M. H. Hendrix, R. Manica, E. Klaseboer, D. Y. Chan and C.-D. Ohl, _Phys. Rev. Lett._ , 2012, 108, 247803.
  * Pack _et al._ 2019 M. Y. Pack, A. Yang, A. Perazzo, B. Qin and H. A. Stone, _Phys. Rev. Fluids_ , 2019, 4, 123603.
  * Plesset and Prosperetti 1977 M. S. Plesset and A. Prosperetti, _Annu. Rev. Fluid Mech._ , 1977, 9, 145–185.
  * Day _et al._ 1998 R. F. Day, E. J. Hinch and J. R. Lister, _Phys. Rev. Lett._ , 1998, 80, 704.
  * Leppinen and Lister 2003 D. Leppinen and J. R. Lister, _Phys. Fluids_ , 2003, 15, 568–578.
  * Cox 1986 R. Cox, _J. Fluid Mech._ , 1986, 168, 169–194.
  * Vo and Tran 2021 Q. Vo and T. Tran, _arXiv preprint arXiv:2101.02821_ , 2021.
  * Guo _et al._ 2020 J. Guo, S. Zou, S. Lin, B. Zhao, X. Deng and L. Chen, _Phys. Fluids_ , 2020, 32, 122112.
  * Siddique _et al._ 2020 A. U. Siddique, M. Trimble, F. Zhao, M. M. Weislogel and H. Tan, _Phys. Rev. Fluids_ , 2020, 5, 063606.
  * Yamamoto _et al._ 2018 K. Yamamoto, M. Motosuke and S. Ogata, _Appl. Phys. Lett._ , 2018, 112, 093701.
  * Keller and Miksis 1983 J. B. Keller and M. J. Miksis, _SIAM J. Appl. Math._ , 1983, 43, 268–277.
  * Lister and Stone 1998 J. R. Lister and H. A. Stone, _Phys. Fluids_ , 1998, 10, 2758–2764.
  * Biance _et al._ 2004 A.-L. Biance, C. Clanet and D. Quéré, _Phys. Rev. E_ , 2004, 69, 016301.
  * Mitra and Mitra 2016 S. Mitra and S. K. Mitra, _Langmuir_ , 2016, 32, 8843–8848.