# Nonlinear bosonic Maxwell’s demon

Atirach Ritboon  Research Unit in Energy Innovations and Modern Physics
(EIMP), Thammasat University, Khlong Nueng, Khlong Luang, Pathum Thani, 12120,
Thailand.  Department of Optics, Palacky University, 17. listopadu 12, 771 46
Olomouc, Czech Republic<EMAIL_ADDRESS>Radim Filip<EMAIL_ADDRESS>Department of Optics, Palacky University, 17. listopadu 12, 771 46 Olomouc,
Czech Republic.

###### Abstract

Maxwell’s demon principle of extracting valuable resources through measuring
fluctuations in the system already stimulated modern quantum physics. In
contrast to classical physics, a free coupling to a probe and its free
measurement fundamentally shape the system state. This becomes a new dimension
of the Maxwell demon effect, as in addition to the gained information, the
back action on the system can be exploited and essential for further
applications. We investigate quantum bosonic Maxwell’s demon coupled to a two-
level system to address this issue straightforwardly. The deterministic
multiple subtractions of energy quanta by an energetically conservative
Jaynes-Cummings interaction leads to an out-of-equilibrium state. Although
still super-Poissonian, it can resonantly excite another two-level system
better than any thermal state. To further reduce the super-Poissonian
statistics close to a Poissonian by a Maxwell’s demon operation and increase
the excitation rate, we suggest subsequent use of still energetically
conservative multiphonon subtractions performed by an available nonlinear
Jaynes-Cummings interaction. The optimal combination of both deterministic
subtractions leads to statistics that approaches a Poissonian distribution
otherwise produced by shot-noise-limited sources as an ideal laser requiring
extreme bosonic nonlinear saturations.

##  1 Introduction

Maxwell’s demon, a thought experiment, demonstrates that, if ones could access
to the information about the state of a system through a classical
measurement, then they can exploit such information to gain mechanical work or
energy from the system through a classical control over it. This thought
experiment leads to the generalization of the second law of thermodynamics by
emphasizing the possibility of information-work conversion. It is one of the
vital principles that rectify thermal fluctuations, without using strong
nonlinearity, simply by measurement and classical control [1]. In classical
Maxwell’s demon, the classical measurement is ideally arbitrarily precise
without any back action on the system, as a measured quantity of a system is
treated as a hidden variable. It changes dramatically in the quantum domain,
as the different couplings to a probe and its subsequent measurement form new
states of the system. Such events often turn out to be destructive; however,
they sometimes can conditionally distill the system into a more useful
resource [2, 3], as broadly explored in quantum information, especially in
resource theories. If the distillation fails, it can still be repeated
continuously until obtaining the resource. In this way, the chance of failure,
in principle, can be minimized to zero at the cost of the speed or protocol
multiplexing. The initial thermal fluctuations are used at the output with a
negligible probability. Practically, it does not significantly influence the
generated resource. The quantum transformation to a better resource,
therefore, depends on the coupling and measurement back action, not only on
information like in an ideal classical measurement where any back action is
left unconsidered [4, 5, 6], and it is regarded as the key point of this work.
First, we need to define free states, free couplings, and free measurements
[7, 8]. The free states are naturally thermal equilibrium states diagonal in
the energy eigenbasis. The essential free unitary coupling of the system to
the probe is then energy-conserving; the free probe energy measurement
commutes with the probe’s energy. Free controlled operations are also energy
conserving couplings with an ancilla in thermal states that maximally have the
energy of the input ones. It defines the most basic but nontrivial playground
to explore and compare Maxwell demon methods mutually and also with other
noise rectification strategies.

The quantum Maxwell demon is more involved and diverse for a bosonic system
representing a single mode of photons, phonons or other bosonic particles.
Here, the simplest case of a free coupling is an energy-conserving beam-
splitter type of resonant coupling. After this beam-splitter coupling,
macroscopic measurement integrating energy already allows conditional
manipulation with continuous energy statistics, for example, used in [9, 10].
Microscopic single-quanta detection opens space for subtracting individual
energy quanta conditionally [11, 12, 13, 15, 16, 17, 18, 19], even for
macroscopic thermal states, and charging the macroscopic battery by average
energy [20, 21].

However, for microscopic phononic states with few quanta on average, the
statistics after subtraction becomes crucial for charging a microscopic
battery. Such battery is represented by a two-level system coupled to the
phonons, light or microwave fields. Multiple subtractions increase mean energy
and reduce autocorrelation between quanta causing them to be more
statistically independent [14, 15]. They mainly increase the mean-to-deviation
ratio of the system’s energy, which is essential for information theory and
thermodynamics [22]. A deterministic, repeat-until-success version also
deterministically excites a two-level system beyond any thermal state with
stimulating thermodynamical consequences [23]. Moreover, as recently
demonstrated, the correlations between two thermal baths allow a Maxwell’s
demon based protocol to extract more work [24] and measurement strategies have
been used in quantum memristors [25].

In this work, we propose a nonlinear bosonic Maxwell’s demon working at
quantum level through simple and deterministic protocol, which is expected be
straightforwardly realized in various quantum platforms. We first investigates
this deterministic Maxwell demon method for a broadly feasible energy-
conserving coupling, a linear Jaynes-Cummings coupling [26], probing a bosonic
system sequentially by two-level systems to reach an out-of-equilibrium state.
We then prove that the output state can excite another two-level system better
than any thermal states. Differently from photonic Maxwell’s demon, we
consider phononic systems represented, for example, by the extensively used
mechanical modes of a single atom [27] or, recently, a macroscopic oscillator
[28]. Alternatively, microwave superconducting experiments can also be
considered for the experimental tests [29]. In these cases, usually, the mean
thermal occupation per mode can be much higher than that of thermal light
sources. Despite the low dimension of the probes, the deterministic linear
subtraction increases both the mean-to-deviation ratio of energy and the
probability of exciting two-level systems higher than that from thermal
states. It proves the power of such operations beyond a conventional Fock
state lowering [30]. We further involve still energy-conserving nonlinear
Jayness-Cumming coupling [31, 32, 33] available at trapped-ion platforms [35],
cavity quantum electrodynamics [36] and superconducting circuits [37]
mentioned above to perform a nonlinear subtraction of more quanta at once.
Trilinear interactions, additionally, can also considered as alternative
options [38, 39, 40, 41]. Remarkably, we prove that optimally implementing
nonlinear subtractions after the linear ones increase both the mean-to-
deviation ratio of energy and the probability of excitation of atoms. It is
the first example of using two-quanta processes to bring the statistics of
phononic mode closer to a Poissonian, without any classical external drive and
intense nonlinear saturation typical for such processes in laser [35]. It
proves that Maxwell’s demon, based on available nonlinear energy-conserving
couplings, can open a new territory for quantum statistical and
thermodynamical investigations.

We study the deterministic enhancement of population inversion of two-level
systems coupled with harmonic oscillators via the well-known Jaynes-Cummings
(JC) interaction. This improvement is originated from the successive change in
the population distribution of the harmonic oscillators as a result of
sequential subtractions of the oscillators’ excitation. If the subtractions
fail, we classically replace the oscillators with thermal oscillators or ones
obtained from the previous round of subtraction to keep the overall procedure
deterministic. The phonon distribution, thereby, is modified into a bell-
shaped so that the end-product oscillator can be used as a resource for better
exciting a two-level system.

##  2 Results

Figure 1: The diagram represents the overview of the deterministic nonlinear
bosonic protocol for resonantly charging a quantum battery, represented by a
two-level system. The protocol is divided into two parts: phonon subtractions
and the two-level charging process. The excitation subtraction also consists
of two parts: linear (green) and nonlinear (blue) subtractions using the
protocol II, illustrated in the bottom inset and described in detail in the
later part of section 3. Several linear subtractions are performed to
gradually shape the probability distribution of a harmonic oscillator,
initially in a thermal state $\rho_{\rm th}$, into a bell-shaped before
nonlinear subtractions take place to trim and squeeze the probability
distribution further. Such harmonic oscillators, hence, are used to charge the
quantum batteries, initially being set to be in the ground state $|{\rm
g}\rangle$, through Jaynes-Cummings interaction. The excitation probability
$p_{\rm e,out}$ is expected to exceed its corresponding thermal bound.

