scaletikzpicturetowidth[1]

# Fast Optomechanical Photon Blockade

Yuxun Ling  Physics Department, Blackett Laboratory, Imperial College London,
Prince Consort Road, SW7 2BW, United Kingdom  Sofia Qvarfort  Nordita, KTH
Royal Institute of Technology and Stockholm University, Hannes Alfvéns väg 12,
SE-106 91 Stockholm, Sweden  Department of Physics, Stockholm University,
AlbaNova University Center, SE-106 Stockholm, Sweden  Florian Mintert  Physics
Department, Blackett Laboratory, Imperial College London, Prince Consort Road,
SW7 2BW, United Kingdom  Helmholtz-Zentrum Dresden-Rossendorf, Bautzner
Landstraße 400, 01328 Dresden, Germany

###### Abstract

The photon blockade effect is commonly exploited in the development of single-
photon sources. While the photon blockade effect could be used to prepare
high-fidelity single-photon states in idealized regimes, practical
implementations in optomechanical systems suffer from an interplay of
competing processes. Here we derive a control scheme that exploits destructive
interference of Fock state amplitudes of more than one photon. The resulting
preparation time for photon-blockaded quantum states is limited only by the
optomechanical interaction strength and can thus be orders of magnitude
shorter than in existing schemes that achieve photon-blockaded state in the
steady state.

##  I Introduction

Single-photon sources are important resources in many fields, including
quantum computation [1, 2] and quantum communication [3]. Therefore, the
photon blockade effect [4] has been a research topic of great interest because
it can be utilized to realize single-photon sources [5, 6, 7]. In analogy with
the Coulomb blockade of electrons, the photon blockade effect occurs if a
single photon in a cavity provides an energetic barrier for a second photon to
enter the cavity. This phenomenon has been predicted in several systems, such
as Kerr-type nonlinear cavity [8, 9], cavity-quantum-electrodynamics systems
[5, 10, 11, 12], superconducting circuits [13, 6, 7] and optomechanical
systems [14, 15, 16].

This work specifically considers a generic optomechanical system [17]
consisting of an optical cavity and a mechanical oscillator where the
displacement of the mechanical oscillator is linearly coupled to the intensity
of the light field inside the cavity. The system can be experimentally
realized in several different ways such as suspended mirror [18, 19, 20],
whispering-gallery-mode microresonators [21, 22, 23, 24, 25], levitated
particles [26, 27, 28, 29] and clamped systems [30, 31, 32].

There is substantial literature on the photon blockade effect in
optomechanical systems [14, 33, 34, 35, 36, 37, 38, 39, 15, 16]. However,
existing schemes achieve the photon blockade effect only in the steady state,
where pumping of the cavity and photon decay out of the cavity reach
equilibrium. The time scale required for the preparation of a photon-blockaded
state is thus at least as long as the lifetime of a photon in the cavity. In
particular, in devices with a high-quality cavity, this poses a severe
limitation to the repetition rate at which photon-blockaded quantum states can
be prepared.

In order to exploit the photon blockade on shorter time scales, we explore the
use of a suitably shaped bichromatic driving of the cavity. Rather than
waiting for the occupation of two-photon states (or states with more photons)
to decay, as is typically the case in schemes based on monochromatic driving,
the dynamics resulting from bichromatic driving can lead to destructive
interference of different amplitudes of the two-photon state. We show that
optimized bichromatic driving profiles can be used to realize quantum states
with a value of the equal-time, second-order autocorrelation function that is
three orders of magnitude smaller than what can be achieved with monochromatic
driving.

The paper is structured as follows. In Sec. II, the model and background are
introduced, and in Sec. III the theoretical methods are summarized. Then,
numerical results are shown in Sec. IV. In particular, the autocorrelation
function when a monochromatic driving function is adopted is examined in Sec.
IV.1, where it is shown that the monochromatic driving is insufficient for
strong photon blockade at time scales shorter than the dissipation time scale.
In Sec. IV.2, monochromatic driving is then compared to bichromatic driving
with which strong photon blockade is predicted within the same amount of time.
Crucially, in Sec. IV.3 and Sec. IV.4, it is shown that it is possible to
increase the duration in which the autocorrelation function stays around its
minimum value, as well as increasing the single-photon occupation at the cost
of weaker photon blockade. Finally, limitations in the theory methods are
discussed in Sec. V.

##  II Background

###  II.1 Optomechanical Photon Blockade

In this work, we consider an optomechanical system driven by an external
coherent light source. The bare system can be modelled by the Hamiltonian

$\displaystyle\hat{H}_{OM}$
$\displaystyle=\omega_{c}\hat{a}^{\dagger}\hat{a}+\omega_{m}\hat{b}^{\dagger}\hat{b}-g_{0}\hat{a}^{\dagger}\hat{a}(\hat{b}^{\dagger}+\hat{b})\
,$ (1)

with the cavity (mechanical) resonance frequency $\omega_{c}$ ($\omega_{m}$),
the annihilation and creation operators $\hat{a}$ ($\hat{b}$) and
$\hat{a}^{\dagger}$ ($\hat{b}^{\dagger}$) for the optical cavity field (the
mechanical oscillator) and the single-photon optomechanical coupling strength
$g_{0}$. The driving Hamiltonian reads

$\displaystyle\hat{H}_{d}=\xi^{*}(t)\hat{a}+\xi(t)\hat{a}^{\dagger}\ ,$ (2)

with a complex-valued driving function $\xi(t)$ describing the pumping of the
cavity through an external light source. The temporal shape of the driving
function $\xi(t)$ can be chosen with regular pulse-shaping techniques, and the
central goal of this paper is the optimization of the function $\xi(t)$ for
the realization of single-photon states.

The spectrum of the full system Hamiltonian in Eq.(1) is anharmonic because of
the cubic interaction term (quadratic in optical creation/annihilation
operators and linear in mechanical creation/annihilation operators). This
spectrum is most easily understood in the frame defined by the polaron
transformation [14]

$\displaystyle\hat{U}_{p}=\exp(\frac{g_{0}}{\omega_{m}}\hat{a}^{\dagger}\hat{a}(\hat{b}^{\dagger}-\hat{b}))\
.$ (3)

In this frame, the undriven system Hamiltonian
$\tilde{H}_{OM}=\hat{U}_{p}^{\dagger}\hat{H}_{OM}\hat{U}_{p}$ reads

$\displaystyle\tilde{H}_{OM}=\left(\omega_{c}-\frac{g_{0}^{2}}{\omega_{m}}\hat{a}^{\dagger}\hat{a}\right)\hat{a}^{\dagger}\hat{a}+\omega_{m}\hat{b}^{\dagger}\hat{b}\
.$ (4)

Instead of the optomechanical interaction term, the Hamiltonian now contains a
Kerr-type nonlinearity in the optical subsystem. Since the eigenstates of the
transformed Hamiltonian $\tilde{H}_{OM}$ are given by the direct product of
photonic and phononic Fock states, the spectrum of the transformed Hamiltonian
$\tilde{H}_{OM}$, and thus of the original Hamiltonian $\hat{H}_{OM}$ is given
by

$\displaystyle
E_{n,m}=n\left(\omega_{c}-n\frac{g_{0}^{2}}{\omega_{m}}\right)+m\omega_{m}\ ,$
(5)

where $n$ and $m$ denote the number of photons and phonons respectively.

Crucially for the photon blockade, the optical part of the spectrum is not
evenly spaced. A photon entering an empty cavity is thus resonant with the
system if it has the frequency $\omega_{c}-g_{0}^{2}/\omega_{m}$, while a
second photon entering the cavity that is already occupied with one photon
would need to have the frequency $\omega_{c}-3g_{0}^{2}/\omega_{m}$ to be on
resonance. An optomechanical system driven by a light field of frequency
$\omega_{c}-g_{0}^{2}/\omega_{m}$ thus allows a single photon to enter the
cavity, but energetically forbids a second photon from entering.

While this qualitative picture provides the correct intuition for photon
blockade, the exact mechanism is more involved. In order to understand the
workings in more detail, one also has to appreciate that the polaron
transformation not only makes the undriven system non-interacting but also
modifies the nature of the driving. Although the driving in the laboratory
frame (Eq. (2)) results in the excitation and de-excitation of only photons,
the driving Hamiltonian in the polaron frame

$\displaystyle\hat{U}_{p}^{\dagger}\hat{H}_{d}(t)\hat{U}_{p}=\xi^{*}(t)\hat{a}\exp(\frac{g_{0}}{\omega_{m}}(\hat{b}^{\dagger}-\hat{b}))+h.c.\
,$ (6)

is expressed in terms of photon creation and annihilation operators dressed
with a displacement of the mechanical degree of freedom. The drive can thus
not only (de)excite the quantum mechanical light field in the cavity, but it
can also create or annihilate phonons. The energetic violation in creating a
second photon that is at the core of the photon blockade can thus be
jeopardized by the possibility of depositing excess energy into the mechanical
degree of freedom.

The excitation or de-excitation of individual phonons during the excitation of
a cavity photon is only negligible if the corresponding energetic violation is
large as compared to the driving amplitude, i.e. if
$\absolutevalue{\omega_{m}-2g_{0}^{2}/\omega_{m}}\gg\absolutevalue{\xi(t)}$.

So far, the discussion has assumed an idealized cavity with a perfectly sharp
resonance frequency. However, since any cavity has a finite linewidth
$\kappa$, which can affect the suppression of otherwise energetically
forbidden processes, the additional condition
$\absolutevalue{\omega_{m}-2g_{0}^{2}/\omega_{m}}\gg\kappa$ is required for
the photon blockade effect.

The driving Hamiltonian in the polaron frame (Eq. (6)) includes processes with
the creation or annihilation of any number of phonons, and any such process
can be energetically suppressed if conditions like the ones formulated above
for the creation or annihilation of a single phonon are satisfied. In
particular, to achieve a strong interaction between the cavity field and the
mechanical oscillator, an unrealistically large number of such constraints
needs to be satisfied. It is much more practical to require a weak
interaction, i.e. $g_{0}\ll\omega_{m}$, so that processes including the
creation or annihilation of more than a single phonon in Eq. (6) are
negligible.

The realization of photon blockade thus requires sufficiently strong coupling
$g_{0}$ in order to realize an anharmonic spectrum, but at the same time
sufficiently weak coupling in order to suppress the creation or annihilation
of more than one phonon.

Since the experimental realization of a system with a clear separation of the
relevant scales – (i) the mechanical resonance frequency $\omega_{m}$, (ii)
the Kerr-type nonlinearity $g_{0}^{2}/\omega_{m}$, and (iii) the cavity
linewidth $\kappa$ or the driving amplitude $|\xi(t)|$ – is extremely
challenging, any practical driving scheme needs to work also outside such an
idealized regime. Existing schemes relying on monochromatic driving can
realize a sizeable occupation of the single-photon Fock state and simultaneous
suppression of the two-photon Fock state only in the steady state when the
effect of external driving and system losses have come to an equilibrium. The
necessity to wait for such an equilibration results in a slow process that
limits the efficiency of state preparation.

The goal of this paper is the design of bichromatic driving schemes that help
to realize photon-blockaded states of the light field in short time scales.
The use of bichromatic driving allows us to exploit the interference of
different absorption processes. A suitably chosen driving function $\xi(t)$
can thus lead to destructive interference of different paths towards the
population of the two-photon Fock state. This maximizes the strength of photon
blockade in an optomechanical system.

###  II.2 Photon Blockade and Autocorrelation

It is helpful to identify a quantifier for the strength of photon blockade.
Such a quantifier is the equal-time second-order autocorrelation function,
which is given by

$\displaystyle
g^{(2)}(t)=\frac{\langle\hat{n}^{2}\rangle-\langle\hat{n}\rangle}{\langle\hat{n}\rangle^{2}}\
,$ (7)

with the photon-number operator $\hat{n}=\hat{a}^{\dagger}\hat{a}$ and the
expectation operator $\langle\cdot\rangle$ with respect to the system state
$\hat{\rho}(t)$. Intuitively, the autocorrelation function $g^{(2)}(t)$
describes how likely it is to find more than one photon in the system at the
time $t$. Therefore, the stronger the photon blockade, the lower the value of
the autocorrelation function $g^{(2)}(t)$ is, and for full photon blockade,
the autocorrelation vanishes.