To understand the overall picture and procedure of this work, we devote this
section to explain the overview of the proposed protocol employed to gradually
shape the probability distribution of a harmonic oscillator into a bell-shaped
for a better probability of exciting a two-level system, which is regarded as
a quantum battery, via JC interaction. Figure 1 displays the overall processes
of the scheme. A harmonic oscillator in thermal equilibrium with a thermal
bath undergoes a linear excitation subtraction using a linear JC coupling for
several times before further extracting its excitations through the nonlinear
interaction so that its probability distribution becomes even more squeezed
from both sides, at low and high quanta. The output harmonic oscillator is
then used to charge a two-level battery through the linear JC coupling to
examine if its performance exceeds the thermal bound, the maximum probability
of exciting a qubit with a thermal state $\rho_{\rm th}$, the initial state it
associated with.

The lower inset depict the excitation subtraction procedure of both linear and
nonlinear subtractions. A harmonic oscillator in a motional state $\rho^{\rm
mo}_{N}$ resonantly interacts with a two-level system through a JC
interaction, with the interaction Hamiltonian linear in oscillator variables,
or a nonlinear JC interaction [32] having higher powers of the oscillator
variables in the interaction Hamiltonian depending on which type of excitation
subtraction being performed at that stage. After they interact for a
strategically chosen period of time, the state of the qubit is measured. If
the measurement results the excited state, the subtraction is thus performed
successfully, otherwise it fails. The measured outcome then feeds forward to
decide weather the harmonic oscillator could be kept or replaced by its
previous successful version before repeating every steps again in the next
round of subtraction.

##  3 Linear subtraction

Figure 2: Protocol I: the diagram illustrates the proposed phonon subtraction
protocol. In the $N$th round of subtraction, an oscillator in a state
$\rho^{\rm mo}_{N}$ couples with a two-level system at its ground state for
the optimal interaction time $t^{\rm op}_{N}$. Then a measurement in the
energy basis of the qubit is performed. If the excited state $|{\rm e}\rangle$
is obtained, we keep the oscillator (blue), as its energy is successfully
subtracted, for the next round, otherwise replacing it with a new oscillator
in the initial thermal state.

The linear subtraction is performed using a resonant JC interaction whose
Hamiltonian in the interaction picture can be written as

$\displaystyle\hat{H}_{\rm
int}=\hbar\lambda\left(\hat{\sigma}_{+}\hat{a}+\hat{\sigma}_{-}\hat{a}^{\dagger}\right),$
(1)

where $\lambda$ is the coupling strength, $\hat{\sigma}_{+}$
($\hat{\sigma}_{-}$) is the rising (lowering) operator of the qubit, and
$\hat{a}^{\dagger}$ ($\hat{a}$) is the creation (annihilation) operator of the
oscillators. We then can express the unitary operator associated with this
coupling running for an interval time $t$ as

$\displaystyle\hat{U}_{\rm JC}(t)=$ $\displaystyle|{\rm e}\rangle\langle{\rm
e}|\cos(\lambda t\sqrt{\hat{n}+1})+|{\rm g}\rangle\langle{\rm g}|\cos(\lambda
t\sqrt{\hat{n}})$ $\displaystyle\quad-{\rm
i}\hat{\sigma}_{-}\hat{a}^{\dagger}\frac{\sin(\lambda
t\sqrt{\hat{n}+1})}{\sqrt{\hat{n}+1}}-{\rm
i}\hat{\sigma}_{+}\hat{a}\frac{\sin(\lambda
t\sqrt{\hat{n}})}{\sqrt{\hat{n}}},$ (2)

where $|{\rm g}\rangle\langle{\rm g}|$ and $|{\rm e}\rangle\langle{\rm e}|$
denote the projections on the ground and excited states of the qubits
respectively, and $\hat{n}=\hat{a}^{\dagger}\hat{a}$ is the number operator of
the harmonic oscillators.

The first proposed protocol, denoted as protocol I, for linearly subtracting
the motional excitations is schematically illustrated in figure 2. At the
beginning, harmonic oscillators are in thermal equilibrium with a thermal bath
at temperature $T$, while the probes, two-level systems, are prepared in the
ground state $|{\rm g}\rangle$. The initial composite state of them thus can
be expressed as

$\displaystyle\rho_{0}=|{\rm g}\rangle\langle{\rm g}|\otimes\rho_{\rm th},$
(3)

where $\rho_{\rm th}$ denotes the state of a harmonics oscillator in thermal
equilibrium with a mean number of excitations $\bar{n}$,

$\displaystyle\rho_{\rm th}$
$\displaystyle=\sum^{\infty}_{m=0}\frac{\bar{n}^{m}}{(\bar{n}+1)^{n+1}}|{m}\rangle\langle{m}|$
$\displaystyle\equiv\rho^{\rm mo}_{0},$ (4)

and is also regarded as the initial motional state of the oscillators.

The mean number of motional excitations $\bar{n}$ is related to the
temperature $T$ by $\bar{n}=(\exp(\hbar\omega/k_{\rm B}T)-1)^{-1}$, where
$\omega$ is the angular frequency of the oscillators. The interaction between
the qubits and the oscillators is run for the optimal time, $t^{\rm op}_{0}$,
chosen to maximize the probability of exciting the qubits, approximately
related to $\bar{n}$ as $\lambda t^{\rm op}_{0}\approx\pi/(2\sqrt{\bar{n}+1})$
(see the appendix for more details). Subsequently, the measurement in the
eigenbasis $\\{|{\rm g}\rangle,|{\rm e}\rangle\\}$ on these qubits is
performed. We then postselect only those oscillators with the probes in the
excited state $|{\rm e}\rangle$ to be used in the further steps of the
protocol. The measurement and postselection project the state of the qubits
and the harmonic oscillators onto

$\displaystyle\rho^{\prime}_{1}(t^{\rm op}_{0})=$
$\displaystyle\frac{1}{P^{(0)}_{\rm e}(t^{\rm
op}_{0})}\left(\hat{\sigma}_{+}\hat{a}\frac{\sin(\lambda t^{\rm
op}_{0}\sqrt{\hat{n}})}{\sqrt{\hat{n}}}\rho_{0}\frac{\sin(\lambda t^{\rm
op}_{0}\sqrt{\hat{n}})}{\sqrt{\hat{n}}}\hat{a}^{\dagger}\hat{\sigma}_{-}\right)$
$\displaystyle=$ $\displaystyle|{\rm e}\rangle\langle{\rm
e}|\otimes\sum^{\infty}_{m=0}\frac{\bar{n}^{m+1}}{(\bar{n}+1)^{m+2}}\frac{\sin^{2}(\lambda
t^{\rm op}_{0}\sqrt{m+1})}{P^{(0)}_{\rm e}(t^{\rm
op}_{0})}|{m}\rangle\langle{m}|$ $\displaystyle=$ $\displaystyle|{\rm
e}\rangle\langle{\rm e}|\otimes\rho^{\rm mo,e}_{0},$ (5)

where $P^{(0)}_{\rm e}(t^{\rm op}_{0})$ is the probability of observing the
excited state $|{\rm e}\rangle$ at the optimal time, which acts as the
normalization factor of the term in the bracket and $\rho^{\rm e}_{0}$ is the
state of the oscillators after the postselection. The qubits of those
postselected systems are then reset back to their ground state $|{\rm
g}\rangle$ by the dissipation of their energy to the environment. We consider
that the qubit dissipation is much slower than the JC interaction and,
therefore, do not decohere the process. After that, we replace those failed
systems with new systems in the initial state $\rho_{0}$. All mentioned
processes are then repeated again, but this time the initial motional state of
the ensemble for the new round has changed from $\rho^{\rm mo}_{0}$ due to the
measurement back action of the first subtraction. For the $N$th round of
repeat-until-success subtraction by the JC interaction, the state of the
ensemble can be expressed as

$\displaystyle\rho^{\rm mo}_{N}=\hat{a}\frac{\sin(\lambda t^{\rm
op}_{N-1}\sqrt{\hat{n}})}{\sqrt{\hat{n}}}\rho^{\rm mo}_{N-1}\frac{\sin(\lambda
t^{\rm
op}_{N-1}\sqrt{\hat{n}})}{\sqrt{\hat{n}}}\hat{a}^{\dagger}+(1-P^{(N-1)}_{\rm
e}(t^{\rm op}_{N-1}))\rho^{\rm mo}_{0},$ (6)

where $\rho^{\rm mo}_{N-1}$ is the achieved motional state of the previous
subtraction, $P^{(N-1)}_{\rm e}$ is the probability of getting the excited
state $|{\rm e}\rangle$ in the previous round, and $t^{\rm op}_{N-1}$ is the
interaction time that maximizes $P^{(N-1)}_{\rm e}$. Note that the motional
state in Eq. 6 when $N=1$ differs from the state in Eq. 3 by additional terms,
associated with the repeat-until-success subtractions as studied in [42]. The
probability $P^{(N)}_{\rm e}$ of getting excited state at the optimal
interaction time $t^{\rm op}_{N}$ can be written as

$\displaystyle P^{N}_{\rm e}(t^{\rm
op}_{N})=\sum^{\infty}_{m=0}p_{m}^{(N)}\sin^{2}(\lambda t^{\rm
op}_{N}\sqrt{m}).$ (7)

where the probability distribution $p_{m}^{(N)}$ of the harmonic oscillators
in the state $\rho_{N}$ can be expressed as

$\displaystyle p^{(N)}_{m}=p^{(N-1)}_{m+1}\sin^{2}(\lambda t^{\rm
op}_{N-1}\sqrt{m+1})+(1-P^{(N-1)}_{\rm e}(t^{\rm op}_{N-1}))p^{(0)}_{m},$ (8)

We denote $p^{(N-1)}_{m}$ to be the achieved probability distribution in the
previous round with $p^{(0)}_{m}=\bar{n}^{m}/(\bar{n}+1)^{m+1}$, the initial
probability distribution at thermal equilibrium. This equation describes how
each subtraction gradually shapes the probability distribution of the
oscillators in each round.

Let us consider the semiclassical case when the average excitation is very
large, $\bar{n}\gg 1$. The square of the sine function, the first term in Eq.
(8), acts as a population filter. After the first subtraction, $N=1$, with the
optimal time $\lambda t^{\rm op}_{0}\approx\pi/(2\sqrt{\bar{n}+1})$, the
probabilities $p^{(1)}_{n}$ with the value of $n$ very different from the
initial average excitation $\bar{n}$ are suppressed as the values of
$\sin^{2}(\pi\sqrt{m+1}/(2\sqrt{\bar{n}+1}))$ is considerably smaller than
unity, while those probabilities $p^{(1)}_{m}$ with $m$ close to $\bar{n}$
dominate the new probability distribution. This filtering effect still holds
true for the subsequent subtractions, and it also makes the optimal
interaction times of several further rounds are approximately the same as the
first one: $\lambda t^{\rm op}_{N}\approx\lambda t^{\rm
op}_{0}\approx\pi/(2\sqrt{\bar{n}+1})$ for $N\ll\bar{n}$. As the probability
of having the qubit in the excited state grows progressively with the number
of performed subtractions, as shown in figure 5, the last term in Eq. (8)
gradually becomes a smaller contribution. Each subtraction of the motional
excitation, with optimized coupling, thus gradually modifies the probability
distribution $p^{(N)}_{m}$ into a bell-shaped centered around $\bar{n}$, as
depicted in figure 4. The center of the probability distribution is shifted to
the left noticeably, when the number of subtractions becomes comparable with
the initial average excitation number $\bar{n}$, indicating that the average
excitation number decreases slightly each time we we perform a subtraction.

Figure 3: Protocol II: the diagram of the improved protocol is similar to that
of protocol I shown in figure 2 but in this case the failed systems are
replaced by the successful systems obtained in the previous round instead of
using systems in the initial thermal state.

Protocol I can be further improved, if in the last step of each subtraction,
instead of replacing the failed systems with systems in thermal equilibrium,
$\rho^{\rm mo}_{0}$, we replace them with the successfully achieved systems of
the previous round. The diagram of the second protocol, named protocol II, is
depicted in figure 3. The achieved state of the ensemble after the $N$th round
of subtraction becomes

$\displaystyle\rho^{\rm mo}_{N}=\hat{a}\frac{\sin(\lambda t^{\rm
op}_{N-1}\sqrt{\hat{n}})}{\sqrt{\hat{n}}}\rho^{\rm mo}_{N-1}\frac{\sin(\lambda
t^{\rm
op}_{N-1}\sqrt{\hat{n}})}{\sqrt{\hat{n}}}\hat{a}^{\dagger}+(1-P^{(N-1)}_{\rm
e}(t^{\rm op}_{N-1}))\rho^{\rm mo}_{N-1},$ (9)

where the population distribution of the harmonic oscillators is modified as

$\displaystyle p^{(N)}_{m}=p^{(N-1)}_{m+1}\sin^{2}(\lambda t^{\rm
op}_{N-1}\sqrt{m+1})+(1-P^{(N-1)}_{\rm e}(t^{\rm op}_{N-1}))p^{(N-1)}_{m}.$
(10)

This modification can suppress both tails of the bell-shaped probability
distribution faster and better than protocol I at the cost of collecting and
storing the outcomes of the previous steps.

Figure 4: The change in the population distributions after five linear
subtractions using (a) protocol I and (b) protocol II is demonstrated. The
black solid lines represent the initial population distribution of the
oscillators being in thermal equilibrium with the mean excitation number
$\bar{n}=30$. Sequential subtractions gradually forms a bell-shaped
probability distributions (orange) with its peak located slightly lower than
the initial mean excitation $\bar{n}=30$. The distribution obtained from
protocol II, in figure 4b, is noticeably narrower than that obtained from
protocol I but has slightly lesser mean excitation, $\langle\hat{n}\rangle$.
These two distributions are then compared to Poisson distributions (blue) with
the same average phonon numbers. Their relevant information are also given,
including their mean excitation $\langle\hat{n}\rangle$, their second-order
correlation functions $g^{2}(0)$, their mean-to-deviation ratios of excitation
$\mathcal{R}$, defined in Eq. 13 and their Fano factors $F$.

We assume that the thermalization time of the motional state is very long
compared to the total time spent in all processes of protocols I and II, so
that the heat transferred from the thermal bath to the considered oscillators
is very small and negligible. The thermalization effect, as a result, can be
ignored.

##  4 Population inversion

A population inversion of qubits happens when the probability $P_{\rm e}$ of
finding the qubits in the excited state $|{\rm e}\rangle$ exceeds the
probability $P_{\rm g}$ of finding them in the ground state $|{\rm g}\rangle$,
_i.e._ $P_{\rm e}-P_{\rm g}>0$ or $P_{\rm e}>1/2$. We then devote this section
to demonstrate and explain the performance of sequential linear subtractions.
Let us first discuss the relation between the population distribution of a
harmonic oscillator and the maximum excitation probability of a two-level
system. For a qubit initially being in the ground state $|{\rm g}\rangle$, the
probability $P_{\rm e}$ of getting the excited state $|{\rm e}\rangle$ after
it is coupled an oscillator via the JC interaction for $t$ is

$\displaystyle P_{\rm e}(t)=\sum^{\infty}_{m=0}p_{m}\sin^{2}(\lambda
t\sqrt{m}),$ (11)