Despite its simple form, the time evolution of the autocorrelation function is
not easy to analytically determine for optomechanical systems. This is because
the system dynamics cannot be solved exactly due to the cubic interaction
term. To address the problem, this work considers the photon blockade effect
in a perturbative regime in which the occupation probability of Fock states
containing three or more photons is negligible as compared to the single-
photon and two-photon occupations. When such higher-photon occupations are
suppressed, the autocorrelation function in Eq. (7) can be expressed as

$\displaystyle
g^{(2)}(t)\approx\frac{2p_{2}(t)}{(p_{1}(t)+2p_{2}(t))^{2}}\mathrel{\mathop{\ordinarycolon}}=\tilde{g}^{(2)}(t)\
,$ (8)

in terms of the occupation probabilities $p_{1}(t)$ and $p_{2}(t)$ of the
single- and two-photon Fock state respectively. In the remainder, it is shown
that the approximate autocorrelation function $\tilde{g}^{(2)}(t)$ can be
analytically determined for optomechanical systems, and the photon blockade
effect are then studied based on the results. As we show in the upcoming
sections, optimal bichromatic driving schemes can help to prepare quantum
states with a value of the exact $g^{(2)}$-function as low as $10^{-4}$ and a
time scale as short as five periods of mechanical motion.

###  II.3 Dissipative Dynamics in Optomechanical Systems

All quantum systems are coupled to their environment, which gives rise to
dissipative processes. There has been impressive progress in the reduction of
mechanical dissipation in many optomechanical experiments [32, 22, 40, 41],
but the necessity to open the cavity to inject light implies a finite rate of
photon loss from the cavity. The following discussion thus takes into account
optical loss but neglects mechanical dissipation. The effects of mechanical
dissipation are explored in Sec. V.

Optical loss for any state $\hat{\rho}$ at rate $\kappa$ can be modelled with
the Lindbladian $\hat{\mathcal{L}}_{c}$ satisfying the relation

$\displaystyle\hat{\mathcal{L}}_{c}\hat{\rho}=\frac{\kappa}{2}\left(2\hat{a}\hat{\rho}\hat{a}^{\dagger}-\left(\hat{a}^{\dagger}\hat{a}\hat{\rho}+\hat{\rho}\hat{a}^{\dagger}\hat{a}\right)\right)\
.$ (9)

The complete Master equation $\dot{\hat{\rho}}=\hat{\mathcal{L}}\hat{\rho}$ of
the system (in the original frame instead of the polaron frame) including
drive and dissipation is then given by the generator $\hat{\mathcal{L}}$ of
the system dynamics that satisfies the equation

$\displaystyle\hat{\mathcal{L}}\hat{\rho}=-i[\hat{H}_{OM}+\hat{H}_{d},\hat{\rho}]+\hat{\mathcal{L}}_{c}(\hat{\rho})\
,$ (10)

with the system Hamiltonian $\hat{H}_{OM}$ as defined in Eq. (1) and the
driving Hamiltonian $\hat{H}_{d}$ as defined in Eq. (2). This can be taken as
a starting point to solve for the time dependence of the photon statistics in
optomechanical systems.

##  III Solving the system dynamics

While steady-state properties can be inferred without solving the system
dynamics, the goal of this work requires an explicit solution of the system’s
equations of motion. The dynamics of the undriven and dissipation-less system
can be solved exactly [42, 43]. However, even though the dynamics of some
observables can be constructed exactly also in the presence of drive [44, 45]
and dissipation [46], this is not the case for the photon occupations of
interest here. It is thus necessary to construct these quantities
perturbatively. Perturbative solutions can then enable the identification of
realizing photon-blockaded states as quickly as possible, which is the goal of
this work.

To solve the dynamics, it is convenient to note that the generator
$\hat{\mathcal{L}}$ in Eq. (10) can be separated into a photon-number-
conserving term $\hat{\mathcal{L}}_{0}$ with

$\displaystyle\hat{\mathcal{L}}_{0}\hat{\rho}$
$\displaystyle=-i[\hat{H}_{OM},\hat{\rho}]-\frac{\kappa}{2}(\hat{a}^{\dagger}\hat{a}\hat{\rho}+\hat{\rho}\hat{a}^{\dagger}\hat{a})\
,$ (11)

and a term $\hat{\mathcal{L}}_{1}$ satisfying the equation

$\displaystyle\hat{\mathcal{L}}_{1}\hat{\rho}$
$\displaystyle=-i[\hat{H}_{d},\hat{\rho}]+\kappa\hat{a}\hat{\rho}\hat{a}^{\dagger}\
,$ (12)

which captures processes in which the photon number is not conserved. In the
interaction picture, the time-evolution of the system can then be described by
the following equation

$\displaystyle\hat{\rho}(t)=\hat{\mathcal{S}}_{0}(t)\hat{\mathcal{S}}_{1}(t)\hat{\rho}_{0}\
,$ (13)

with the propagator $\hat{\mathcal{S}}_{0}$ describing the dynamics induced by
the photon-conserving term $\hat{\mathcal{L}}_{0}$, the propagator
$\hat{\mathcal{S}}_{1}$ describing the photon-varying processes induced by the
photon-varying term $\hat{\mathcal{L}}_{1}$ in the frame defined by
$\hat{\mathcal{S}}_{0}$, and $\hat{\rho}_{0}$ being the initial state of the
system.

The dynamics induced by the photon-conserving term $\hat{\mathcal{L}}_{0}$ can
be constructed exactly by solving the differential equation

$\displaystyle\frac{d}{dt}\hat{\mathcal{S}}_{0}(t)\hat{\rho}_{0}=\hat{\mathcal{L}}_{0}\hat{\mathcal{S}}_{0}(t)\hat{\rho}_{0}\
,$ (14)

for any initial state $\hat{\rho}_{0}$ of the system. The full analytic
expression of the propagator $\hat{\mathcal{S}}_{0}$ is lengthy but interested
readers can find it together with its derivation in Appendix. A.1.

The remainder of the dynamics must be solved perturbatively. The generator
$\tilde{\mathcal{M}}$ of the rest of the system dynamics in the frame defined
by $\hat{\mathcal{S}}_{0}$ reads

$\displaystyle\tilde{\mathcal{M}}(t)=\hat{\mathcal{S}}_{0}^{-1}(t)\hat{\mathcal{L}}_{1}(t)\hat{\mathcal{S}}_{0}(t)\
.$ (15)

The propagator $\hat{\mathcal{S}}_{1}$ induced by the generator
$\tilde{\mathcal{M}}$ can be constructed perturbatively in the driving
strength $\absolutevalue{\xi(t)}$ and the optical decay $\kappa$ using the
Neumann series [47] which read as
$\hat{\mathcal{S}}_{1}=\sum_{k=0}^{\infty}\hat{\mathcal{S}}^{(k)}_{1}$. The
zeroth order term $\hat{\mathcal{S}}^{(0)}_{1}(t)=\hat{\mathbbm{1}}$ is the
identity map, and the $k$-th order term $\hat{\mathcal{S}}^{(k)}_{1}(t)$ (for
$k>0$) reads

$\displaystyle\hat{\mathcal{S}}^{(k)}_{1}(t)=\int_{\bm{t^{(k)}}\in\mathbb{T}}\left(\prod_{j=1}^{k}\tilde{\mathcal{M}}(t_{j})\right)d\bm{t^{(k)}}\
,$ (16)

where $\bm{t^{(k)}}=(t_{1},\ldots,t_{k})$ is the vector of time-ordered time
variables $t_{1}>\ldots>t_{k}$ and $\mathbb{T}$ is the hyperpyramidal domain
$t>t_{1}>\ldots>t_{k}>0$.

Due to the time dependence of the propagator $\hat{\mathcal{S}}_{0}$, it is
generally not possible to evaluate the integrals in Eq. (16) analytically.
However, in order to construct the time-dependence of the probabilities
$p_{1}$ and $p_{2}$ with corrections up to their leading order in
$\absolutevalue{\xi}/\omega_{m}$, only the terms
$\hat{\mathcal{S}}^{(k)}_{1}(t)$ with $k\leq 4$ are required. Furthermore, all
terms in the $\hat{\mathcal{S}}^{(k)}_{1}(t)$ that are at least quartic in
$g_{0}/\omega_{m}$ can be neglected. In these approximations, the
probabilities $p_{1}$ and $p_{2}$ can indeed be constructed analytically for
any driving function $\xi(t)$ that is given as a Fourier series. This is
discussed in more detail in Appendix. A.3.

##  IV Results

We strive in the following for the identification of driving patterns such
that the system evolves towards a state with a minimal value of the
autocorrelation $g^{(2)}$ at a desired point in time. This is achieved using
the analytic dependence of the approximate autocorrelation function
$\tilde{g}^{(2)}$ on the driving profile $\xi(t)$ derived in Sec. III,

Since most of the earlier works focused on monochromatic driving

$\displaystyle\xi_{1}(t)=\epsilon e^{-i(\omega_{c}+\Delta)t}\ ,$ (17)

with driving strength $\epsilon$ and detuning $\Delta$ from the cavity
resonance frequency, this situation is relevant for comparison. The main focus
of the subsequent discussion, however, is on bichromatic driving

$\displaystyle\xi_{2}(t)=\epsilon_{1}e^{-i(\omega_{c}+\Delta_{1})t}+\epsilon_{2}e^{i\psi}e^{-i(\omega_{c}+\Delta_{2})t}\
,$ (18)

with two driving amplitudes $\epsilon_{j}$ and corresponding detunings
$\Delta_{j}$ as well as a phase difference $\psi$ between the two driving
components.

The construction of the optimized driving patterns is based on the system
ground state $\hat{\rho}_{0}=\outerproduct{0}{0}\otimes\outerproduct{0}{0}$ as
the initial state. Even though the approximate autocorrelation function
$\tilde{g}^{(2)}$ is given as an analytic function of all system parameters
including the Fourier components of the driving profile $\xi$, the actual
optimization of the drive cannot be performed analytically. In all the
subsequent examples, numerical optimization is performed with Mathematica
[48].

The validity of all results obtained based on perturbative treatments is
verified through comparison with numerically exact solutions. Numerical
simulations are performed with Qutip [49, 50] in a Hilbert space truncated to
include up to $6$ photons and $15$ phonons.

Since the explicit shape of the driving patterns identified as optimal is not
necessarily insightful on its own, it is not specified directly in the
subsequent discussion. For completeness, however, the optimized numerical
values of the parameters in the driving functions $\xi_{1}$ and $\xi_{2}$ can
be found in Appendix. B.

###  IV.1 Photon Blockade for Monochromatic Driving

For comparison with the results obtained with bichromatic driving discussed
below in Sec. IV.2, this section is focused on the photon-blockaded states
that can be realized with monochromatic driving. In particular, it is can be
found that, within the perturbative regime
($(g_{0}/\omega_{m})^{4},(\epsilon/\omega_{m})^{2}\ll 1$), monochromatic
driving in Eq. (17) is insufficient for strong photon blockade in the time
limit $t\ll 1/\kappa$. Figure 1 plots the values of the autocorrelation
function $g^{(2)}(t_{op})$ from the full numerical simulation as a function of
the target time $t_{op}$. For each data point, the driving detuning $\Delta$
is optimized for minimal approximate autocorrelation function
$\tilde{g}^{(2)}(t_{op})$ at the target time $t_{op}$. In addition, a range of
coupling strength ($g_{0}=0.3\omega_{m},\ldots,0.6\omega_{m}$) is also
considered while other parameters are kept constant
$(\epsilon,\kappa)=(\omega_{m}/200,\omega_{m}/50)$.