where $p_{m}$ is the population distribution of the oscillator. As an
oscillator being in its motional ground state $|{0}\rangle$ cannot be coupled
with a qubit in the ground state $|{\rm g}\rangle$, the probability $P_{\rm
e}$ thus must be smaller than $1-p_{0}$, where $p_{0}$ is the probability of
finding the oscillator in its ground state. This means the desired probability
distribution should have a small probability $p_{0}$. For a weak coupling case
in which $\lambda t\ll 1$, the probability approaches a linear rule $(\lambda
t)^{2}\langle\hat{n}\rangle$ and the statistics of the oscillator do not
matter in the classical excitation limit. For stronger coupling, however, this
simple approximation breaks. Each term in the summation oscillates in time
with different frequency depending on its index $m$. Narrower probability
distributions thus give constructive interference of the oscillations, as they
cause smaller mis-match between the oscillating frequencies,
$\lambda\sqrt{m}$, of the dominant probabilities $p_{m}$ and provide a higher
chance of getting the excited state $|{\rm e}\rangle$. For example, the
perfect scenario in which the excited state $|{\rm e}\rangle$ is obtained via
the JC interaction for certain is when the oscillator is in an arbitrary Fock
state $|{n}\rangle$ as we can just choose the interaction time $t$ precisely
to match $\lambda t\sqrt{n}=\pi/2$, leading to $P_{\rm e}=1$. Another factor
to be considered for a bell-shaped distribution, with probabilities $p_{m}$
falling shapely when being far from the peak, is the mean excitation number,
$\langle\hat{n}\rangle$. If two distributions have an identical bell-like
shape but different mean excitation numbers, the one with a larger mean
excitation can give larger $P_{\rm e}$. A larger $\langle\hat{n}\rangle$
provides smaller mis-match between the oscillating frequencies
$\lambda\sqrt{m}$ of the terms in the summation. When $m$ differs from the
mean excitation number $\langle\hat{n}\rangle$ by $\Delta m$ such that
$m=\langle\hat{n}\rangle+\Delta m$ and $|\Delta m|\ll\langle\hat{n}\rangle$,
the oscillating frequency of the probability $p_{m}$ can be approximated as

$\displaystyle\lambda\sqrt{m}$
$\displaystyle=\lambda\sqrt{\langle\hat{n}\rangle+\Delta m}$
$\displaystyle\approx\lambda\sqrt{\langle\hat{n}\rangle}+\frac{\lambda\Delta
m}{2\sqrt{\langle\hat{n}\rangle}}.$ (12)

The difference between the oscillating frequencies of the dominant
probabilities in the summation of Eq. 11 is, therefore, inversely proportional
to $\sqrt{\langle\hat{n}\rangle}$. Of course, the mean excitation becomes
irrelevant when it comes to the case of an excited Fock state, as demonstrated
earlier that, with a single oscillating term in the summation, $P_{\rm e}=1$
can be obtained for certain regardless of the mean excitation number. However,
we need to bear in mind that the statistics will immediately play a crucial
role once there exists small deviation from Fock states. From these discussed
facts, among the parameters commonly used for analyzing the statistics of
excitation, such as second-order correlation function,
$g^{2}(0)=\langle\hat{a}^{\dagger
2}\hat{a}^{2}\rangle/\langle\hat{a}^{\dagger}\hat{a}\rangle^{2}$, and Fano
factors $F=\langle(\Delta\hat{n})^{2}\rangle/\langle\hat{n}\rangle$, the
appropriate parameter indicating the desirable phonon statistics, motivated by
Eq. 4, should be the mean-to-deviation ratio (MDR) of the population, denoted
by $\mathcal{R}$, which is defined as

$\displaystyle\mathcal{R}=\frac{\langle\hat{n}\rangle}{\sqrt{\langle(\Delta\hat{n})^{2}\rangle}}$
(13)

where
$\langle(\Delta\hat{n})^{2}\rangle=\langle\hat{n}^{2}\rangle-\langle\hat{n}\rangle^{2}$
represents the phonon fluctuation. It is more likely that the atom is excited
better by a phonon with a greater value of $\mathcal{R}$.

Figure 5: The figure demonstrates the increase of the maximum probability
$P^{(N)}_{\rm e}$ of getting the excited state $|{\rm e}\rangle$ through the
JC coupling after the $N$th subtraction, when using (a) protocol I, depicted
in figure 2, and (b) protocol II, shown in figure 3. The probability
$P^{(N)}_{\rm e}$ increases sharply with the initial mean excitation $\bar{n}$
and reaches its plateau around $\bar{n}=5$. Protocol II apparently provides a
better chance of exciting a two level-system compared to protocol I. The
dotted lines are set to be at 0.5 to mark population inversion $P^{(N)}_{\rm
e}>0.5$, while the dash-dotted lines display the excitation probability if the
oscillator was in a coherent state. A thermal state with a very large mean
excitation, $\bar{n}\gg 1$, can excites a qubit with probability of $P_{\rm
e}\approx 0.6411$.

Figure 5 shows the maximum probability of having the excited state $|{\rm
e}\rangle$ increases after each linear subtraction. This increase originates
from the fact that the population distribution becomes narrower after each
subtraction. The population distribution after five subtractions is shown in
figure 4. However, for a small initial mean excitation number, around
$\bar{n}\sim 1$-$2$, the subtractions do not always increase $P^{(N)}_{\rm e}$
as the excitations of the oscillators are almost exhausted, i.e. most
oscillators are in their motional ground state and no longer coupled with the
two-level systems. As expected, since protocol II gives a smaller probability
of being in the motional ground state, $p_{0}$, and a narrower probability
distribution, it then gives higher probabilities $P^{(N)}_{\rm e}$, for $N\geq
2$. From the figure, we can clearly see that the increase of $P^{(N)}_{\rm e}$
gradually becomes saturated, as the value of $\Delta P^{(N)}_{\rm
e}=P^{(N)}_{\rm e}-P^{(N-1)}_{\rm e}$ becomes smaller. Further subtractions
barely increase the excitation probability. It is obvious from the figure that
the saturated value of $P^{(N)}_{\rm e}$ obtained from protocol II is slightly
greater than that from protocol I. This is because protocol II shapes the
distribution in such a way that the probability $p_{0}$ and its neighborhood
become very small, as shown in figure 4b, compared to the distribution
obtained from protocol I.

##  5 Nonlinear subtraction

As we pointed out in figure 4, the distribution still has a long decaying tail
for larger populations resembling the thermal statistics. To remove this
limitation and shape the population distribution even faster and better,
linear subtraction alone is no longer sufficient. From figure 5, the
performance of the linear subtractions eventually will reach its saturation,
but there is still a way to break through it by utilizing a nonlinear
interaction whose interaction Hamiltonian is of the form,

$\displaystyle\hat{H}_{\rm
non}=\hbar\lambda^{\prime}\left(\hat{\sigma}_{+}\hat{a}^{2}+\hat{\sigma}_{-}\left(\hat{a}^{\dagger}\right)^{2}\right),$
(14)

where $\lambda^{\prime}$ denotes the coupling strength of the interaction.
This interaction Hamiltonian would have the same form as the Hamiltonian in
Eq. (1) if the annihilation and creation operators in Eq. (1) were replaced by
their squares, $\hat{a}\rightarrow\hat{a}^{2}$ and
$\hat{a}^{\dagger}\rightarrow\left(\hat{a}^{\dagger}\right)^{2}$. The unitary
operator describing the time evolution of this nonlinear coupling is given as

$\displaystyle\hat{U}_{\rm non}(t)=$ $\displaystyle|{\rm e}\rangle\langle{\rm
e}|\cos\left(\lambda^{\prime}t\sqrt{(\hat{n}+1)(\hat{n}+2)}\right)+|{\rm
g}\rangle\langle{\rm
g}|\cos\left(\lambda^{\prime}t\sqrt{\hat{n}(\hat{n}-1)}\right)$
$\displaystyle\quad-{\rm
i}\sigma_{-}\left(\hat{a}^{\dagger}\right)^{2}\frac{\sin\left(\lambda^{\prime}t\sqrt{(\hat{n}+1)(\hat{n}+2)}\right)}{\sqrt{(\hat{n}+1)(\hat{n}+2)}}-{\rm
i}\sigma_{+}\hat{a}^{2}\frac{\sin\left(\lambda^{\prime}t\sqrt{\hat{n}(\hat{n}-1)}\right)}{\sqrt{\hat{n}(\hat{n}-1)}}.$
(15)

As mentioned, nonlinear repeat-until-success subtractions are performed in the
same way as linear subtractions explained in section 3, except this time the
employed interaction becomes nonlinear. However, unlike the linear case, a
nonlinear subtraction cannot trim the tails of the population distribution
desirably, see Appendix for more details. Therefore, the prior population
distribution should be in a bell-shaped to some extend, and the probabilities
$p_{m}$ associated with high-energy levels, $m\gg\bar{n}$, must be already
very small. Otherwise, the nonlinear subtraction will form a ripple in the
probability distribution of the harmonic oscillators, which is an undesired
effect. Due to normalization condition, the probability distribution with a
ripple is more dispersed compared to those without it. The probability $P_{\rm
e}$, as a result, is not as large as it potentially should be.

Consequently, in order to use nonlinear subtractions properly, we have to
perform several linear subtractions using protocol II, displayed in figure 3,
so that the population distribution is modified into a proper bell-shaped with
sufficiently short tails. After that, nonlinear subtractions using the
procedure of protocol II can then be performed. The achieved state after a
nonlinear subtraction can be expressed as

$\displaystyle\rho^{\rm f}=\hat{a}^{2}S(\hat{n})\rho^{\rm
i}S(\hat{n})\left(\hat{a}^{\dagger}\right)^{2}+\left(1-P^{\prime}_{\rm
e}(\tau^{\rm op,})\right)\rho^{\rm i},$ (16)

with

$\displaystyle S(\hat{n})=\frac{\sin\left(\lambda^{\prime}\tau^{\rm
op}\sqrt{\hat{n}(\hat{n}-1)}\right)}{\sqrt{\hat{n}(\hat{n}-1)}}$ (17)

where $\rho^{\rm i}$ is the former state of the oscillators before the
nonlinear subtraction, $\tau^{\rm op}$ is the interaction time that gives the
first locally optimal probability of successful subtraction. The probability
of successful subtraction, on the other hand, reads

$\displaystyle P^{\prime}_{\rm e}(t)=\sum^{\infty}_{m=2}p^{\rm
i}_{m}\sin^{2}\left(\lambda^{\prime}t\sqrt{m(m-1)}\right),$ (18)

where $p^{\rm i}_{m}$ is the prior probability distribution of the harmonic
oscillators, the diagonal elements of $\rho^{\rm i}$. In contrast to linear
subtractions, the suitably chosen interaction time for nonlinear subtractions
is not the time that optimizes the probability $P^{\prime}_{\rm e}$. From the
equation, it is easy to notice that the optimal interaction time is
approximately at $\lambda^{\prime}t=\pi/2$, resulting in
$\sin^{2}\left(\lambda^{\prime}t\sqrt{m(m-1)}\right)\approx 1$ for $m>1$.
However, with this interaction time, a nonlinear subtraction almost does not
modified the probability distribution at all as the sine function barely shows
its influence. We instead need to choose the interaction time $\tau^{\rm op}$
that gives the first locally optimal $P^{\prime}_{\rm e}$ such that
$\lambda^{\prime}\tau^{\rm op}<\pi/2$. The interaction time $\tau^{\rm op}$
for this case is approximately related to the mean phonon number
$2\langle\hat{n}\rangle$ as

$\displaystyle\lambda^{\prime}\tau^{\rm
op}\approx\frac{\pi}{2\langle\hat{n}\rangle}.$ (19)

This means we can run a nonlinear interaction even faster than a linear one. A
nonlinear subtraction using protocol II manipulates the probability
distribution as

$\displaystyle p^{\rm f}_{m}=p^{\rm
i}_{m+2}\sin^{2}\left(\lambda^{\prime}\tau^{\rm
op}\sqrt{(m+2)(m+1)}\right)+(1-P^{\prime}_{\rm e}(\tau^{\rm op}))p^{\rm
i}_{m}.$ (20)

Figure 6 compares the probabilities of success and their change after each
subtraction of two different schemes: the scheme that employ only linear
subtractions (the upper bar chart) and the scheme that use linear subtractions
followed by nonlinear subtractions (the lower bar chart) to boost the
performance before charging the quantum battery. From the upper bar chart,
when only linear subtractions being used, the sequential increase of the
probability of success gradually reaches its saturation. Therefore, additional
linear subtractions just barely improves the probability of success and the
charging performance, denoted by the last green bar. The lower bar chart, on
the other hand, demonstrates that the saturated performance can be further
boosted with the help of nonlinear subtractions, whose probabilities of
success are represented by the two red bars. The success probability of a
nonlinear subtraction is noticeably lower than that of the previous linear
subtractions, but it is still sufficiently large to make the nonlinear
subtraction protocol practical. After six linear and two nonlinear
subtractions, the charging performance or the probability of getting the
excited state through the Jaynes-Cumming coupling, denoted by the green bar,
becomes even larger than its previous version when all eight subtractions are
linear. This increase lies in the change in the shape of the population
distribution after the nonlinear subtractions. As shown in the inset of figure
6, the distribution becomes more squeezed which is more desirable for exciting
a two-level system.

Figure 6: The upper bar chart represents the probabilities of success,
obtaining the excited state, in $N$th subtractions using protocol II with an
initial average phonon number of $\bar{n}=30$. On the other hand, the lower
bar chart also shown such probabilities, but this time the 7th and 8th
subtractions are performed by two consecutive nonlinear subtractions, with
their probabilities of success denoted with the red bar. The actual
probabilities of success, $P_{\rm e}$, are written in white on these bars so
that we can see their fractional differences. Above each bar, a probability
distribution of the harmonic oscillator contributing to the probability of
success is displayed together with the value $\mathcal{R}$ of the mean-to-
deviation ratio (MRD) of phonons. The last green bars of both bar charts
represent the probabilities of getting the excited state through the JC
coupling after eight subtractions, which can be regarded as the charging
performance. To compare the performance of both schemes, we mark the height of
the upper green bar on the lower one with the white dashed line. In the inset,
the probability distributions of the two schemes of subtraction are compared.
The solid black line represents the distribution obtained from eight linear
subtraction while the orange histogram shows the distribution obtained when
the two nonlinear subtractions are introduced. The later is then compared with
a Poissonian distribution of the same mean excitation $\langle\hat{n}\rangle$,
displayed by the blue histogram.

We note here that a nonlinear subtraction cannot be used with the procedure of
protocol I, as at the end of the protocol the failed systems are replaced by
systems in a thermal state, which makes high-energy populations not
sufficiently small. As a result, nonlinear subtractions cause the population
distribution to be even more dispersed and a small ripple in the distribution
to form.

##  6 Discussion

Figure 7: The bar chart illustrates the probabilities $(P_{\rm e})^{100}$ of a
hundred out of a hundred quantum batteries being successfully charged after
eight phonon subtractions have been performed on an oscillator in a thermal
state with $\bar{n}=30$ using different strategies including all eight
subtractions being linear (the blue bar) and different combinations of six
linear and two consecutive nonlinear subtractions. The labels $n$th+$(n+1)$th
identify nonlinear subtractions in the sequence of subtractions for each
strategy. It is apparent that the later nonlinear subtractions being performed
the better the probability of success we achieve. The numbers above these bars
are the actual values of such probabilities $(P_{\rm e})^{100}$.