The strength of photon blockade depends on the optomechanical coupling
strength $g_{0}$. This can be confirmed by the data points corresponding to
different coupling strength $g_{0}$ for each fixed target time $t_{op}$ in
Fig. 1. The minimized values of the autocorrelation function decrease as the
coupling strength $g_{0}$ increases until it reaches the value
$g_{0}=0.5\omega_{m}$. The main reason for the increase in the autocorrelation
is that the energy difference between different transitions in the unequally
spaced energy spectrum (Eq. (5)) is proportional to the square of the coupling
strength (i.e. $E_{n+1,m}-E_{n,m}\propto g_{0}^{2}$). As the coupling strength
increases further (e.g. $g_{0}\geq 0.6\omega_{m}$, starred points in Fig. 1),
the energy required to pump a single photon into the cavity is eventually
approached by the energy required to pump a second photon and an extra phonon
into the cavity (i.e. $E_{2,1}-E_{1,0}=E_{1,0}-E_{0,0}$), and thus the photon
blockade effect vanishes [14]. This can be confirmed by the increase in the
autocorrelation function for larger coupling strength ($g_{0}>0.5\omega_{m}$).
The increase in the autocorrelation for larger coupling strength is also
consistent with the fact that the assumption of small coupling strength
($(g_{0}/\omega_{m})^{4}\ll 1$) is violated close to the value
$g_{0}=0.5\omega_{m}$, but the errors only contribution to a small portion of
the increase.111The error of the theoretical values relative to the numerical
values increases from $1.4\%$ to $5.1\%$ as the target time $t_{op}$ increases
from $3T$ to $15T$. However, the relative difference between the
autocorrelation for coupling strength with values $g_{0}=0.5$ and $g_{0}=0.6$
ranges from $11\%$ to $32\%$ which are one order of magnitude larger than the
relative difference between analytical results and numerical results.

Furthermore, Fig. 1 also shows that monochromatic driving is insufficient for
strong photon blockade for any target time within the first $15$ mechanical
periods. By looking at the data points corresponding to different target time
$t_{op}$ for each fixed coupling strength $g_{0}$, it can be observed that the
minimized values of the autocorrelation function $g^{(2)}(t_{op})$ decrease
monotonically as the target time $t_{op}$ increases. However, the best value
($\sim 10^{-2}$) of the autocorrelation function that can be obtained within
the first $15$ mechanical periods ($T=2\pi/\omega_{m}$) is still far from the
value ($g^{(2)}\lesssim 10^{-4}$) that can be achieved using a bichromatic
drive.

Figure 1: The value of the autocorrelation function $g^{(2)}$ as a function of
the target time $t_{op}$ in terms of mechanical periods ($T=2\pi/\omega_{m}$)
when the approximate autocorrelation function $\tilde{g}^{(2)}(t_{op})$ is
minimized with respect to the driving detuning $\Delta$ at the considered time
$t_{op}$. Different values of the optomechanical coupling strength
$g_{0}=0.3\omega_{m}$ (cross), $g_{0}=0.4\omega_{m}$ (circle), and
$g_{0}=0.5\omega_{m}$ (triangle) and $g_{0}=0.6\omega_{m}$ (star) are
considered. The driving strength is $\epsilon=\omega_{m}/200$ and cavity decay
rate is $\kappa=0.02\omega_{m}$. The minimized autocorrelation function
decreases when either the driving strength or the evolution time increases,
but strong photon blockade ($g^{(2)}\sim 10^{-4}$) cannot be observed within
the range of parameters considered.

###  IV.2 Fast Photon Blockade for Bichromatic Driving

It is possible to induce a stronger fast photon blockade given an optimized
bichromatic drive. To compare a bichromatic drive with a monochromatic drive,
the time-evolution of the autocorrelation function $g^{(2)}$ is plotted in
Fig. 2 as a function of time $t$. Figure 2(a) shows the autocorrelation
$g^{(2)}$ for the bichromatic drive, while Fig. 2(b) shows the autocorrelation
$g^{(2)}$ for a monochromatic drive. In addition to the autocorrelation
$g^{(2)}$, the single- and two-photon occupation numbers $p_{1}$ and $p_{2}$
are plotted to demonstrate the difference in the time evolution of the optical
subsystem between the two cases. The cavity decay rate is set to
$\kappa=\omega_{m}/50$, as in Sec. IV.1. The coupling strength is fixed to
$g_{0}=0.3\omega_{m}$, which is less than the optimal value ($g_{0}\approx
0.5\omega_{m}$) in Fig. 1.222This is to ensure that the assumption of small
coupling strength ($(g_{0}/\omega_{m})^{4}\ll 1$) is satisfied with an error
less than $0.01$. The driving strength for each component is kept constant as
$\epsilon_{1}=\epsilon_{2}=\omega_{m}/200$. In both plots, the detunings take
the optimized values $\Delta=-0.040\omega_{m}$, $\Delta_{1}=-0.019\omega_{m}$
and $\Delta_{2}=-0.028\omega_{m}$ and the relative phase of the bichromatic
drive takes the optimized value $\psi=2.8$ so that the approximate
autocorrelation $\tilde{g}^{(2)}$ reaches its minimum around $t_{op}=5T$. This
target time was chosen because it is close to the minimal time required to
achieve the autocorrelation value $g^{(2)}(t_{op})\lesssim 10^{-4}$. The
optimized values can be found in Appendix. B.

Figure 2: Single-photon occupation $p_{1}$ (dashed), two-photon occupation
$p_{2}$ (dotted) and autocorrelation function $g^{(2)}$ (solid) as functions
of time for parameters that minimize the autocorrelation function at $5$-th
mechanical period when the parameters have the values
$\epsilon=\omega_{m}/200$, $g_{0}=0.3\omega_{m}$, $\kappa=0.02\omega_{m}$, and
the driving detuning and phases are optimized. Inset (a) corresponds to the
dynamics when the driving is bichromatic and strong photon blockade
($g^{(2)}\sim 10^{-4}$) can be observed at $5T$. Inset (b) corresponds to the
dynamics when the driving is monochromatic and the lowest value of the
autocorrelation function only reaches $10^{-1}$.

From Fig. 2, it can be seen that much lower values of the autocorrelation
$g^{(2)}$ can be achieved with the bichromatic drive. The evolution plotted in
Fig. 2(a) shows a second-order autocorrelation of value $g^{(2)}(5T)\approx
5.8\times 10^{-5}$ at the end of the fifth mechanical period, which is three
orders of magnitude smaller than the best possible value for a monochromatic
driving (see Sec. IV.1) $g^{(2)}(5T)\approx 8.9\times 10^{-2}$ shown in Fig.
2(b). The significant enhancement in the strength of the photon blockade
effect demonstrates the benefit of adopting an optimized bichromatic driving.
Furthermore, it can be observed that the single-photon occupation $p_{1}$ also
increases monotonically in both subplots but the two-photon occupation only
shows a rapid decrease near the $5$-th mechanical period when a bichromatic
driving is adopted. This supports the argument in Sec. II that the entry of
the second drive induces strong destructive interference only between
different transition paths to the two-photon Fock state without affecting the
single-photon occupation. It is also worth noticing that the low
autocorrelation value $g^{(2)}(5T)\approx 5.8\times 10^{-5}$ is achieved given
a coupling strength $g_{0}$ weaker than the optimal value indicated by Fig. 1,
which suggests a further advantage of adopting a bichromatic driving.

It can also be noticed that in Fig. 2, the single-photon population $p_{1}$
for the two subplots are different. However, the comparison of the intracavity
photon numbers is not insightful. This is because the compound driving
strength $\absolutevalue{\xi_{2}(t)}$ (see Eq. (18)) of a bichromatic drive is
dynamic in contrast to the time-independent driving strength
$\absolutevalue{\xi_{1}(t)}=\epsilon$ of a monochromatic drive. Therefore, the
compound driving strength $\absolutevalue{\xi_{2}(t)}$ of a bichromatic drive
can be much weaker than that of a monochromatic drive even if the driving
strengths of the components are equal (i.e.
$\epsilon_{1}=\epsilon_{2}=\epsilon$). On the other hand, the autocorrelation
$g^{(2)}$ does not suffer from the same problem. As long as the driving
strengths for the two driving components are equal and weak (i.e.
$\epsilon_{1}=\epsilon_{2}\ll\omega_{m}$), the value of the approximate
autocorrelation function $\tilde{g}^{(2)}$ to leading order becomes
independent of the driving strength. Therefore, the value of the exact
autocorrelation function $g^{(2)}$ is independent of the driving strength
given the validity of the perturbative method. Further discussion of this
point can be found in Sec. V.3.

Figure 3: The minimal value of autocorrelation function against the relative
cavity decay rate $\kappa/\omega_{m}$ when $t=5T$ (circle) and $t=10T$
(cross). The driving strength is $\epsilon=\omega_{m}/200$ and optomechanical
coupling rate is $g_{0}=0.3\omega_{m}$ The minimal values of the
autocorrelation function increase as the cavity decay increases, but do not
have a monotonic dependence on the evolution time.

Although the driving function is the only degree of freedom to be optimized,
other system parameters also need to be kept within a suitable range for the
strong photon blockade to occur. Especially, the cavity decay rate has a
strong influence on $g^{(2)}$. To explore the behaviour of the
$g^{(2)}$-function as the cavity decay $\kappa$ increases, the minimized
values of the autocorrelation function $g^{(2)}(t_{op})$ are plotted in Fig.
3. Two different evolution time are considered, namely $t_{op}=5T$ (blue
circles) and $t_{op}=10T$ (yellow crosses). The optomechanical coupling is
kept at $g_{0}=0.3\omega_{m}$ and the driving is
$\epsilon_{1}=\epsilon_{2}=\omega_{m}/200$ as before. For each point,
$\Delta_{1}$, $\Delta_{2}$ and $\psi$ have been optimized such that $g^{(2)}$
reaches its lowest value.

From studying Fig. 3, it is clear that the minimal value of the
autocorrelation function $g^{(2)}$ stays below $10^{-4}$ as long as the cavity
linewidth remains below $\kappa\sim 0.03\omega_{m}$. However, as $\kappa$
becomes comparable to the optomechanical coupling $g_{0}$, the unequal
spacings between different energy levels, as shown in Eq. (5), are no longer
distinguishable by the cavity, and thus the photon blockade effect eventually
vanishes for large $\kappa$.

It is also noticeable from Fig. 3 that a longer evolution time does not
necessarily lead to a lower minimal $g^{(2)}$. This is in sharp contrast with
the monochromatic case, as seen in Fig. 2. Such a feature can be explained by
the fact that the occurrence of the minimum of $g^{(2)}$ depends on the
destructive interference between the two driving fields and it is not
guaranteed that a longer evolution time enhances the destructive interference
between the two driving tones.

While it is generally important to minimize the autocorrelation function,
there are other aspects that can be further improved. They are explored in the
following sections.

###  IV.3 Extending the Lifetime of Photon Blockade

The results presented in the previous section show the desired low values of
the autocorrelation function $g^{(2)}$, which indicate strong photon blockade.
Besides strong photon blockade, it is also helpful to extend the duration for
which the photon-blockaded state can be accessed or equivalently, the duration
in which the autocorrelation function stays at its minimum value. In this
section, it is shown that the maximization of this duration can be achieved at
the same time when the strength of photon blockade is maximized.

A broader minimum of the autocorrelation can be achieved by minimizing the
following target function

$\displaystyle\eta_{1}=\biggl{(}1+w_{d}\frac{d}{dt}+w_{s}\frac{d^{2}}{dt^{2}}\biggr{)}\tilde{g}^{(2)}(t_{op})\
,$ (19)

with scalar weights $w_{d}$ and $w_{s}$. The minimization of the first-order
derivative then guarantees that the local minimum of the autocorrelation
function is located in the vicinity of the given time $t_{op}$. A minimal
second-order derivative then assures that the vanishing first-order derivative
stays at a low value around the given time $t_{op}$.

Figure 4 demonstrates that this choice of target function produces the desired
result. A bichromatic driving is optimized for a flat minimum with fixed
weights $w_{d}=1$ and $w_{s}=10$ and the other system parameters are chosen to
be
$(\epsilon_{1},\epsilon_{2},g_{0},\kappa,t_{op})=(\omega_{m}/200,\omega_{m}/200,3\omega_{m}/10,\omega_{m}/50,5T)$
in accordance with Fig. 2. The weight $w_{s}$ for the second-order derivative
is chosen to be larger in order to prioritize the maximization of the width of
the minimum. As can be noted from Fig. 4, the optimized driving in Fig. 4
makes the autocorrelation stay within $\pm 5\%$ of its value at the target
time $t_{op}=5T$ for $0.38$ mechanical periods as compared to $1.1\times
10^{-2}$ mechanical periods in Fig. 2. The improvement is more than one order
of magnitude, which suggests that the target function $\eta_{1}$ in Eq. (19)
correctly predicts the desired driving function that widens the minimum.