The previous section compared the charging performance obtained from the two
subtraction strategies: the linear-subtractions-only strategy and the
combination of linear and nonlinear subtractions. It emerged that the later
provides a better charging performance compared to the first. The remaining
question is whether should nonlinear subtractions be performed at an earlier
stage or at the very end of a subtraction sequence, or even something in
between, to get the optimal charging performance. To answer this, we then
examined different combinations of linear and nonlinear subtractions, where
the coupling is optimized for each subtraction depending on the previous
measurement outcomes. The result turns out to be that the later the nonlinear
subtractions take place the greater the performance can be. The probabilities
of successful charging an individual battery associated with these strategies
of subtraction may, at first, look insignificantly different, but at the scale
of mass production the differences become eventually magnified. This fact is
demonstrated in figure 7 by comparing the probabilities of successfully
charging a hundred out of a hundred quantum batteries using different
combinations of linear and nonlinear subtractions. To relate it with the
result in section 5, we consider only the cases in which eight subtractions,
including six linear and two consecutive nonlinear subtractions, are performed
before the charging stage and compare their charging performances to that of
the linear-subtractions-only strategy. From the figure, the best performance
is obtained if the last two subtractions are nonlinear. Remarkably, its
increase is even several time larger than the order of the excitation
probability $(P_{\rm e})^{100}$ obtained through the only-linear-subtraction
case. The underlying reason originates from the fact that nonlinear
subtractions are better at squeezing the phonon population but improper for
trimming its tails. Earlier use of nonlinear subtractions, as a result,
reduces the mean phonon number of the end-product motional state and leads to
a lower mean-to-deviation ratio $\mathcal{R}$. Several linear subtractions
thus prepare a properly trimmed phonon distribution to be squeezed by the
following nonlinear subtractions.

##  7 Conclusions

The idea of classical Maxwell’s demon initiates the reformation of the
classical thermodynamics and generates a connection between information and
thermodynamic work, and provides the fundamental idea of the conversion
between these two quantities. With the information of a system obtained
through measurements and precise control, the system then can be manipulated
into an out-of-equilibrium state, and, in return, its energy can be extracted.
The idea of such conversion is carried on to its quantum version with some
fundamental differences. In contrast to the classical case, in which a
measurement is treated to be arbitrarily sharp without any back action on the
measured system, in the quantum domain, both measurement outcomes and their
back actions unavoidably affect the way we conditionally control and
manipulate the system in order to generate a useful resource. Each performed
measurement not only extracts the system’s information but also transforms its
state accordingly.

We have proposed a simple but deterministic protocol to realize a bosonic
Maxwell’s demon at a quantum level exploiting the free coupling between the
mechanical modes of a single atom and its internal electronic state. A
measurement on the qubit with its outcome implying absorption of phonons by
the qubit is regarded as phonon subtraction. It is shown that linear
subtractions from both protocols I and II transforms the phonon state from an
initial thermal state into an out-of-equilibrium state with a bell-shaped
phonon distribution. This transformed motional state can eventually be used to
charge a microscopic battery, another qubit, by exciting it through a linear
JC coupling. The charging performance of such out-of-equilibrium states is
higher than that of its initial thermal state, which can be indicated by its
increased mean-to-deviation ratio, $\mathcal{R}$. The performance is enhanced
each time a subtraction being performed but becomes saturated eventually. To
break through this limitation, a nonlinear subtraction, using a nonlinear JC
coupling to absorb more phonons at once, must be exploited using the procedure
of subtraction protocol II. The nonlinear interaction can boost the charging
performance further, at the cost of its speed in the repeat-until-success
protocols. It can further squeeze the phonon distribution better than the
linear version, which increases $\mathcal{R}$ as a result. Nonetheless, it
still has a drawback as it cannot trim the tails of the phonon population
properly, making it better used as a final performance booster. Although, we
use a trapped ion as an example quantum platform, this proposed protocol can
also be realized easily in other platforms in which a nonlinear JC coupling is
available, such as superconducting circuits and cavity quantum
electrodynamics.

Optimized nonlinear subtractions after several linear subtractions can shape
the end-product phonon statistics closer to, but not yet reach, a Poissonian
without the help of an external classical drive. We, therefore, believe that
nonlinear-based Maxwell’s demon can potentially pave the way for a new area of
theoretical and experimental research in quantum statistics and
thermodynamics.

##  8 Author contributions

R.F. conceived and supervise the project. A.R. and R.F. jointly build theory
and analyze the results. A.R. performed numerical simulations and wrote the
paper with feedback from R.F..

## Acknowledgements

A.R. acknowledges funding support from the NSRF via the Program Management
Unit for Human Resources & Institutional Development, Research and Innovation
(grant number B05F650024). R.F. acknowledges support of project 22-27431S of
the Czech Science Foundation. The research was further supported by Thammasat
University Research Unit in Energy Innovations and Modern Physics (EIMP) and
also the European Union’s 2020 research and innovation programme (CSA-
Coordination and support action, H2020-WIDESPREAD-2020-5) under grant
agreement No. 951737 (NONGAUSS).