There are however drawbacks to the optimization for broader minima. As a
trade-off, the value of the autocorrelation function $g^{(2)}(5T)=1.2\times
10^{-4}$ at is slightly higher as compared to the result in Fig. 2
($g^{(2)}(5T)\approx 5.8\times 10^{-5}$). However, this still suggests strong
photon blockade (i.e. $g^{(2)}\sim 10^{-4}$). The trade-off also indicates
that the increase in the length of time in which the autocorrelation function
stays around its minimum is not unlimited, but rather depends on the threshold
value of the autocorrelation function for strong photon blockade.

Figure 4: Single-photon occupation $p_{1}$ (dashed), two-photon occupation
$p_{2}$ (dotted) and autocorrelation function $g^{(2)}$ (solid) as functions
of time for parameters that minimize the target function Eq. (19) at $5$-th
mechanical period when the parameters have the values
$\epsilon=\omega_{m}/200$, $g_{0}=0.3\omega_{m}$ and $\kappa=0.02\omega_{m}$.
The autocorrelation function stays around its minimum for approximately $0.2$
mechanical periods.

###  IV.4 Optimization for Higher Single-Photon Occupation

While it is important to maximize the strength of the photon blockade effect
by minimizing the $\tilde{g}^{(2)}$ function, maximizing the single-photon
population $p_{1}$ at the same time is also important if a single-photon
source is targeted. This can be done by minimizing the following target
function

$\displaystyle\eta_{2}=\biggl{(}1+w_{d}\frac{d}{dt}+w_{s}\frac{d^{2}}{dt^{2}}\biggr{)}\tilde{g}^{(2)}(t_{op})+\frac{w_{1}}{p_{1}(t_{op})}\
,$ (20)

with scalar weight $w_{1}$. The target function $\eta_{2}$ in Eq. (20) is the
target function $\eta_{1}$ in Eq. (19) plus a second term that maximizes the
single-photon population. Although the multi-photon populations are
automatically minimized whenever the single-photon population approaches
unity, this is not the case in the perturbative regime of weak driving.
Therefore, the autocorrelation function still requires minimization at the
same time when the single-photon population is maximized. The target function
$\eta_{2}$ also assumes that the single-photon population only varies slowly
over the interval in which the autocorrelation function stays near its
minimum. This assumption is justified by the weak driving strength.

For better comparison with the results presented in Sec. IV.3, it is helpful
to assume the system parameters
$(\kappa,g_{0},t_{op})=(\omega_{m}/50,3\omega_{m}/10,5T)$ and the weights
$(w_{d},w_{s})=(1,10)$. The weight $w_{1}$ for the term that is inversely
proportional to the single-photon population $p_{1}$ is taken to be the value
$w_{1}=\epsilon_{1}^{2}/10$ which is proportional to the square strength so
that the term becomes unitless to leading order. By minimizing the target
function $\eta_{2}$ with respect to the driving strengths
($\epsilon_{1},\epsilon_{2}$), the detunings ($\Delta_{1},\Delta_{2}$) and the
phase $\psi$, the single-photon occupation at fifth mechanical period
increases to $p_{1}(5T)\approx 9.4\times 10^{-3}$ as compared to the value
$p_{1}(5T)\approx 6.6\times 10^{-3}$ in Fig. 4. As a trade-off, the value of
the autocorrelation function $g^{(2)}(5T)$ also increases from $1.2\times
10^{-4}$ to $1.6\times 10^{-4}$ but both values are of the same order of
magnitude which suggests strong photon blockade. Similarly, the length of time
at which the autocorrelation function stays within $\pm 5\%$ of the latter
value only slightly reduces to $0.37$ mechanical periods from $0.38$
mechanical periods. Within this interval, the change in the single-photon
population is also less than $5\%$ of its value at the fifth mechanical
period, which supports the slow-varying assumption discussed earlier in the
section.

The photon population in the cavity is ultimately bounded by the maximal
pumping power, and thus the single-photon population can be further enhanced
by increasing the driving strength in addition to minimizing the target
function Eq. (20). For example, with the optimized driving parameters, the
single-photon occupation can be enhanced to exceed the vacuum occupation (i.e.
$p_{1}>0.5$) if the pumping power is increased by $10$-times. However, such
modification can lead to a decrease in the accuracy of the perturbative model
as is discussed in Sec. V.3.

##  V Discussion

In this section, we discuss the assumptions made when deriving the main
result. For example, the interaction between the mechanical oscillator and the
environment can be non-negligible and the perturbative assumptions in the
theory fail when the system is driven too strongly. These differences are
quantified by numerically computing the autocorrelation function and comparing
it with the analytic approximate function $\tilde{g}(t)$ shown in Eq. (8). For
consistency and unless otherwise specified, this section adopts the system
parameters
$(\epsilon_{1},\epsilon_{2},g_{0},\kappa,t_{op})=(\omega_{m}/200,\omega_{m}/200,3\omega_{m}/10,\omega_{m}/50,5T)$
and the driving parameters as used in Fig. 4 for which the autocorrelation
function reaches the value $g^{(2)}(5T)\approx 1.2\times 10^{-4}$ at the end
of the fifth mechanical period and stays within $\pm 5\%$ of this value for
$0.38$ mechanical periods. All numerical simulations in Sec. V are also
performed in a larger Hilbert space including up to $10$ photons and $25$
phonons so that the increasing excitation due to errors can be accurately
calculated.

###  V.1 Thermal Initial States

Thus far, the mechanical initial state has been assumed to be the ground
state. However, undesired thermal excitations in the mechanical oscillator are
common in many experiments. It is therefore helpful to study a thermal initial
state of the mechanical oscillator, which reads

$\displaystyle\hat{\rho}=\frac{1}{1+\bar{n}_{m}}\sum_{n=0}^{\infty}\left(\frac{\bar{n}_{m}}{1+\bar{n}_{m}}\right)^{n}\ket{n}\bra{n}\
,$ (21)

with $\bar{n}_{m}$ being the average phonon number initially in the mechanical
mode. Therefore, in this section, the error of the theoretical model is
analyzed when the mechanical initial state is cooled to a near-ground state
instead of a perfect ground state.

Through numerical simulation, it can be concluded that when the initial state
is cooled to a near-ground state with the phonon number below the threshold
$\bar{n}_{m}<0.1$, the value of the autocorrelation function at the end of the
fifth mechanical period remains at a low value $g^{(2)}(5T)<1.3\times
10^{-4}$. However, for thermal states with higher phonon population
$\bar{n}_{m}=10$, the autocorrelation function $g^{(2)}(5T)$ increases by an
order of magnitude to $3.7\times 10^{-3}$ which suggests a significant
decrease in the strength of the photon blockade. This decrease can be
explained by the larger uncertainty in the position of the oscillator due to
the nature of the thermal state. When this uncertainty becomes comparable to
the displacement of the oscillator induced by the optical field, there is a
non-vanishing probability for another photon to enter the cavity and thus the
strength of the photon blockade effect is diminished.

Although it is not a trivial task to perform a near-ground-state cooling,
there are extensive research works [51, 52, 53] on this topic with less-than-
unity average phonon number $\bar{n}_{b}$.

###  V.2 Mechanical Dissipation

In Sec. IV, we assumed that the interactions between the mechanical oscillator
and the environment are negligible. In any realistic setting, however,
mechanical noise is likely to play an important role. Here, we examine the
impact of thermal noise and dissipation on our results. For optomechanical
systems, mechanical dissipation at a rate $\gamma$ in an environment with mean
thermal occupation $\bar{n}_{b}$ can be modelled by the Lindbladian
$\hat{\mathcal{L}}_{m}+\hat{\mathcal{L}}_{p}$ [54]. The term
$\hat{\mathcal{L}}_{m}$ describes the thermalization of the mechanical mode in
a frame displaced by the optical field and reads

$\displaystyle\hat{\mathcal{L}}_{m}=$
$\displaystyle\gamma(\bar{n}_{b}+1)\hat{\mathcal{D}}[\hat{b}-\frac{g_{0}}{\omega_{m}}\hat{a}^{\dagger}\hat{a}]$
(22)
$\displaystyle+\gamma\bar{n}_{b}\hat{\mathcal{D}}[\hat{b}^{\dagger}-\frac{g_{0}}{\omega_{m}}\hat{a}^{\dagger}\hat{a}]\
,$

where the superoperator $\hat{\mathcal{D}}[\hat{o}]$ satisfying the equation

$\displaystyle\hat{\mathcal{D}}[\hat{o}]\hat{\rho}=\hat{o}\hat{\rho}\hat{o}^{\dagger}-\frac{1}{2}(\hat{o}^{\dagger}\hat{o}\hat{\rho}+\hat{\rho}\hat{o}^{\dagger}\hat{o})$
(23)

is the Lindblad superoperator. Due to the optomechanical coupling, there is
also an extra optical dephasing term

$\displaystyle\hat{\mathcal{L}}_{p}=\frac{4\gamma
g_{0}^{2}}{\omega_{m}^{2}\ln(1+1/\bar{n}_{b})}\hat{\mathcal{D}}[\hat{a}^{\dagger}\hat{a}]$
(24)

induced by the mechanical dissipation.

To determine the effect of mechanical noise on the photon blockade, we plot
the autocorrelation function $g^{(2)}(5T)$ numerically for finite mechanical
loss rate $\gamma$ for different values of the environmental thermal
populations (i.e. $\bar{n}_{b}=0,0.1,1,10$). Figure 5 shows the change in the
autocorrelation function at the end of the fifth mechanical period as a
function of the mechanical decay $\gamma/\omega_{m}$ in units rescaled by the
mechanical frequency. The values of $\gamma$ range from the negligible value
$\gamma=2\times 10^{-6}\omega_{m}$ to $\gamma=2\times 10^{-2}\omega_{m}$ which
equals to the optical decay $\kappa$.

From Fig. 5, we see that a higher dissipation rate $\gamma$ leads to higher
values of autocorrelation $g^{(2)}$, which implies a weaker photon blockade.
In a zero-temperature environment (i.e. $\bar{n}_{b}=0$, indicated by the
crosses), the autocorrelation function appears robust to noise (i.e., it stays
at the same order of magnitude $g^{(2)}\sim 10^{-4}$) even when $\gamma$
becomes comparable to the cavity decay $\kappa=2\times 10^{-2}\omega_{m}$.
This suggests that, in the zero-temperature environment, the standard of
negligible mechanical decay can be met by a variety of experimental setups
[17, 20, 25, 41]. On the other hand, for environmental thermal occupation
larger than the value $\bar{n}_{b}=1$, an extremely low mechanical decay rate
$\gamma\approx 10^{-6}\omega_{m}$ is required for the autocorrelation function
to deviate by less than $10\%$ from its original value. Furthermore, the
autocorrelation function increases by more than one order of magnitude when
the mechanical decay rate $\gamma$ becomes comparable to the optical decay
rate $\kappa$. Although an autocorrelation of value $\sim 10^{-3}$ does not
meet the criterion of strong photon blockade used in this work, it is still
more than one order of magnitude smaller than the best value that can be
achieved with a monochromatic driving as in Sec. IV.1.

Figure 5: The autocorrelation function $g^{(2)}$ as a function of the relative
mechanical decay rate $\gamma/\omega_{m}$ for average phonon number in the
heat bath being $\bar{n}_{b}=0$ (cross), $\bar{n}_{b}=0.1$ (circle),
$\bar{n}_{b}=1$ (triangle), and $\bar{n}_{b}=10$ (square). The dashed line
specifies the value of the autocorrelation function that is $10\%$ larger than
the value without the mechanical loss (i.e. $\gamma=0$). The assumptions of
vanishing mechanical decay rate $\gamma\ll\omega_{m}$ and environmental phonon
occupation $\bar{n_{b}}\approx 0$ are required in order to achieve strong
photon blockade.

The minimum value of the autocorrelation function $g^{(2)}$ has a super-linear
dependence on both the mean phonon occupation $\bar{n}_{b}$ in the heat bath
and the mechanical decay rate $\gamma$. The increasing dependence on the
temperature and the coupling rate of the thermal bath is expected because a
strongly coupled, hot bath can create excitations in the cavity through
optomechanical coupling, which dominates the system dynamics and renders the
optimization invalid.

###  V.3 Validity of Perturbative Methods

The optimization performed on the approximate autocorrelation function
$\tilde{g}^{(2)}$ (Eq. (8)) is only valid as long as the perturbative
solutions presented in Sec. III hold. Their validity is based on the following
two assumptions: the driving strength remains small (i.e.
$(\epsilon/\omega_{m})^{2}\ll 1$), and the coupling strength remains small
(i.e. $(g_{0}/\omega_{m})^{4}\ll 1$). Here, we investigate how the accuracy of
$\tilde{g}^{(2)}$ changes as these two conditions break down.

The points at which the two perturbative assumptions break down can be found
by studying the plot in Fig. 6 where the values of the autocorrelation
function $g^{(2)}(5T)$ and the approximate autocorrelation function
$\tilde{g}^{(2)}(5T)$ are plotted as functions of the driving strength. The
coupling strength $g_{0}$ is kept fixed throughout333This is because a change
in the coupling strength $g_{0}$ requires a re-optimization of the driving
parameters, and the balance between the strength of the autocorrelation
function and the width of the local minimum can be different for each local
minimum of the target function $\eta_{1}$. Therefore, instead of varying the
coupling strength $g_{0}$, the perturbative expansion that assumes a small
coupling strength is truncated to a higher order such that only the assumption
$(g_{0}/\omega_{m})^{12}\ll 1$ is required to hold.  The dashed and dotted
lines show $\tilde{g}^{(2)}$ when terms up to orders $(g_{0}/\omega_{m})^{2}$
and $(g_{0}/\omega_{m})^{6}$ are included, respectively. The circular dots and
the crosses show numerical simulations of $g^{(2)}$ and $\tilde{g}^{(2)}$,
respectively.

Firstly, it can be observed that the theoretical values for
$\tilde{g}^{(2)}(5T)$ (dashed and dotted lines) stay constant in Fig. 6. This
can be explained by the following argument. Because the approximate
autocorrelation function $\tilde{g}^{(2)}$ omits the contribution from the
three and more-photon populations, and because the two-photon population is
negligible near its minimum, the function $\tilde{g}^{(2)}$ is approximately
proportional to the ratio $p_{2}/p_{1}^{2}$. According to the Neumann series
in Eq. (16), the single-photon occupation is dominated by the second-order
term $\hat{\mathcal{S}}^{(2)}_{1}$ which is proportional to $\epsilon^{2}$.
The two-photon occupation is dominated by the fourth-order term
$\hat{\mathcal{S}}^{(4)}_{1}$, which is proportional to $\epsilon^{4}$. Thus,
to leading order, the approximate autocorrelation function is independent of
the driving strength. As a result, the magnitude of $\epsilon$ does not
influence $g^{(2)}$. In contrast, the three-photon occupation is included in
the full numerical calculation of the autocorrelation function $g^{(2)}$, and
thus its value changes as a function of the driving strength.

Next, the numerically obtained values for $g^{(2)}$ and $\tilde{g}^{(2)}$
(circles and crosses) can be compared with the theoretical values (lines). For
a very small driving strength $\epsilon=\omega_{m}/1000$ (Fig. 6 inset), the
majority of the error around the minimum can be attributed to the assumption
of small coupling strength ($(g_{0}/\omega_{m})^{4}\ll 1$) as there is a wider
disagreement between the results when the relevant perturbative expansion is
truncated to different orders (dashed and dotted lines) than between the
theoretical values and the numerical values (lines and scattered data).
However, the overall error of the theoretical model is negligible ($<0.01$)
even when the terms quartic in the coupling strength are truncated. On the
other hand, for stronger driving in Fig. 6, the error is dominated by
inaccuracy due to the truncated Hilbert space in the theoretical model. Both
higher-order correction to the single-photon and two-photon populations (i.e.
$p_{1}$ and $p_{2}$) and the three-photon population can explain the diversion
of the exact autocorrelation values (circular dots) from the theoretical
values (lines). However, the relatively minor difference between the exact
values of the autocorrelation function $g^{(2)}$ and the numerically simulated
values of the approximate autocorrelation function $\tilde{g}^{(2)}$ suggests
that the higher order corrections to the single- and two-photon populations
are responsible for the majority of the error. In the presence of these
errors, the optimized driving parameters may not lead to the true minimum of
the autocorrelation function, but still, lead to a fast decrease in the value
of the autocorrelation function near the specified time $t=5T$.

Overall, according to Fig. 6, it is safe to conclude that the analytic
approximations only start breaking down for values larger than
$\epsilon/\omega_{m}\approx 0.005$, which is the value used in Sec. IV.

Figure 6: The autocorrelation $g^{(2)}$ and $\tilde{g}^{(2)}$ as functions of
the driving strength $\epsilon/\omega_{m}$ rescaled by the mechanical
frequency. The circular data points plot the exact autocorrelation values
$g^{(2)}(5T)$ obtained from a full numerical simulation. The cross data points
plot the approximate autocorrelation values $\tilde{g}^{(2)}(5T)$ obtained
from numerical simulation. The dashed curve plots the case when the
perturbative expansion is truncated up to the order $(g_{0}/\omega_{m})^{2}$
(i.e. $(g_{0}/\omega_{m})^{4}\ll 1$). The dotted curve plots the case when the
relevant perturbative expansion is truncated to a higher order which imposes a
less strict condition (i.e. $(g_{0}/\omega_{m})^{8}\ll 1$) on the coupling
strength. The inset zooms into the results for small driving strengths. The
figure shows that a weak driving strength is necessary for obtaining the
results in this work.

##  VI Conclusions and Outlook

While it is elusive to realize a strong photon blockade effect without the
clear separation of scales, the present control scheme can substantially
improve the fidelity of photon-blockaded quantum states under imperfect
conditions. Since the question of what can be practically achieved with
suitably shaped driving depends on the actual system parameters, a brief
overview over system parameters that have been realized in existing
optomechanical platforms is in order. Table 1 shows the relevant parameters
from a selection of works, together with the calculated autocorrelation
function. The autocorrelation values are computed by minimizing the target
function $\tilde{g}^{(2)}(5T)$ at target time $t_{op}=5T$ with respect to the
driving parameters shown for each case. For each of the examples, two
different temperatures (i.e. $\bar{n}_{b}=0,1$) of the thermal bath
surrounding the mechanical oscillator are shown.

| $g_{0}/\omega_{m}$ | $\kappa/\omega_{m}$ | $\gamma/\omega_{m}$ | $g_{0}/\kappa$ |  $g^{(2)}(5T)$ ($\bar{n}_{b}=0$) |  $g^{(2)}(5T)$ ($\bar{n}_{b}=1$)
---|---|---|---|---|---|---
Nano-acoustic resonator [40] | $1.4\times 10^{-4}$ | $2.3\times 10^{-1}$ | $4.2\times 10^{-11}$ | $6.3\times 10^{-4}$ | $1.0$ | $1.0$
Hybrid system (theory) [55] | $4.0\times 10^{-1}$ | $4.6\times 10^{-2}$ | $1.0\times 10^{-2}$ | $8.7$ | $1.5\times 10^{-3}$ | $1.5\times 10^{-2}$
Parameters used in this work | $3.0\times 10^{-1}$ | $2.0\times 10^{-2}$ | $0$ | $15$ | $5.8\times 10^{-5}$ | $5.8\times 10^{-5}$
Table 1: System parameters in experiments in comparison with the parameters
adopted in this work and the strength of the resulting photon blockade at the
target time $t_{op}=5T$. The coupling strength assumed in this work is too
large in comparison with experimental data in the work on nano-acoustic
resonators. The parameters in the theoretical work [55] suggest that hybrid
systems potentially constitute good platforms for realizing photon blockade.

While it is generally difficult to achieve strong photon blockade, many
experimental proposals are moving towards the values assumed in this work.
Most difficulty originates from the requirement of strong coupling strength
$g_{0}>\kappa$. Especially for solid-state resonators, the best values [32,
22, 40, 41] of the coupling strength are still three orders of magnitude
smaller than the cavity linewidth. As shown in Table 1, the strength of photon
blockade is negligible in such a case. However, there are several theoretical
proposals [56, 57, 55] for enhancing the optomechanical coupling strength. For
example, as listed in Table 1, the coupling strength and the cavity decay rate
in the hybrid optomechanical system [55] are of the same order of magnitude as
considered in this work. Thus, the hybrid system could potentially be a good
platform to realize the results in Sec. IV. Some optomechanical systems with
levitated particles [26, 27, 29, 28] also have strong values of coupling
strength $g_{0}$ as compared to the mechanical frequency $\omega_{m}$.
However, in such a system, it is usually true that the inequality
$g_{0}>\omega_{m}$ holds, which is not consistent with the approximations used
in the present approach.

Even though the present work thus does not directly apply to systems with
strong coupling strength ($g_{0}>\omega_{m}$), one would expect that suitably
designed bichromatic (or polychromatic) driving offers increased capabilities
for the preparation of quantum mechanical states also in systems with stronger
optomechanical interactions. In particular, strong optomechanical interactions
could be used for the preparation of broader classes of quantum states, and
the driving techniques developed in this work can be generalized
straightforwardly to any other quantum state with desirable properties.
Stronger phonon-photon coupling would certainly make it generally more
difficult to find driving patterns such that the system evolves towards the
desired states. However, the potential of the non-Gaussian quantum states that
might be realized in such systems make such an endeavor a worthwhile goal.

##  VII Acknowledgements

We thank Markus Rademacher for helpful discussions. S.Q. is funded in part by
the Wallenberg Initiative on Networks and Quantum Information (WINQ) and in
part by the Marie Skłodowska–Curie Action IF programme “Nonlinear
optomechanics for verification, utility, and sensing” (NOVUS) – Grant-Number
101027183.