## References

  * [1] K. Maruyama, F. Nori, and V. Vedral, Colloquium : The physics of Maxwell’s demon and infomation, Rev. Mod. Phys. 81, 1 (2009).
  * [2] C.H. Bennett, H.J. Bernstein, S. Popescu and B. Schumacher, Concentrating partial entanglement by local operations. Phys. Rev. A 53, 2046–2052 (1996).
  * [3] P.G. Kwiat, S. Barraza-Lopez, A. Stefanov and N. Gisin, Experimental entanglement distillation and ‘hidden’ non-locality, Nature 409, 1014 (2001).
  * [4] J. Parrondo, J. Horowitz, and T. Sagawa, Thermodynamics of information, Nature Phys. 11, 131 (2015).
  * [5] S. Toyabe, T. Sagawa, M. Ueda, E. Muneyuki, and M. Sano, Experimental demonstration of information-to-energy conversion and validation of the generalized Jarzynski equality, Nature Phys. 6, 988 (2010).
  * [6] M. Debiossac, D. Grass, J. J. Alonso, E. Lutz, and N. Kiesel, Thermodynamics of continuous non-Markovian feedback control, Nature Comm. 11, 1360 (2020).
  * [7] E. Chitambar and G. Gour, Quantum resource theories, Rev. Mod. Phys. 91, 025001 (2019).
  * [8] F.G.S.L. Brandão, M. Horodecki, J. Oppenheim, J.M. Renes, and R.W. Spekkens, Resource Theory of Quantum States Out of Thermal Equilibrium, Phys. Rev. Lett. 111, 250404 (2013).
  * [9] T. Sh. Iskhakov, V. C. Usenko, U. L. Andersen, R. Filip, M. V. Chekhova, and G. Leuchs, Heralded source of bright multi-mode mesoscopic sub-Poissonian light, Optics Letters 41, 2149-2152 (2016).
  * [10] T.Sh. Iskhakov, V.C. Usenko, R. Filip, M.V. Chekhova, and G. Leuchs, Low-noise macroscopic twin beams, Phys. Rev. A 93, 043849 (2016).
  * [11] V. Parigi, A. Zavatta, M. Kim, and M. Bellini, Probing quantum commutation rules by addition and subtraction of single photons to/from a light field, Science 317, 1890 (2007).
  * [12] A. Zavatta, V. Parigi, M. S. Kim, and M. Bellini, Subtracting photons from arbitrary light fields: experimental test of coherent state invariance by single-photon annihilation, New J. Phys. 10, 123006 (2008).
  * [13] I. A. Fedorov, A. E. Ulanov, Y. V. Kurochkin, and A. I. Lvovsky, Quantum vampire: collapse-free action at a distance by the photon annihilation operator, Optica 2, 112 (2015).
  * [14] S. M. Barnett, G. Ferenczi, C. R. Gilson, and F. C. Speirits, Statistics of photon-subtracted and photon-added states, Phys. Rev. A 98, 013809 (2018).
  * [15] Y. I. Bogdanov, K. G. Katamadze, G. V. Avosopiants, L. V. Belinsky, N. A. Bogdanova, A. A. Kalinkin, and S. P. Kulik, Multiphoton subtracted thermal states: Description, preparation, and reconstruction, Phys. Rev. A 96, 063803 (2017).
  * [16] K. G. Katamadze, G. V. Avosopiants, Y. I. Bogdanov, and S. P. Kulik, How quantum is the quantum vampire effect?: testing with thermal light, Optica 5, 723 (2018).
  * [17] K. G. Katamadze, E. V. Kovlakov, G. V. Avosopiants, and S. P. Kulik, Direct test of the quantum vampire’s shadow absence with use of thermal light, Opt. Lett. 44, 3286 (2019).
  * [18] K. G. Katamadze, G. V. Avosopiants, N. A. Bogdanova, Yu. I. Bogdanov, and S. P. Kulik, Multimode thermal states with multiphoton subtraction: Study of the photon-number distribution in the selected subsystem, Phys. Rev. A 101, (2020).
  * [19] G. Enzian, J. J. Price, L. Freisem, J. Nunn, J. Janousek, B. C. Buchler, P. K. Lam, and M. R. Vanner, Single-phonon addition and subtraction to a mechanical thermal state, Phys. Rev. Lett. 126, 033601 (2021).
  * [20] M. D. Vidrighin, O. Dahlsten, M. Barbieri, M. S. Kim, V. Vedral, and I. A. Walmsley, Photonic Maxwell’s demon, Phys. Rev. Lett. 116, 050401 (2016).
  * [21] A. Shu, J. Dai, and V. Scarani, Power of an optical Maxwell’s demon in the presence of photon-number correlations, Phys. Rev. A 95, 022123 (2017).
  * [22] J. Hloušek, M. Ježek, and R. Filip, Work and information from thermal states after subtraction of energy quanta, Sci. Rep. 7, 13046 (2017).
  * [23] J. Hloušek, T. Denzler, V. Švarc, M. Ježek, E. Lutz, and R. Filip, Experimental realization of a quantum photonic Maxwell demon, will be on Arxive.
  * [24] G. L. Zanin, M. Antesberger, M. J. Jacquet, P. H. S. Ribeiro, L. A. Rozema, and P. Walther, Enhanced photonic Maxwell’s demon with correlated baths, Quantum 6, 810 (2022).
  * [25] M. Spagnolo, J. Morris, S. Piacentini, M. Antesberger, F. Massa, A. Crespi, F. Ceccarelli, R. Osellame, and P. Walther, Experimental photonic quantum memristor, Nat. Photon. 16, 318 (2022).
  * [26] D. Leibfried, R. Blatt, C. Monroe, and D. Wineland, Quantum dynamics of single trapped ions, Rev. Mod. Phys. 75, 281 (2003).
  * [27] Y. Chu, P. Kharel, T. Yoon, L. Frunzio, P.T. Rakich and R.J. Schoelkopf, Creation and control of multi-phonon Fock states in a bulk acoustic-wave resonator, Nature 563, 666 (2018).
  * [28] B.W. Shore and P.L. Knight, The Jaynes-Cummings Model, J. Mod. Opt. 40, 1195 (1993).
  * [29] J. M. Fink, M. Göppl, M. Baur, R. Bianchetti, P. J. Leek, A. Blais and A. Wallraff, Climbing the Jaynes-Cummings ladder and observing its nonlinearity in a cavity QED system, Nature 454, 315-318 (2008).
  * [30] Mark Um, Junhua Zhang, Dingshun Lv, Yao Lu, Shuoming An, Jing-Ning Zhang, Hyunchul Nha, M. S. Kim and Kihwan Kim, Phonon arithmetic in a trapped ion system, Nature Communications 7, 11410 (2016).
  * [31] C. V. Sukumar and B. Buck, Multi-phonon generalisation of the Jaynes-Cummings model, Phys. Lett. 83A, 211 (1981).
  * [32] S. Singh, Field statistics in some generalized JaynesCummings models, Phys. Rev. A 25, 3206 (1982).
  * [33] Celso J. Villas-Boas and Daniel Z. Rossatto, Multiphoton Jaynes-Cummings Model: Arbitrary Rotations in Fock Space and Quantum Filters, Phys. Rev. Lett. 122, 123604 (2019).
  * [34] K.C. McCormick, J. Keller, S.C. Burd, D.J. Wineland, A.C. Wilson and D. Leibfried, Quantum-enhanced sensing of a single-ion mechanical oscillator, Nature 572, 86 (2019).
  * [35] M.O. Scully and M.S. Zubairy, Quantum Optics, Cambridge University Press; 1st edition (1997).
  * [36] M. Brune, J. M. Raimond, P. Goy, L. Davidovich and S. Haroche, Realization of a two-photon maser oscillator, Phys. Rev. Lett. 59, 1899.
  * [37] L. Garziano, R. Stassi, V. Macrì, A. F. Kockum, S. Savasta and F. Nori, Multiphoton quantum Rabi oscillations in ultrastrong cavity QED, Phys. Rev. A 92, 063830 (2015).
  * [38] S. Ding, G. Maslennikov, R. Hablützel, H. Loh, and D. Matsukevich, Quantum parametric oscillator with trapped ions, Phys. Rev. Lett. 119, 150404 (2017).
  * [39] S. Ding, G. Maslennikov, R. Hablützel, and D. Matsukevich, Cross-Kerr nonlinearity for phonon counting, Phys. Rev. Lett. 119, 193602 (2017).
  * [40] G. Maslennikov, S. Ding, R. Hablützel, J. Gan, A. Roulet, S. Nimmrichter, J. Dai, V. Scarani, and D. Matsukevich, Quantum absorption refrigerator with trapped ions, Nat. Commun. 10, 202 (2019).
  * [41] S. Ding, G. Maslennikov, R. Hablützel, and D. Matsukevich, Quantum simulation with a trilinear Hamiltonian, Phys. Rev. Lett. 121, 130502 (2018).
  * [42] P. Marek, J. Provazník and R. Filip, Loop-based subtraction of a single photon from a traveling beam of light, Opt. Express 26, 29837 (2018).
  * [43] P. L. Knight and P. M. Radmore, Quantum revivals of a two-level system driven by chaotic radiation, Phys. Lett 90A 342 (1982).

##  Appendix A Optimal interaction times

In semi-classical treatment $\bar{n}\gg 1$, the probability of a qubit in its
ground state being excited by interacting with a thermal oscillator having a
large mean phonon number $\bar{n}\gg 1$ as can be approximated as [43]

$\displaystyle P_{\rm e}(t)\approx\lambda t\sqrt{\bar{n}}D(\lambda
t\sqrt{\bar{n}})$ (21)

where $\lambda$ is the coupling strength between the qubit and phonons, and
$D(x)$ is Dawson’s integral, defined by

$\displaystyle D(x)=e^{-x^{2}}\int^{x}_{0}e^{x^{\prime 2}}{\rm d}x^{\prime}.$
(22)

Before reaching the quantum-revival region, this semi-classical approximation
agrees well with the calculation obtained from rigorous quantum treatment. The
probability $P_{\rm e}$ is, therefore, maximized when $\lambda
t\sqrt{\bar{n}}\approx 1.502$.

On the other hand, a thermal oscillator with small mean phonon number,
$\bar{n}\ll 1$ dominantly occupies in its ground and the first excited states,
$|0\rangle$ and $|1\rangle$. The population of all other motional states
becomes negligible compared to that of the two states. It is the phonon in the
Fock state $|{1}\rangle$ that mainly excites a qubit. The probability of
successful excitation of a qubit in this case then becomes