## References

  * Lounis and Orrit [2005] B. Lounis and M. Orrit, Single-photon sources, Rep. Prog. Phys. 68, 1129–1179 (2005).
  * Wang  _et al._ [2020] J. Wang, F. Sciarrino, A. Laing, and M. G. Thompson, Integrated photonic quantum technologies, Nat. Photonics 14, 273–284 (2020).
  * O’Brien  _et al._ [2009] J. L. O’Brien, A. Furusawa, and J. Vučković, Photonic quantum technologies, Nat. Photonics 3, 687–695 (2009).
  * Imamoḡlu  _et al._ [1997] A. Imamoḡlu, H. Schmidt, G. Woods, and M. Deutsch, Strongly interacting photons in a nonlinear cavity, Phys. Rev. Lett. 79, 1467–1470 (1997).
  * Faraon  _et al._ [2010] A. Faraon, A. Majumdar, and J. Vučković, Generation of nonclassical states of light via photon blockade in optical nanocavities, Phys. Rev. A 81 (2010).
  * Wang  _et al._ [2016] X. Wang, A. Miranowicz, H.-R. Li, and F. Nori, Multiple-output microwave single-photon source using superconducting circuits with longitudinal and transverse couplings, Phys. Rev. A (Coll. Park.) 94 (2016).
  * Gu  _et al._ [2017] X. Gu, A. F. Kockum, A. Miranowicz, Y.-X. Liu, and F. Nori, Microwave photonics with superconducting quantum circuits, Phys. Rep. 718-719, 1–102 (2017).
  * Huang  _et al._ [2018] R. Huang, A. Miranowicz, J.-Q. Liao, F. Nori, and H. Jing, Nonreciprocal photon blockade, Phys. Rev. Lett. 121 (2018).
  * Ghosh and Liew [2019] S. Ghosh and T. C. H. Liew, Dynamical blockade in a single-mode bosonic system, Phys. Rev. Lett. 123 (2019).
  * Bin  _et al._ [2018] Q. Bin, X.-Y. Lü, S.-W. Bin, and Y. Wu, Two-photon blockade in a cascaded cavity-quantum-electrodynamics system, Phys. Rev. A (Coll. Park.) 98 (2018).
  * Han  _et al._ [2018] Y. F. Han, C. J. Zhu, X. S. Huang, and Y. P. Yang, Electromagnetic control and improvement of nonclassicality in a strongly coupled single-atom cavity-QED system, Phys. Rev. A (Coll. Park.) 98 (2018).
  * Trivedi  _et al._ [2019] R. Trivedi, M. Radulaski, K. A. Fischer, S. Fan, and J. Vučković, Photon blockade in weakly driven cavity quantum electrodynamics systems with many emitters, Phys. Rev. Lett. 122, 243602 (2019).
  * Liu  _et al._ [2014] Y.-X. Liu, X.-W. Xu, A. Miranowicz, and F. Nori, From blockade to transparency: Controllable photon transmission through a circuit-QED system, Phys. Rev. A 89 (2014).
  * Rabl [2011] P. Rabl, Photon blockade effect in optomechanical systems, Phys. Rev. Lett. 107, 063601 (2011).
  * Wang  _et al._ [2015] H. Wang, X. Gu, Y.-X. Liu, A. Miranowicz, and F. Nori, Tunable photon blockade in a hybrid system consisting of an optomechanical device coupled to a two-level system, Phys. Rev. A 92 (2015).
  * Zhu  _et al._ [2018] G.-L. Zhu, X.-Y. Lü, L.-L. Wan, T.-S. Yin, Q. Bin, and Y. Wu, Controllable nonlinearity in a dual-coupling optomechanical system under a weak-coupling regime, Phys. Rev. A (Coll. Park.) 97 (2018).
  * Aspelmeyer  _et al._ [2014] M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, Cavity optomechanics, Rev. Mod. Phys. 86, 1391–1452 (2014).
  * Verlot  _et al._ [2009] P. Verlot, A. Tavernarakis, T. Briant, P.-F. Cohadon, and A. Heidmann, Scheme to probe optomechanical correlations between two optical beams down to the quantum level, Phys. Rev. Lett. 102, 103601 (2009).
  * Gröblacher  _et al._ [2009] S. Gröblacher, K. Hammerer, M. R. Vanner, and M. Aspelmeyer, Observation of strong coupling between a micromechanical resonator and an optical cavity field, Nature 460, 724–727 (2009).
  * Kleckner  _et al._ [2011] D. Kleckner, B. Pepper, E. Jeffrey, P. Sonin, S. M. Thon, and D. Bouwmeester, Optomechanical trampoline resonators, Opt. Express 19, 19708–19716 (2011).
  * Ding  _et al._ [2011] L. Ding, C. Baker, P. Senellart, A. Lemaitre, S. Ducci, G. Leo, and I. Favero, Wavelength-sized GaAs optomechanical resonators with gigahertz frequency, Appl. Phys. Lett. 98, 113108 (2011).
  * Verhagen  _et al._ [2012] E. Verhagen, S. Deléglise, S. Weis, A. Schliesser, and T. J. Kippenberg, Quantum-coherent coupling of a mechanical oscillator to an optical cavity mode, Nature 482, 63–67 (2012).
  * Dong  _et al._ [2012] C. Dong, V. Fiore, M. C. Kuzyk, and H. Wang, Optomechanical dark mode, Science 338, 1609–1613 (2012).
  * Jiang  _et al._ [2015] X. Jiang, M. Wang, M. C. Kuzyk, T. Oo, G.-L. Long, and H. Wang, Chip-based silica microspheres for cavity optomechanics, Opt. Express 23, 27260–27265 (2015).
  * Monifi  _et al._ [2016] F. Monifi, J. Zhang, Ş. K. Özdemir, B. Peng, Y.-X. Liu, F. Bo, F. Nori, and L. Yang, Optomechanically induced stochastic resonance and chaos transfer between optical fields, Nat. Photonics 10, 399–405 (2016).
  * Murch  _et al._ [2008] K. W. Murch, K. L. Moore, S. Gupta, and D. M. Stamper-Kurn, Observation of quantum-measurement backaction with an ultracold atomic gas, Nat. Phys. 4, 561–564 (2008).
  * Schleier-Smith  _et al._ [2011] M. H. Schleier-Smith, I. D. Leroux, H. Zhang, M. A. Van Camp, and V. Vuletić, Optomechanical cavity cooling of an atomic ensemble, Phys. Rev. Lett. 107, 143005 (2011).
  * Millen  _et al._ [2020] J. Millen, T. S. Monteiro, R. Pettit, and A. N. Vamivakas, Optomechanics with levitated particles, Rep. Prog. Phys. 83, 026401 (2020).
  * Samutpraphoot  _et al._ [2020] P. Samutpraphoot, T. Đorđević, P. L. Ocola, H. Bernien, C. Senko, V. Vuletić, and M. D. Lukin, Strong coupling of two individually controlled atoms via a nanophotonic cavity, Phys. Rev. Lett. 124 (2020).
  * Massel  _et al._ [2011] F. Massel, T. T. Heikkilä, J.-M. Pirkkalainen, S. U. Cho, H. Saloniemi, P. J. Hakonen, and M. A. Sillanpää, Microwave amplification with nanomechanical resonators, Nature 480, 351–354 (2011).
  * Teufel  _et al._ [2011] J. D. Teufel, T. Donner, D. Li, J. W. Harlow, M. S. Allman, K. Cicak, A. J. Sirois, J. D. Whittaker, K. W. Lehnert, and R. W. Simmonds, Sideband cooling of micromechanical motion to the quantum ground state, Nature 475, 359–363 (2011).
  * Chan  _et al._ [2011] J. Chan, T. P. M. Alegre, A. H. Safavi-Naeini, J. T. Hill, A. Krause, S. Gröblacher, M. Aspelmeyer, and O. Painter, Laser cooling of a nanomechanical oscillator into its quantum ground state, Nature 478, 89–92 (2011).
  * Stannigel  _et al._ [2012] K. Stannigel, P. Komar, S. J. M. Habraken, S. D. Bennett, M. D. Lukin, P. Zoller, and P. Rabl, Optomechanical quantum information processing with photons and phonons, Phys. Rev. Lett. 109, 013603 (2012).
  * Liao and Nori [2013] J.-Q. Liao and F. Nori, Photon blockade in quadratically coupled optomechanical systems, Phys. Rev. A 88 (2013).
  * Liao and Law [2013] J.-Q. Liao and C. K. Law, Correlated two-photon scattering in cavity optomechanics, Phys. Rev. A 87 (2013).
  * Xu  _et al._ [2013] X.-W. Xu, Y.-J. Li, and Y.-X. Liu, Photon-induced tunneling in optomechanical systems, Phys. Rev. A 87 (2013).
  * Kronwald  _et al._ [2013] A. Kronwald, M. Ludwig, and F. Marquardt, Full photon statistics of a light beam transmitted through an optomechanical system, Phys. Rev. A 87 (2013).
  * Kómár  _et al._ [2013] P. Kómár, S. D. Bennett, K. Stannigel, S. J. M. Habraken, P. Rabl, P. Zoller, and M. D. Lukin, Single-photon nonlinearities in two-mode optomechanics, Phys. Rev. A 87 (2013).
  * Lü  _et al._ [2013] X.-Y. Lü, W.-M. Zhang, S. Ashhab, Y. Wu, and F. Nori, Quantum-criticality-induced strong kerr nonlinearities in optomechanical systems, Sci. Rep. 3, 2943 (2013).
  * MacCabe  _et al._ [2020] G. S. MacCabe, H. Ren, J. Luo, J. D. Cohen, H. Zhou, A. Sipahigil, M. Mirhosseini, and O. Painter, Nano-acoustic resonator with ultralong phonon lifetime, Science 370, 840–843 (2020).
  * Ren  _et al._ [2020] H. Ren, M. H. Matheny, G. S. MacCabe, J. Luo, H. Pfeifer, M. Mirhosseini, and O. Painter, Two-dimensional optomechanical crystal cavity with high quantum cooperativity, Nat. Commun. 11, 3373 (2020).
  * Mancini  _et al._ [1997] S. Mancini, V. I. Man’ko, and P. Tombesi, Ponderomotive control of quantum macroscopic coherence, Phys. Rev. A 55, 3042–3050 (1997).
  * Bose  _et al._ [1997] S. Bose, K. Jacobs, and P. L. Knight, Preparation of nonclassical states in cavities with a moving mirror, Phys. Rev. A 56, 4175–4186 (1997).
  * Latmiral and Mintert [2018] L. Latmiral and F. Mintert, Deterministic preparation of highly non-classical macroscopic quantum states, Npj Quantum Inf. 4 (2018).
  * Ling and Mintert [2021] Y. Ling and F. Mintert, Deterministic preparation of nonclassical states of light in cavity optomechanics, Phys. Rev. Research 3 (2021).
  * Qvarfort  _et al._ [2021] S. Qvarfort, M. R. Vanner, P. F. Barker, and D. E. Bruschi, Master-equation treatment of nonlinear optomechanical systems with optical loss, Phys. Rev. A (Coll. Park.) 104 (2021).
  * Neumann [1877] C. G. Neumann, Untersuchungen über das logarithmische und newtonsche potential, Leipzig: Teubner (1877).
  * [48] Wolfram Research, Inc., Mathematica, Version 13.1, champaign, IL, 2022\.
  * Johansson  _et al._ [2013] J. R. Johansson, P. D. Nation, and F. Nori, QuTiP 2: A python framework for the dynamics of open quantum systems, Comput. Phys. Commun. 184, 1234–1240 (2013).
  * Johansson  _et al._ [2012] J. R. Johansson, P. D. Nation, and F. Nori, QuTiP: An open-source python framework for the dynamics of open quantum systems, Comput. Phys. Commun. 183, 1760–1772 (2012).
  * Yang  _et al._ [2019] J.-Y. Yang, D.-Y. Wang, C.-H. Bai, S.-Y. Guan, X.-Y. Gao, A.-D. Zhu, and H.-F. Wang, Ground-state cooling of mechanical oscillator via quadratic optomechanical coupling with two coupled optical cavities, Opt. Express 27, 22855–22867 (2019).
  * Naseem and Müstecaplıoğlu [2021] M. T. Naseem and Ö. E. Müstecaplıoğlu, Ground-state cooling of mechanical resonators by quantum reservoir engineering, Commun. phys. 4 (2021).
  * Brubaker  _et al._ [2022] B. M. Brubaker, J. M. Kindem, M. D. Urmey, S. Mittal, R. D. Delaney, P. S. Burns, M. R. Vissers, K. W. Lehnert, and C. A. Regal, Optomechanical ground-state cooling in a continuous and efficient electro-optic transducer, Phys. Rev. X. 12 (2022).
  * Hu  _et al._ [2015] D. Hu, S.-Y. Huang, J.-Q. Liao, L. Tian, and H.-S. Goan, Quantum coherence in ultrastrong optomechanics, Phys. Rev. A 91 (2015).
  * Dezfouli  _et al._ [2019] M. K. Dezfouli, R. Gordon, and S. Hughes, Molecular optomechanics in the anharmonic cavity-QED regime using hybrid metal–dielectric cavity modes, ACS Photonics 6, 1400–1408 (2019).
  * Xuereb  _et al._ [2012] A. Xuereb, C. Genes, and A. Dantan, Strong coupling and long-range collective interactions in optomechanical arrays, Phys. Rev. Lett. 109, 223601 (2012).
  * Rimberg  _et al._ [2014] A. J. Rimberg, M. P. Blencowe, A. D. Armour, and P. D. Nation, A cavity-cooper pair transistor scheme for investigating quantum optomechanics in the ultra-strong coupling regime, New J. Phys. 16, 055008 (2014).

##  Appendix A Theory

In this Appendix, details on the solution to the Master equation in Eq. (10)
are presented.

###  A.1 System dynamics

It is helpful to define the left-hand and the right-hand ladder superoperators
$\hat{a}_{L}$, $\hat{b}_{L}$ and $\hat{a}_{R}$, $\hat{b}_{R}$ that satisfy the
equations