$\displaystyle P_{\rm e}(t)=\sum^{\infty}_{m=0}p_{m}\sin^{2}(\lambda
t\sqrt{m})\approx\bar{n}\sin^{2}(\lambda t\sqrt{1})+\mathcal{O}(\bar{n}^{2}).$
(23)

This imply that the maximum $P_{\rm e}$ occurs when $\lambda t\approx\pi/2$.
To compromise between these two extreme limits, $\bar{n}\gg 1$ and $\bar{n}\ll
1$, we then approximate the optimal interaction time as $\lambda t^{\rm
op}_{0}\approx\pi/(2\sqrt{\bar{n}+1})$, which approaches $\pi/2$ for small
$\bar{n}$ and still well agrees with the semi-classical treatment for large
$\bar{n}$. We note here that this approximation still hold true even in the
intermediate limit of the average phonon number $\bar{n}$. For example, when
$\bar{n}=2$, the approximation of $\lambda t^{\rm
op}_{0}\approx\pi/(2\sqrt{\bar{n}+1})$ is differed from its actual value only
by $\sim 1\%$.

For the semi-classical case, $\bar{n}\gg 1$, after several linear subtractions
the population distribution of phonons is modified into a bell-shaped with its
peak centered around its average phonon number $\langle\hat{n}\rangle$ and
can, therefore, be qualitatively approximated as a Gaussian distribution as

$\displaystyle
p_{m}\sim\frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{(m-\langle\hat{n}\rangle)^{2}}{2\sigma^{2}}\right),$
(24)

where $\sigma^{2}$ represents the variance of the distribution. The optimal
interaction times for linear subtractions for the approximated Gaussian
distribution is then of the form $\lambda t^{\rm
op}_{N}\approx\pi/(2\sqrt{\langle\hat{n}\rangle+1})$, which is similar to the
previous result. On the other hand, for this bell-shaped distribution, the
probability associated with a successful nonlinear subtraction becomes

$\displaystyle P^{\prime}_{\rm e}\approx\frac{1-e^{-2\lambda^{\prime
2}t^{2}\sigma^{2}}\cos(2\lambda^{\prime}t(\langle\hat{n}\rangle))}{2}.$ (25)

With a sufficiently small variance $\sigma^{2}$, we therefore can
qualitatively approximate the optimal interaction time for a nonlinear
subtraction as $\lambda^{\prime}\tau^{\rm
op}\approx\pi/2\langle\hat{n}\rangle$, as shown in Eq. 19.

##  Appendix B Truncation

The approximated optimal interaction times, presented previously, allow us to
analyze how linear and nonlinear phonon subtractions modify the population
distribution. Let us first begin with linear subtractions. We particularly
consider the first terms of Eqs. 8 and 10

$\displaystyle p^{(N)}_{m}=p^{(N-1)}_{m+1}\sin^{2}(\lambda t^{\rm
op}_{N-1}\sqrt{m+1})+{\rm second\ term},$ (26)

which are the dominant terms of the equations. As $\lambda t^{\rm
op}_{N-1}\approx\pi/(2\sqrt{\langle\hat{n}\rangle+1})$, it is straightforward
to realize that after a linear subtraction the populations $p_{m}$ of the Fock
states $|{m}\rangle$ with $m\sim 4k^{2}\langle\hat{n}\rangle$, where $k$ is an
arbitrary integer, are suppressed the most due to the small values of the sine
function, as

$\displaystyle\sin^{2}(\lambda t^{\rm
op}_{N-1}\sqrt{4k^{2}\langle\hat{n}\rangle+1})$
$\displaystyle\approx\sin^{2}\left(\frac{\pi}{2}\sqrt{\frac{4k^{2}\langle\hat{n}\rangle+1}{\langle\hat{n}\rangle+1}}\right)$
(27) $\displaystyle\sim\sin^{2}\left(k\pi\right)=0.$ (28)

We then can consider the first two most-suppressed points in the population
distribution, $m=0$ and $m\sim 4\langle\hat{n}\rangle$, as the locations at
which linear subtractions trim the phonon distribution. However, the
population of a thermal phonon being in an energy level
$m>4\langle\hat{n}\rangle$ is already relatively very small at the beginning.
Especially in the semi-classical case $\bar{n}\gg 1$, the estimated
probability of a thermal phonon populating in such high energy levels is
already much less than unity as

$\displaystyle
P(m>4\bar{n})=\sum^{\infty}_{m>4\bar{n}}\frac{\bar{n}^{m}}{(\bar{n}+1)^{m+1}}\sim\exp(-4)=0.018.$
(29)

A ripple in the distribution originated from the sine function in Eq. 26,
which in principle can occur in the region of $m>4\bar{n}$, therefore, becomes
relatively insignificant and can be omitted as the population in such region
is already very small. This means linear subtractions can truncate the tails
of population distribution desirably as they do not create a noticeable ripple
in the distribution.

For nonlinear subtractions, we can analyze its effect on the phonon
distribution in a similar manner. The modification of the distribution due to
a nonlinear subtraction is described by Eq. 20. To focus on the first term,
the dominant term of the equation, we then can rewrite this equation as

$\displaystyle p^{\rm f}_{m}=p^{\rm
i}_{m+2}\sin^{2}\left(\lambda^{\prime}\tau^{\rm
op}\sqrt{(m+2)(m+1)}\right)+{\rm second\ term}$ (30)

For a bell-shaped distribution with $\langle\hat{n}\rangle\gg 1$ and the
approximated optimal time$\lambda^{\prime}\tau^{\rm
op}\approx\pi/2\langle\hat{n}\rangle$, the population in the energy levels
around $n\sim 2k\langle\hat{n}\rangle$, where again $k$ is an arbitrary
integer, are suppressed the most. We then again select the first two most
suppressed points, $m=0$ and $m\sim 2\langle\hat{n}\rangle$, to be the
truncated points of the population distribution for this case. The oscillation
in phonon population due to nonlinear subtractions can still apparently appear
in the range $m>2\langle\hat{n}\rangle$, if the population of phonons being in
a Fock state higher than $2\langle\hat{n}\rangle$ is not yet sufficiently
small beforehand. Nonlinear subtractions then can undesirably generate a
ripple in the phonon population. To explicitly visualize it, let us consider
the case when a nonlinear subtraction is performed after just a single linear
subtraction. The population is not yet shaped properly by the linear
subtraction. We, however, can somehow qualitatively approximated the optimal
interaction time of the nonlinear subtraction as $\lambda^{\prime}\tau^{\rm
op}\sim\pi/2\bar{n}$ for $\bar{n}\gg 1$. The nonlinear subtraction thus will
trim the population around $m=0$ and $m\sim 2\bar{n}$. As the population of
phonons in energy levels $m>2\bar{n}$ is not yet sufficiently small, the
oscillation in the phonon population then becomes apparent as shown in Fig 8a.
On the other hand, as the first highest peak noticeably becomes a bell-shaped,
we may also exploit several linear subtractions to suppress the ripple and
then further squeezes it with several nonlinear subtractions. The main peak
can be even more squeezed than a Possionian as shown in figure 8b. We note
here that the shape of the distribution in the figure may resembles a sub-
Possonian, but its second-order coherence function $g^{(2)}(0)$ implies
otherwise as $g^{(2)}(0)>1$. This is because there is still a very small, but
finite probability that very high energy Fock states are populated.

Figure 8: (a)The phonon population distribution obtained from three phonon
subtractions, a linear subtraction followed by two consecutive nonlinear
subtractions, for the initial thermal state with $\bar{n}=30$. We can clearly
see the ripple of the distribution as the populations of the energy levels
around $m\sim 45$ and $m\sim 80$ are greatly suppressed due to the nonlinear
interaction. (b) The distribution after eliminating the ripple in Figure 8a
using 8 linear subtractions and squeezing the main part of the distribution
further by 4 nonlinear subtractions. Its statistics approaches Poissonian. A
sub-Poissonian is prevented by the small, but finite population of very high
energy Fock states.