$\displaystyle\hat{a}_{L}(\hat{\rho})=\hat{a}\hat{\rho},\qquad\hat{a}_{R}(\hat{\rho})=\hat{\rho}\hat{a},\qquad\hat{b}_{L}(\hat{\rho})=\hat{b}\hat{\rho},\qquad\hat{b}_{R}(\hat{\rho})=\hat{\rho}\hat{b}\
.$ (25)

In terms of the ladder superoperators, the photon-conserving evolution
operator $\hat{\mathcal{S}}_{0}$ defined in Sec. III reads

$\displaystyle\hat{\mathcal{S}}_{0}(t)=\hat{\mathcal{S}}_{F}(t)\hat{\mathcal{S}}_{D}(t)\hat{\mathcal{S}}_{I}(t)\hat{\mathcal{S}}_{K}(t)$
(26)

with

$\displaystyle\hat{\mathcal{S}}_{F}(t)=$
$\displaystyle\exp\left(-i\omega_{c}(\hat{a}^{\dagger}_{L}\hat{a}_{L}-\hat{a}^{\dagger}_{R}\hat{a}_{R})t\right)\exp\left(-i\omega_{m}(\hat{b}^{\dagger}_{L}\hat{b}_{L}-\hat{b}^{\dagger}_{R}\hat{b}_{R})t\right)\
,$ (27) $\displaystyle\hat{\mathcal{S}}_{D}(t)=$
$\displaystyle\exp\left(-\frac{\kappa}{2}(\hat{a}^{\dagger}_{L}\hat{a}_{L}+\hat{a}^{\dagger}_{R}\hat{a}_{R})t\right)\
,$ $\displaystyle\hat{\mathcal{S}}_{I}(t)=$
$\displaystyle\exp\left(f_{1}(t)(\hat{a}_{L}^{\dagger}\hat{a}_{L}\hat{b}_{L}^{\dagger}-\hat{a}_{R}^{\dagger}\hat{a}_{R}\hat{b}_{R})-h.c.\right)\
,$ $\displaystyle\hat{\mathcal{S}}_{K}(t)=$
$\displaystyle\exp\left(if_{2}(t)\left((\hat{a}_{L}^{\dagger}\hat{a}_{L})^{2}-(\hat{a}_{R}^{\dagger}\hat{a}_{R})^{2}\right)\right)\
.$

The propagators $\hat{\mathcal{S}}_{F}(t)$, $\hat{\mathcal{S}}_{D}(t)$,
$\hat{\mathcal{S}}_{I}(t)$, $\hat{\mathcal{S}}_{K}(t)$ correspond to the free
evolution, the dephasing, the optomechanical interaction and the Kerr-type
nonlinearity respectively and

$\displaystyle f_{1}(t)$
$\displaystyle=\frac{g_{0}}{\omega_{m}}(e^{i\omega_{m}t}-1)\ ,$ (28)
$\displaystyle f_{2}(t)$
$\displaystyle=\left(\frac{g_{0}}{\omega_{m}}\right)^{2}(\omega_{m}t-\sin(\omega_{m}t))\
,$

are real, time-dependent functions. In the coherent case without driving and
dissipation (i.e. $\kappa=\xi(t)=0$), this solution coincides with results
from earlier works [42, 43].

The photon-varying part $\hat{\mathcal{L}}_{1}(t)$ of the Master equation in
the frame defined by the photon-conserving superoperator
$\hat{\mathcal{S}}_{0}$ then reads

$\displaystyle\tilde{\mathcal{M}}(t)=\hat{\mathcal{S}}_{0}^{-1}(t)\hat{\mathcal{L}}_{1}(t)\hat{\mathcal{S}}_{0}(t)=-i(\xi^{*}\hat{\mathcal{A}}_{L}-\xi\hat{\mathcal{A}}_{R}+\xi\hat{\mathcal{A}}_{L}^{\dagger}-\xi^{*}\hat{\mathcal{A}}_{R}^{\dagger})+\kappa\hat{\mathcal{A}}_{L}\hat{\mathcal{A}}_{R}\
,$ (29)

in which

$\displaystyle\hat{\mathcal{A}}_{L}=$
$\displaystyle\hat{a}_{L}e^{-(\kappa/2+i\omega_{c})t}e^{if_{2}(2\hat{a}_{L}^{\dagger}\hat{a}_{L}-1)}e^{f_{1}\hat{b}^{\dagger}_{L}-f_{1}^{*}\hat{b}_{L}}\
,$ (30) $\displaystyle\hat{\mathcal{A}}_{L}^{\dagger}=$
$\displaystyle\hat{a}_{L}^{\dagger}e^{(\kappa/2+i\omega_{c})t}e^{-if_{2}(2\hat{a}_{L}^{\dagger}\hat{a}_{L}+1)}e^{f_{1}^{*}\hat{b}_{L}-f_{1}\hat{b}^{\dagger}_{L}}\
,$ $\displaystyle\hat{\mathcal{A}}_{R}=$
$\displaystyle\hat{a}_{R}e^{(-\kappa/2+i\omega_{c})t}e^{-if_{2}(2\hat{a}_{R}^{\dagger}\hat{a}_{R}-1)}e^{f_{1}^{*}\hat{b}^{\dagger}_{R}-f_{1}\hat{b}_{R}}\
,$ $\displaystyle\hat{\mathcal{A}}_{R}^{\dagger}=$
$\displaystyle\hat{a}_{R}^{\dagger}e^{(\kappa/2-i\omega_{c})t}e^{if_{2}(2\hat{a}_{R}^{\dagger}\hat{a}_{R}+1)}e^{f_{1}\hat{b}_{R}-f_{1}^{*}\hat{b}_{R}^{\dagger}}$

are the creation and annihilation superoperators in the frame defined by the
superoperator $\hat{\mathcal{S}}_{0}$. Given any initial state
$\hat{\rho}_{i}$, the final state $\hat{\rho}_{f}$ can thus be perturbatively
expanded using the Neumann series (Eq. (16)) into the form of the following
sum of integrals

$\displaystyle\hat{\rho}_{f}=\sum_{k=0}^{\infty}\int_{\bm{t^{(k)}}\in\mathbb{T}}\left(\hat{\mathcal{S}}_{0}(t)\prod_{j=1}^{k}\tilde{\mathcal{M}}(t_{j})\hat{\rho}_{i}\right)d\bm{t^{(k)}}\
,$ (31)

with $\tilde{\mathcal{M}}$ defined in Eq. (29) and the vector $\bm{t^{(k)}}$
of time-ordered variables $t_{k}<\ldots<t_{1}$, and the hyperpyramidal domain
$\mathbb{T}$ being $0<t_{k}<\ldots<t_{1}<t$ as defined in Eq. (16). This can
be taken as a starting point for the calculation of expectation values.

###  A.2 Dynamics of The Vacuum State

One natural choice of the initial state is the vacuum state
$\hat{\rho}_{i}=\outerproduct{0}{0}_{c}\otimes\outerproduct{0}{0}_{m}$.
Whenever the Neumann series in Eq. (31) is truncated to finite order, the
final state of the system given a vacuum initial state is a finite
superposition

$\displaystyle\hat{\rho}_{f}(t)=\sum_{j=0}^{N}g_{j}(t)\hat{\rho}(n_{L,j},n_{R,j},\beta_{L,j},\beta_{R,j})\
,$ (32)

with $N$ being a finite integer, $g_{j}$ being time-dependent, scalar
functions and the density matrices

$\displaystyle\hat{\rho}(n_{L},n_{R},\beta_{L},\beta_{R})=\outerproduct{n_{L}}{n_{R}}_{c}\otimes\outerproduct{\beta_{L}}{\beta_{R}}_{m}\
,$ (33)

with the optical components $\outerproduct{n_{L}}{n_{R}}_{c}$ in the Fock
basis, and the mechanical states $\outerproduct{\beta_{L}}{\beta_{R}}_{m}$ in
the coherent basis. As a result, the photon statistics of the final state can
be extracted from corresponding components that satisfy the relation
$n_{L}=n_{R}$.

To prove that the final state can be written as a finite sum of matrices
$\hat{\rho}(n_{L},n_{R},\beta_{L},\beta_{R})$, it suffices to show that any
matrix in the form $\hat{\rho}(n_{L},n_{R},\beta_{L},\beta_{R})$ is mapped to
another matrix in the same form when acted on by the superoperators
$\hat{\mathcal{A}}_{L}^{({\dagger})}$, $\hat{\mathcal{A}}_{R}^{({\dagger})}$,
and $\hat{\mathcal{S}}_{0}$. This is because the expression
$\hat{\mathcal{S}}_{0}(t)\prod\tilde{\mathcal{M}}(t_{j})$ in each of the
$k$-th order Neumann term (Eq. (31)) is a finite polynomial of the
superoperators $\hat{\mathcal{A}}_{L}^{({\dagger})}$,
$\hat{\mathcal{A}}_{R}^{({\dagger})}$, and $\hat{\mathcal{S}}_{0}$.

The action of the superoperators on the matrix
$\hat{\rho}(n_{L},n_{R},\beta_{L},\beta_{R})$ reads

$\displaystyle\hat{\mathcal{A}}_{L}\hat{\rho}(n_{L},n_{R},\beta_{L},\beta_{R})$
$\displaystyle=g_{L}\hat{\rho}(n_{L}-1,n_{R},\beta_{L}+f_{1},\beta_{R})\ ,$
(34)
$\displaystyle\hat{\mathcal{A}}_{R}\hat{\rho}(n_{L},n_{R},\beta_{L},\beta_{R})$
$\displaystyle=g_{R}\hat{\rho}(n_{L},n_{R}-1,\beta_{L},\beta_{R}+f_{1}^{*})\
,$
$\displaystyle\hat{\mathcal{S}}_{0}\hat{\rho}(n_{L},n_{R},\beta_{L},\beta_{R})$
$\displaystyle=g_{c}\hat{\rho}(n_{L},n_{R},e^{-i\omega_{m}t}(\beta_{L}+f_{1}n_{L}),e^{i\omega_{m}t}(\beta_{R}+f^{*}_{1}n_{R}))$

with

$\displaystyle g_{L}$
$\displaystyle=\sqrt{n_{L}}e^{-(\kappa/2+i\omega_{c})t}e^{if_{2}(2n_{L}-1)}e^{(f_{1}\beta_{L}^{*}-f_{1}^{*}\beta_{L})/2}\
,$ (35) $\displaystyle g_{R}$
$\displaystyle=\sqrt{n_{R}}e^{(-\kappa/2+i\omega_{c})t}e^{-if_{2}(2n_{R}-1)}e^{(f_{1}^{*}\beta_{R}^{*}-f_{1}\beta_{R})/2}\
,$ $\displaystyle g_{c}$
$\displaystyle=e^{-\kappa(n_{L}+n_{R})t/2}e^{-i\omega_{c}(n_{L}-n_{R})t}e^{n_{L}(f_{1}\beta_{L}^{*}-h.c.)/2}e^{n_{R}(f_{1}^{*}\beta_{R}^{*}-h.c.)/2}e^{f_{2}(n_{L}^{2}-n_{R}^{2})}$

being scalar, time-dependent functions, and thus Eq. (32) is proven. In the
next section, the photon statistics obtained by partial tracing out the
mechanical degree of freedom of the matrix elements satisfying the relation
$n_{L}=n_{R}$ is presented.

###  A.3 Photon Occupation

Since the approximate autocorrelation $\tilde{g}^{(2)}$ in Eq. (8) is a
function of photon occupation $p_{1}$ and $p_{2}$ only, it is sufficient to
calculate only the corresponding matrix elements in the density matrix in Eq.
(32).

In the regime of weak driving ($\absolutevalue{\xi}^{2}\ll\omega_{m}^{2}$),
the leading order single-photon occupation $p_{1}(t)$ at time $t$ appears in
the second order of the perturbative series in Eq. (31) and reads

$\displaystyle
p_{1}(t)=2\real\int_{0}^{t}\int_{0}^{t_{1}}\exp\left(\frac{g_{0}^{2}}{\omega_{m}^{2}}(e^{i\omega_{m}(t_{1}-t_{2})}-1)+i\frac{g_{0}^{2}}{\omega_{m}}(t_{2}-t_{1})-\kappa
t+\frac{\kappa}{2}(t_{1}+t_{2})\right)\zeta(t_{1})\zeta^{*}(t_{2})dt_{2}dt_{1}\
,$ (36)

where $\real$ denotes the real part of the expression and
$\zeta(t)=\xi(t)\exp(i\omega_{c}t)$ being the driving function in the frame
rotating at the cavity resonance frequency.

The leading order contribution to two-photon occupation $p_{2}(t)$ at time $t$
appears in the fourth order of the series in Eq. (31) (i.e.
$\propto\absolutevalue{\xi}^{4}$) and can be written as

$\displaystyle
p_{2}(t)=4\real\int_{\bm{t^{(4)}}\in\mathbb{T}}\exp\left(-\frac{2g_{0}^{2}}{\omega_{m}^{2}}+\frac{\kappa}{2}(t_{1}+t_{2}+t_{3}+t_{4}-4t)\right)\sum_{j=1}^{3}h_{j}(t,\bm{t^{(4)}})d\bm{t^{(4)}}\
,$ (37)

with the vector $\bm{t^{(4)}}=(t_{1},t_{2},t_{3},t_{4})$ of time-ordered
variables and the domain $\mathbb{T}$ as defined in Eq. (16), and the time-
dependent functions $h_{j}(t,\bm{t^{(4)}})$ defined as

$\displaystyle h_{1}(t,\bm{t^{(4)}})=$
$\displaystyle\exp\left(\frac{g_{0}^{2}}{\omega_{m}^{2}}\sum_{1\leq j<k\leq
4}s^{(1,1)}_{j,k}\exp\left(s^{(2,1)}_{j,k}(t_{j}-t_{k})\right)-s^{(1,1)}_{j,k}s^{(2,1)}_{j,k}(t_{j}-t_{k})\right)\zeta^{*}(t_{1})\zeta^{*}(t_{2})\zeta(t_{3})\zeta(t_{4})\
,$ (38) $\displaystyle h_{2}(t,\bm{t^{(4)}})=$
$\displaystyle\exp\left(\frac{g_{0}^{2}}{\omega_{m}^{2}}\sum_{1\leq j<k\leq
4}s^{(1,2)}_{j,k}\exp\left(s^{(2,2)}_{j,k}(t_{j}-t_{k})\right)-s^{(1,2)}_{j,k}s^{(2,2)}_{j,k}(t_{j}-t_{k})\right)\zeta^{*}(t_{1})\zeta(t_{2})\zeta^{*}(t_{3})\zeta(t_{4})\
,$ $\displaystyle h_{3}(t,\bm{t^{(4)}})=$
$\displaystyle\exp\left(\frac{g_{0}^{2}}{\omega_{m}^{2}}\sum_{1\leq j<k\leq
4}s^{(1,3)}_{j,k}\exp\left(s^{(2,3)}_{j,k}(t_{j}-t_{k})\right)-s^{(1,3)}_{j,k}s^{(2,3)}_{j,k}(t_{j}-t_{k})\right)\zeta(t_{1})\zeta^{*}(t_{2})\zeta^{*}(t_{3})\zeta(t_{4})\
,$

and the expressions $s^{(1,i)}_{j,k}$ and $s^{(2,i)}_{j,k}$ are either
positive or minus sign as defined in Table 2.

$(j,k)$ $(l,i)$ | $(1,1)$ | $(1,2)$ | $(1,3)$ | $(2,1)$ | $(2,2)$ | $(2,3)$
---|---|---|---|---|---|---
$(1,2)$ | $-$ | $+$ | $+$ | $+$ | $-$ | $+$
$(1,3)$ | $+$ | $-$ | $+$ | $-$ | $+$ | $+$
$(1,4)$ | $+$ | $+$ | $-$ | $-$ | $+$ | $+$
$(2,3)$ | $+$ | $+$ | $-$ | $-$ | $-$ | $-$
$(2,4)$ | $+$ | $-$ | $+$ | $-$ | $-$ | $-$
$(3,4)$ | $-$ | $+$ | $+$ | $-$ | $-$ | $-$
Table 2: The values of the sign variables $s^{(l,i)}_{j,k}$ in Eq. (38).

The integral in Eq. (37) contains double exponentials in the form of
$\exp(\exp(ix))$ which cannot be integrated into a closed form and numerical
tools are inefficient for multiple integrals. However, given a small coupling
strength that satisfies the inequality $g^{4}_{0}/\omega_{m}^{4}\ll 1$, the
double exponential can be expanded into the following form

$\displaystyle\exp(\frac{g_{0}^{2}}{\omega_{m}^{2}}e^{i\omega_{m}(t_{p}-t_{q})})\approx
1+\frac{g_{0}^{2}e^{i\omega_{m}(t_{p}-t_{q})}}{\omega_{m}^{2}}\ ,$ (39)

using the Taylor expansion. The right-hand side of Eq. (39) can be integrated
exactly but the results are too long to be presented in the text. Interested
readers can find the full results online 444Fast-Photon-Blockade-Codes
(Version 1.0.0) [Github Repository]. github.com/lingyx97/Fast-Photon-Blockade-
Codes.

##  Appendix B Parameters Used for Figures

In this section, the optimized parameters used for figures in the main text
are listed. All parameters are normalized in units of the mechanical frequency
(i.e. $\omega_{m}=1$). For example, one unit of time is one mechanical period,
and one unit frequency is the mechanical resonance frequency.

Figure 1 fixes the driving strength $\epsilon=1/200$, and the cavity decay
rate $\kappa=1/50$. The evolution time at which the autocorrelation function
is minimized $t_{op}$, the coupling strength $g_{0}$ and the optimized driving
detuning $\Delta$ are listed in Table 3.

$\bm{t_{op}/T}$ | $\bm{g_{0}}$ | $\bm{\Delta}$ | $\bm{t_{op}/T}$ | $\bm{g_{0}}$ | $\bm{\Delta}$ | $\bm{t_{op}/T}$ | $\bm{g_{0}}$ | $\bm{\Delta}$ | $\bm{t_{op}/T}$ | $\bm{g_{0}}$ | $\bm{\Delta}$ | $\bm{t_{op}/T}$ | $\bm{g_{0}}$ | $\bm{\Delta}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
3 | 0.3 | 0.055491 | 3 | 0.35 | -0.01397 | 3 | 0.4 | -0.09137 | 3 | 0.5 | -0.22977 | 3 | 0.6 | -0.41061
4 | 0.3 | -0.00498 | 4 | 0.35 | -0.07028 | 4 | 0.4 | -0.1376 | 4 | 0.5 | -0.22298 | 4 | 0.6 | -0.38487
5 | 0.3 | -0.03953 | 5 | 0.35 | -0.09842 | 5 | 0.4 | -0.0986 | 5 | 0.5 | -0.24291 | 5 | 0.6 | -0.37614
6 | 0.3 | -0.06016 | 6 | 0.35 | -0.06547 | 6 | 0.4 | -0.1269 | 6 | 0.5 | -0.23397 | 6 | 0.6 | -0.38582
7 | 0.3 | -0.06974 | 7 | 0.35 | -0.08281 | 7 | 0.4 | -0.14577 | 7 | 0.5 | -0.24525 | 7 | 0.6 | -0.37771
8 | 0.3 | -0.04821 | 8 | 0.35 | -0.09852 | 8 | 0.4 | -0.13611 | 8 | 0.5 | -0.23984 | 8 | 0.6 | -0.36937
9 | 0.3 | -0.05612 | 9 | 0.35 | -0.10887 | 9 | 0.4 | -0.14063 | 9 | 0.5 | -0.24643 | 9 | 0.6 | -0.36882
10 | 0.3 | -0.06519 | 10 | 0.35 | -0.10828 | 10 | 0.4 | -0.15001 | 10 | 0.5 | -0.24317 | 10 | 0.6 | -0.37111
11 | 0.3 | -0.07255 | 11 | 0.35 | -0.10351 | 11 | 0.4 | -0.15037 | 11 | 0.5 | -0.24718 | 11 | 0.6 | -0.36774
12 | 0.3 | -0.07763 | 12 | 0.35 | -0.10879 | 12 | 0.4 | -0.14792 | 12 | 0.5 | -0.24516 | 12 | 0.6 | -0.36490
13 | 0.3 | -0.07932 | 13 | 0.35 | -0.11361 | 13 | 0.4 | -0.15261 | 13 | 0.5 | -0.24770 | 13 | 0.6 | -0.36573
14 | 0.3 | -0.07634 | 14 | 0.35 | -0.11482 | 14 | 0.4 | -0.15415 | 14 | 0.5 | -0.24642 | 14 | 0.6 | -0.36584
15 | 0.3 | -0.07642 | 15 | 0.35 | -0.11248 | 15 | 0.4 | -0.15218 | 15 | 0.5 | -0.24808 | 15 | 0.6 | -0.36412
Table 3: Optimized parameters for Fig. 1.

Figure 2 fixes the driving strength
$\epsilon=\epsilon_{1}=\epsilon_{2}=1/200$, the cavity decay rate
$\kappa=1/50$, the coupling strength $g_{0}=0.3$ and the evolution time at
which the autocorrelation function is minimized $t_{op}=5T$. The optimized
monochromatic detuning reads $\Delta=-0.03953$. The optimized bichromatic
detunings read $\Delta_{1}=-0.0192947$, $\Delta_{2}=-0.0282437$, and the phase
difference between the two driving fields reads $\psi=2.83667$.

Figure 3 fixes the driving strength $\epsilon_{1}=\epsilon_{2}=1/200$, the
coupling strength $g_{0}=0.3$, the evolution time $t_{op}=5T$. The cavity
decay $\kappa$, the optimized driving detunings $\Delta_{1}$, $\Delta_{2}$,
and the optimized phase difference between the two driving $\psi$ are listed
in Table 4.

$\bm{t_{op}/T}$ | 5 | 5 | 5 | 5 | 5 | 5 | 5 | 5 | 5 | 5
---|---|---|---|---|---|---|---|---|---|---
$\bm{\kappa}$ | $0.01$ | $0.02$ | $0.03$ | $0.04$ | $0.05$ | $0.06$ | $0.07$ | $0.08$ | $0.09$ | $0.1$
$\bm{\Delta_{1}}$ | $-0.03492$ | $-0.01929$ | $-0.01616$ | $-0.0092$ | $-0.00843$ | $-0.00877$ | $-0.00756$ | $-0.00822$ | $-0.00804$ | $-0.00752$
$\bm{\Delta_{2}}$ | $-0.02443$ | $-0.02824$ | $-0.01193$ | $-0.00477$ | $-0.00496$ | $-0.00433$ | $-0.00465$ | $-0.00389$ | $-0.00375$ | $-0.00379$
$\bm{\psi}$ | $3.51878$ | $2.83667$ | $3.27736$ | $3.27906$ | $3.24947$ | $3.27937$ | $3.23194$ | $3.27625$ | $3.27532$ | 3.25781
$\bm{t_{op}/T}$ | 10 | 10 | 10 | 10 | 10 | 10 | 10 | 10 | 10 | 10
$\bm{\kappa}$ | 0.01 | 0.02 | 0.03 | 0.04 | 0.05 | 0.06 | 0.07 | 0.08 | 0.09 | 0.1
$\bm{\Delta_{1}}$ | -0.04765 | -0.11406 | -0.03877 | -0.04014 | -0.04077 | -0.04413 | -0.04394 | 0.95322 | -0.04514 | -0.04499
$\bm{\Delta_{2}}$ | -0.05129 | -0.11768 | -0.04053 | -0.04133 | -0.04378 | -0.04292 | -0.04459 | 0.95523 | -0.04852 | -0.04600
$\bm{\psi}$ | 2.89786 | 2.90915 | 3.032 | 3.06708 | 2.95241 | 3.21805 | 3.10049 | 3.26766 | 2.92939 | 3.0779
Table 4: Optimized parameters for Fig. 3.

Figure 4 fixes the driving strength $\epsilon_{1}=\epsilon_{2}=1/200$, the
cavity decay rate $\kappa=1/50$, the coupling strength $g_{0}=0.3$, and the
evolution time at which the autocorrelation function is minimized $t_{op}=5$.
The optimized driving detunings are $\Delta_{1}=0.0165539$ and
$\Delta_{2}=-0.0364786$, and the optimized phase difference is $\psi=1.4571$.