# Impact of the Central Frequency of Environment on Non-Markovian Dynamics in
Piezoelectric Optomechanical Devices

Quanzhen Ding  Department of Physics, and Center for Quantum Science and
Engineering, Stevens Institute of Technology, Hoboken, NJ 07030, USA  Peng
Zhao  Department of Physics, New York Institute of Technology, Old Westbury,
NY 11568, USA  Department of Engineering and Computing Sciences, New York
Institute of Technology, Old Westbury, NY 11568, USA  Yonghong Ma  School of
Science, Inner Mongolia University of Science and Technology, Baotou 014010,
China<EMAIL_ADDRESS>Yusui Chen  Department of Physics, New York Institute
of Technology, Old Westbury, NY 11568, USA<EMAIL_ADDRESS>
###### Abstract

The piezoelectric optomechanical devices supply a promising experimental
platform to realize the coherent and effective control and measurement on
optical circuits working in Terahertz (THz) frequencies via superconducting
electron devices typically working in Radio (MHz) frequencies. However,
quantum fluctuations are unavoidable when the size of mechanical oscillators
enter into the nanoscale. The consequences of the noisy environment is still
challenging due to the lack of analytical tools. In this paper, a semi-
classical and full-quantum model of piezoelectric optomechanical systems
coupled to a noisy bosonic quantum environment are introduced and solved in
terms of quantum-state diffusion (QSD) trajectories in the non-Markovian
regime. We show that the noisy environment, particularly the central frequency
of environment, can enhance the entanglement generation between optical
cavities and LC circuits in some parameter regimes. Moreover, we observe the
critical points in the coefficient functions, which can lead the different
behaviors in the system. Besides, we also witness the entanglement transfers
between macroscopic objects due to the memory effect of the environment. Our
work can be applied in the fields of electric/ optical switches, and long-
distance distribution in large-scale quantum network.

###### keywords:

Piezo Optomechanics, Non-Markovian dynamics, Master Equation, Entanglement

## Introduction

Optomechanical systems have received substantial interests as a promising
experimental platform to improve the resolution and precision of measurement
to beat the standard quantum limit, to observe macroscopic quantum phenomena,
and to realize sideband cooling and parametric amplification[1, 2, 3, 4, 5, 6,
7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17]. Its applications, ranging from
quantum information to quantum sensing, have been thoroughly studied.
Recently, conventional optomechanical systems consisting of the coupled
mechanical mode and cavity field have been extended to piezoelectric
optomechanical devices in which radio frequency (RF) and cavity field are
parametrically coupled through the coupling with the mechanical oscillator[18,
19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37,
38, 39]. Such hybrid systems supply an experimental platform in which the THz
optical field and GHz acoustic waves are combined. In addition, the mechanical
oscillator can be driven by RF signals and be measured precisely by
interferometers in the optical domain, which raises significant interests for
classical signal processing.

However, it remains a tremendous challenge to suppress the impacts of the
noisy environment[40, 41]. Various schemes have been investigated to address
these challenges, including eliminating the back-action noise by using the
coherent quantum noise cancellation (CQNC) scheme and reducing the shot noise
by lowering the input field power[1, 10, 42, 43, 44]. Thermal noise is another
source of noise unless the temperature is very low. But when the system is
cooled down to few milli-Kelvin, quantum fluctuations exceed the thermal noise
and become the major noise. So a fully understanding of the dynamics of the
piezoelectric optomechanical (POM) system in the present of noisy environment
is crucial. This paper addresses the evolution of POM systems starting from a
complete microscopic Hamiltonian. Our proposed model consists of a mechanical
oscillator coupled to an optical cavity and an LC circuit simultaneously, as
shown in Fig. 1.

Figure 1: A schematic of our model, consisting of a piezo-mechanical resonator
coupled with an optical cavity and a LC circuit simultaneously.

Traditional methods to deal with our proposed model are limited in Markovian
regime[41], assuming that the environment is too large to back-act on the
system and the relaxation time of the system is much shorter comparing to the
operational time. However, many future applications require high-frequency
operation, for example the quantum network and quantum computing. Without loss
of generalization, we need to initiate our discussion in non-Markovian regime.
Moreover, the non-Markovian environment can parametrically tune the dynamics
of the system and enhance the entanglement in some particular regime, unlike
the Markovian environment in which only an exponential decay can be observed.

In this paper, we investigate the POM model that a mechanical oscillator,
coupled to an optical cavity and an LC circuit simultaneously and embedded in
a non-Markovian environment. The aim of the paper is twofold. Firstly, we
derive a non-Markovian stochastic Schrödinger equation (SSE) to characterize
the dynamics of the system by using the fully quantized Hamiltonian[45, 46].
Each numerical solution of the SSE is so-called a trajectory and the reduced
density operator of the system can be reproduced by taking average of all
trajectories. Moreover, the corresponding master equation can be obtained
formally[47, 48, 42]. Secondly, we investigate how the entanglement is
influenced by various non-Markovian environment[43, 44]. By adjusting the
frequencies of both optical and mechanical mode, coupling strength and central
frequency of the environment, we carefully study the entanglement transfer
inside the POM system. We discuss the physical mechanism in the model and the
experimental feasibility of our proposed model.

## Background

In this section, we briefly review the relevant background of non-Markovian
dynamics, particularly using the quantum-state diffusion (QSD) approach and
the corresponding master equation (MEQ) approach[47, 42].

### General non-Markovian dynamics of open quantum systems

A general open quantum system can be studied in the system-environment
framework and the formal Hamiltonian is expressed as

$\displaystyle H_{tot}=H_{sys}+H_{int}+H_{env}.$

The quantized bosonic environment Hamiltonian, and the interaction Hamiltonian
can be expressed as

$\displaystyle H_{env}$
$\displaystyle=\sum_{k}\omega_{k}d_{k}^{\dagger}d_{k},$ $\displaystyle
H_{int}$
$\displaystyle=\sum_{k}(g_{k}Ld_{k}^{\dagger}+g_{k}^{*}L^{\dagger}d_{k}).$

where $\omega_{k}$ is the eigen frequency of the $k$th mode in the
environment. The coupling between system and environment is characterized by
the Lindblad operator $L$. Under some conditions, e.g. the environment is too
large to allow the existence of back flows, and the system-environment
coupling is relatively weak, the state of the system can be solely depended on
the state at the last time instant. Due to the so-called Born-Markov
approximation, the dynamics of the system can be described by using Lindblad
master equations, Redfield master equations, quantum jump, Heisenberg approach
etc. However, generalizing a non-Markovian quantum dynamics is difficult.

The QSD approach is a method to describe the non-Markovian dynamics of the
system by tuning it into a set of continuous time stochastic processes. Each
of the stochastic evolution of the state is called a trajectory. And the
reduced density operator of the system can be regenerated numerically as the
ensemble average of all trajectories. In the interaction picture, the QSD
equation is given as

$\displaystyle\partial_{t}\psi_{t}=(-iH_{sys}+Lz_{t}^{*}-L^{\dagger}\int_{0}^{t}ds\
\alpha(t,s)\frac{\delta}{\delta z_{s}^{*}})\psi_{t},$ (1)

where $\psi_{t}=\langle z|\Psi_{tot}\rangle$ is the continuous time stochastic
trajectory of the state of the system.
$z_{t}^{*}=-i\sum_{k}g_{k}z_{k}^{*}e^{-i\omega_{k}t}$ is a complex Gaussian
stochastic process. $|z\rangle=\otimes_{k}|z_{k}\rangle$ is the collective
basis of all modes in the environment, where the $|z_{k}\rangle$ is the
Bargamann coherent state of the $k$th mode, defined as
$d_{k}|z_{k}\rangle=z_{k}|z_{k}\rangle$. With the completeness identity
$I=\int\frac{d^{2}z}{\pi}e^{-|z|^{2}}|z\rangle\langle z|$, the reduced density
operator can formally expressed as

$\displaystyle\rho_{sys}=\mathcal{M}[|\psi_{t}\rangle\langle\psi_{t}|]=\int_{0}^{t}\frac{d^{2}z}{\pi}e^{-|z|^{2}}|\psi_{t}\rangle\langle\psi_{t}|,$
(2)

where the symbol $\mathcal{M}[\cdot]$ represents the ensemble average.
Moreover, Gaussian stochastic processes satisfy that $\mathcal{M}[z_{t}]=0$,
and are characterized by the correlation functions:
$\mathcal{M}[z_{t}z_{s}]=0$, and $\mathcal{M}[z_{t}^{*}z_{s}]=\alpha(t,s)$. In
particular, $\alpha(t,s)$ is defined as

$\displaystyle\alpha(t-s)=\int d\omega\ J(\omega)e^{-i\omega(t-s)},$ (3)

and depends solely on the spectral function
$J(\omega)=\sum_{k}|g_{k}|^{2}\delta(\omega-\omega_{k})$, which describes how
the system-environment coupling strength is frequency depended. For instance,
the correlation function $\alpha(t,s)$ will take a form of Dirac delta
function when the distribution of $J(\omega)$ is flat. By considering
different coupling spectral functions and the corresponding correlation
functions, the QSD approach can be applied to study a variety of types of
environments.

After replacing the functional derivative $\frac{\delta\psi_{t}}{\delta
z_{s}^{*}}$ by the product of an operator $O(t,s,z^{*})$ in the system’s
Hilbert space and the stochastic wave function, $\frac{\delta\psi_{t}}{\delta
z_{s}^{*}}\ =O(t,s,z^{*})\psi_{t}$, the time-local QSD equation (1) can be
transformed to

$\displaystyle\partial_{t}\psi_{t}=(-iH_{sys}+Lz_{t}^{*}-L^{\dagger}\int_{0}^{t}ds\
\alpha(t,s)O(t,s))\psi_{t},$ (4)

where the $O$ operator can be determined by an evolution equation that

$\displaystyle\partial_{t}O$
$\displaystyle=\left[-iH_{sys}+Lz_{t}^{*}-L^{\dagger}\bar{O},O\right]-L^{\dagger}\frac{\delta\bar{O}}{\delta
z_{s}^{*}},$ (5)

together with its initial condition that $O(t,s=t)=L$, where
$\bar{O}=\int_{0}^{t}ds\alpha(t,s)O(t,s)$.

With the QSD equation, the corresponding non-Markovian master equations can be
derived as

$\displaystyle\partial_{t}\rho_{sys}=-i[H_{sys},\>\rho_{sys}]+[L,\>\mathcal{M}(|\psi_{t}\rangle\langle\psi_{t}|\bar{O}^{\dagger})]-[L^{\dagger},\>\mathcal{M}(\bar{O}|\psi_{t}\rangle\langle\psi_{t}|)].$
(6)

### Macroscopic entanglement

To study the entanglement dynamics between components in the system, we employ
the negativity to measure the scale of entanglement,

$\displaystyle\mathcal{N}(\rho)=\frac{||\rho^{\Gamma_{A}}||_{1}-1}{2}=\sum_{i}\frac{|\lambda_{i}|-\lambda_{i}}{2},$

where $\rho^{\Gamma_{A}}$ is the partial transpose of $\rho$ with respect to
the subsystem $A$. $||\cdot||_{1}$ represents the trace norm. And
$\\{\lambda_{i}\\}$ are all negative eigenvalues of $\rho^{\Gamma_{A}}$. The
quantity $\mathcal{N}(\rho)$ is monotonic to the upper boundary of the
distillable entanglement.

## The Model

The model we proposed is illustrated in Fig. 1. An optomechanical system,
consisting of an optical cavity and a movable mirror, is coupled to a high-
frequency piezo-mechanical resonator. As a result, the mechanical oscillator
is coupled to the optical cavity and the capacitor of the superconducting LC
resonator simultaneously, and achieves the coupling between the optical field
and the microwave field. Moreover, a quantum environment coupled to the
mechanical oscillator may be introduced by a variety of imperfections, e.g.,
due to diffraction, external fields, lattice vibrations etc. Thus, the system
Hamiltonian can be expressed as (setting $\hbar=1$)

$\displaystyle H_{sys}$
$\displaystyle=-\omega_{o}a^{\dagger}a+\omega_{m}b^{\dagger}b+\omega_{e}c^{\dagger}c-g_{om}a^{\dagger}a\hat{x}-g_{me}\hat{x}\left(c+c^{\dagger}\right)$
(7)

where $a(a^{\dagger})$, $b(b^{\dagger})$, and $c(c^{\dagger})$ are the
annihilation (creation) operators of the optical field, the mechanical mode,
and the L-C resonator respectively, with the eigen frequency of $\omega_{o}$,
$\omega_{m}$ and $\omega_{e}$. The displacement operator
$\hat{x}=x_{ZPF}(b+b^{\dagger})$, where $x_{ZPF}=\sqrt{1/2m\omega_{m}}$ is the
zero point frequency of the mechanical mode. The optomechanical interaction
Hamiltonian $g_{om}\hat{x}a^{\dagger}a$ is nonlinear for the displacement of
the movable mirror $\hat{x}$, which only is valid when the mechanical
oscillation frequency is relatively smaller than the free spectral range of
the cavity. The piezoelectric interaction Hamiltonian
$g_{me}\hat{x}(c+c^{\dagger})$ describes a linear coupling, where the coupling
strength depends on the coupling between the strain and the electric field.

It is difficult to analytically solve the dynamics of the system, particularly
when the optomechanical coupling in nonlinear. However, in the semi-classical
limit, the optical mode operator $a$ can be linearized and expended as the sum
of the average amplitude $\alpha$ and the fluctuation $\delta a$, that
$a=\alpha+\delta a$. As a result,the nonlinear interaction
$\hat{x}a^{\dagger}a$ can be rewritten as $\hat{x}\left(\alpha^{*}+\delta
a^{\dagger}\right)\left(\alpha+\delta a\right)$. Assuming that the fluctuation
$\delta a$ is much smaller than the average amplitude $\alpha$, the final
linearization of the interaction Hamiltonian can be expressed as

$\displaystyle H_{om}=g_{om}\hat{x}(a^{\dagger}+a).$

### System-Environment couplings

To consider a dissipative environment, the interaction Hamiltonian is given

$\displaystyle
H_{int}=b\sum_{k}g_{k}d_{k}^{\dagger}+b^{\dagger}\sum_{k}g_{k}^{*}d_{k}.$

There are some typical spectral functions such as Lorentzian, Ohmic, sub-
Ohmic, and super-Ohmic for the structure of the environment. In order to
observe how the dynamics of the system gradually transits from the Markovian
regime to the non-Markovian regime, we consider the Lorentzian spectral
density as

$\displaystyle
J(\omega)=\frac{1}{2\pi}\frac{\Gamma\gamma^{2}}{(\omega-\Omega)^{2}+\gamma^{2}},$

where $\Omega$ is the central frequency of the environment. The parameter
$\gamma$ defines the spectral width of the coupling and $1/\gamma$ indicates
the memory time of the environment. $\Gamma$ is the global coupling strength.
By substituting $J(\omega)$ into Eq. (3), the correlation function
$\alpha(t,s)$ can be determined analytically, in the form of the
Ornstein–Uhlenbeck $(O-U)$ correlation function,

$\displaystyle\alpha(t,s)=\frac{\Gamma\gamma}{2}e^{-\gamma|t-s|-i\Omega(t-s)}.$

The analytical solution of the dissipative dynamics of the system depends on
the system Hamiltonian and the spectral density of the environment
simultaneously. In this section, we introduce two classes of coupling in the
system, and in each case, we study the non-Markovian effects influenced by the
central frequency of the environment.

### Case I: Weak coupling

For the piezoelectric interaction, if the coupling strength $g_{me}$ is weak,
comparing to the eigen frequencies of the mechanical mode and the LC resonator
$\omega_{m}$ and $\omega_{e}$, then it is reasonable to use the rotating-wave
approximation to simplify

$\displaystyle H_{me}=g_{me}(b+b^{\dagger})(c+c^{\dagger})\approx
g_{me}(bc^{\dagger}+b^{\dagger}c).$

And consequently, the interaction of the optical field and the mechanical
oscillator can approximate to

$\displaystyle
g_{om}\left(a+a^{\dagger}\right)\left(b+b^{\dagger}\right)\approx
g_{om}\left(ab^{\dagger}+a^{\dagger}b\right).$

From Eq. (5), the ansatz of $O$ operator must take the form of

$\displaystyle
O\left(t,s\right)=f_{1}\left(t,s\right)b+f_{2}\left(t,s\right)a+f_{3}\left(t,s\right)c,$

wither three to-be determined coefficient functions $f_{1}(t,s)$,
$f_{2}(t,s)$, and $f_{3}(t,s)$. Substituting the ansatz into Eq. (5), the
three coefficient functions can be determined by a group of partial
differential equations

$\displaystyle\partial_{t}f_{1}$
$\displaystyle=i\omega_{m}f_{1}-ig_{om}f_{2}-ig_{me}f_{3}+F_{1}f_{1},$
$\displaystyle\partial_{t}f_{2}$
$\displaystyle=-i\Delta_{o}f_{2}-ig_{om}f_{1}+F_{2}f_{1},$
$\displaystyle\partial_{t}f_{3}$
$\displaystyle=i\omega_{e}f_{1}-ig_{me}f_{1}+F_{3}f_{1},$

with the intial conditions

$\displaystyle f_{1}\left(t,t\right)$ $\displaystyle=1,$ $\displaystyle
f_{2}\left(t,t\right)$ $\displaystyle=f_{3}\left(t,t\right)=0.$

Due to the exponential $O-U$ style correlation function
$\alpha\left(t,s\right)=\frac{\Gamma\gamma}{2}e^{-\gamma\left|t-s\right|-i\Omega\left(t-s\right)}$,
we can continue to analytically derive the evolution function of the operator
$\bar{O}=F_{1}\left(t\right)b+F_{2}\left(t\right)a+F_{3}\left(t\right)c$,
where $F_{j}(t)=\int_{0}^{t}ds\alpha(t,s)f_{j}(t,s)$. Hence the $\bar{O}$
operator can be numerically found by the following group of expressions,

$\displaystyle\partial_{t}F_{1}$
$\displaystyle=\frac{\Gamma\gamma}{2}-\left(\gamma+i\Omega-i\omega_{m}\right)F_{1}-ig_{om}F_{2}-ig_{me}F_{3}+F_{1}^{2},$
$\displaystyle\partial_{t}F_{2}$
$\displaystyle=-\left(\gamma+i\Omega+i\Delta_{o}\right)F_{2}-ig_{om}F_{1}+F_{1}F_{2},$
$\displaystyle\partial_{t}F_{3}$
$\displaystyle=-\left(\gamma+i\Omega-i\omega_{e}\right)F_{3}-ig_{me}F_{1}+F_{1}F_{3}.$

Figure 2: Evolution of the real ($a,b$) and imaginary ($c,d$) value of the
complex coefficient function $F_{1}(t)$ in the $\bar{O}$ operator for Case I.
In (a) and (c), $\gamma=0.5$. In (b) and (d), $\gamma=5.0$. Other parameters
are chosen as: $\omega_{o}=\omega_{e}=\omega_{m}$,
$g_{om}=g_{me}=0.1\omega_{m}$, $\Gamma=1$.

In Fig. 2, the dynamics of the coefficient $F_{1}(t)$ in the $\bar{O}$
operator is illustrated when the central frequency of environment $\Omega$
equals $0,1$ and $1.6$. We only plot $F_{1}$ as it is one order of magnitude
larger than $F_{2}$ and $F_{3}$, and takes the dominant role in the $\bar{O}$
operator. In Fig. 2 (a, c), $F_{1}$’s real and imaginary values oscillate in
the short time interval and gradually approach to a stable value in the long-
time limit for $\Omega=0$ (blue dashed line), $1.8$ (green dash-dotted line).
When $\Omega\approx 1$ (red solid line), singularities emerge in the real and
imaginary parts of $F_{1}(t)$, which also show a periodic behavior. In Fig.
2(b, d), when the decay factor $\gamma=5$, $F_{1}$ approaches to the stable
value $\frac{\Gamma}{2}$ quickly. In addition, since the $\bar{O}$ operator is
noise free, thus
$\mathcal{M}[|\psi_{t}\rangle\langle\psi_{t}|\bar{O}^{\dagger}]=\rho_{sys}\bar{O}^{\dagger}$
and consequently the formal master equation can be explicitly written as

$\displaystyle\partial_{t}\rho_{sys}$
$\displaystyle=-i\left[H_{sys},\rho_{sys}\right]+\left[L,\rho_{sys}\bar{O}^{\dagger}\right]-\left[L^{\dagger},\>\bar{O}\rho_{sys}\right]$
$\displaystyle=-i\left[-\omega_{o}a^{\dagger}a+\omega_{m}b^{\dagger}b+\omega_{e}c^{\dagger}c-g_{om}\left(ab^{\dagger}+a^{\dagger}b\right)-g_{me}\left(bc^{\dagger}+b^{\dagger}c\right),\rho_{sys}\right]$
$\displaystyle\quad+\left[b,\rho\left(F_{1}^{*}b^{\dagger}+F_{2}^{*}a^{\dagger}+F_{3}^{*}c^{\dagger}\right)\right]+\left[\left(F_{1}b+F_{2}a+F_{3}c\right)\rho_{sys},b^{\dagger}\right].$

### Case II: Strong coupling

In this section, we will consider a more general coupling that the mechanical
oscillator is strongly coupled to the optical cavity and the LC circuit,
therefore the system Hamiltonian is

$\displaystyle H_{sys}$
$\displaystyle=-\omega_{o}a^{\dagger}a+\omega_{m}b^{\dagger}b+\omega_{e}c^{\dagger}c-g_{om}\left(a^{\dagger}+a\right)\hat{x}-g_{me}\hat{x}\left(c+c^{\dagger}\right).$

The anstz of $O$ operator satisfying Eq. (5) may be found

$\displaystyle O(t,s)$
$\displaystyle=O_{0}(t,s)+\int_{0}^{t}ds^{\prime}O_{1}(t,s,s^{\prime})z_{s^{\prime}},$
$\displaystyle O_{0}\left(t,s\right)$
$\displaystyle=f_{1}\left(t,s\right)b+f_{2}\left(t,s\right)b^{\dagger}+f_{3}\left(t,s\right)a+f_{4}\left(t,s\right)a^{\dagger}+f_{5}\left(t,s\right)c+f_{6}\left(t,s\right)c^{\dagger},$
$\displaystyle O_{1}\left(t,s,s^{\prime}\right)$
$\displaystyle=f_{7}\left(t,s,s^{\prime}\right)I,$

where $O_{0}(t,s)$ and $O_{1}(t,s,s^{\prime})$ represent the noise-free and
first order noise terms in the $O$ operator respectively. Similarly, by
substituting the anstz of the $O$ operator into Eq. 5, all coefficient
functions can be numerically found from the following expressions

$\displaystyle\partial_{t}f_{1}$
$\displaystyle=i\omega_{m}f_{1}-ig_{om}f_{3}+ig_{om}f_{4}-ig_{me}f_{5}+ig_{me}f_{6}+F_{1}f_{1},$
$\displaystyle\partial_{t}f_{2}$
$\displaystyle=-i\omega_{m}f_{2}-ig_{om}f_{3}+ig_{om}f_{4}-ig_{me}f_{5}+ig_{me}f_{6}+2F_{2}f_{1}-F_{1}f_{2}+F_{4}f_{3}-f_{3}F_{4}+F_{6}f_{5}-F_{5}f_{6}-F_{7}\left(t,s\right),$
$\displaystyle\partial_{t}f_{3}$
$\displaystyle=-i\Delta_{o}f_{3}-ig_{om}f_{1}+ig_{om}f_{2}+F_{3}f_{1},$
$\displaystyle\partial_{t}f_{4}$
$\displaystyle=i\Delta_{o}f_{4}-ig_{om}f_{1}+ig_{om}f_{2}+F_{4}f_{1},$
$\displaystyle\partial_{t}f_{5}$
$\displaystyle=i\omega_{e}f_{5}-ig_{me}f_{1}+ig_{me}f_{2}+F_{5}f_{1},$
$\displaystyle\partial_{t}f_{6}$
$\displaystyle=-i\omega_{e}f_{6}-ig_{me}f_{1}+ig_{me}f_{2}+F_{6}f_{1},$
$\displaystyle\partial_{t}f_{7}$
$\displaystyle=F_{7}\left(t,s^{\prime}\right)f_{1},$

with the initial and boundary conditions

$\displaystyle f_{1}\left(t,t\right)$ $\displaystyle=1,$ $\displaystyle
f_{2}\left(t,t\right)$
$\displaystyle=f_{3}\left(t,t\right)=f_{4}\left(t,t\right)=f_{5}\left(t,t\right)=f_{6}\left(t,t\right)=0,$
$\displaystyle f_{7}\left(t,s,s^{\prime}=t\right)$ $\displaystyle=f_{2}(t,s),$
$\displaystyle f_{7}\left(t,s=t,s\right)$ $\displaystyle=0.$

Consequently, the $\bar{O}$ may be determined by

$\displaystyle\partial_{t}F_{1}$
$\displaystyle=\frac{\Gamma\gamma}{2}-\left(\gamma+i\Omega-i\omega_{m}\right)F_{1}-ig_{om}F_{3}+ig_{om}F_{4}-ig_{me}F_{5}+ig_{me}F_{6}+F_{1}^{2},$
$\displaystyle\partial_{t}F_{2}$
$\displaystyle=-\left(\gamma+i\Omega+i\omega_{m}\right)F_{2}-ig_{om}F_{3}+ig_{om}F_{4}-ig_{me}F_{5}+ig_{me}F_{6}+F_{1}F_{2}-\tilde{F}_{7},$
$\displaystyle\partial_{t}F_{3}$
$\displaystyle=-\left(\gamma+i\Omega+i\Delta_{o}\right)F_{3}-ig_{om}F_{1}+ig_{om}F_{2}+F_{1}F_{3},$
$\displaystyle\partial_{t}F_{4}$
$\displaystyle=-\left(\gamma+i\Omega-i\Delta_{o}\right)F_{4}-ig_{om}F_{1}+ig_{om}F_{2}+F_{1}F_{4},$
$\displaystyle\partial_{t}F_{5}$
$\displaystyle=-\left(\gamma+i\Omega-i\omega_{e}\right)F_{5}-ig_{me}F_{1}+ig_{me}F_{2}+F_{1}F_{5},$
$\displaystyle\partial_{t}F_{6}$
$\displaystyle=-\left(\gamma+i\Omega+i\omega_{e}\right)F_{6}-ig_{me}F_{1}+ig_{me}F_{2}+F_{1}F_{6},$

and

$\displaystyle\partial_{t}\tilde{F}_{7}=\frac{\Gamma\gamma}{2}F_{2}-2\left(\gamma+i\Omega\right)\tilde{F}_{7}+F_{1}\tilde{F}_{7},$

where $\tilde{F}_{7}(t)=\int_{0}^{t}ds\alpha(t,s)F_{7}(t,s)$. Numerical
estimations for $F_{1}(t)$ and $F_{2}(t)$ are illustrated in Fig. 3.
Singularities are observed in Fig. 3 (a,c), the non-Markovian evolution of
both $F_{1}$ and $F_{2}$, when $\Omega$ is close to $1.18$. Although the exact
$O$ operator is obtained, the master equation is still hard to derive, as the
$O$ operator consists of the first order noise. According to a systematic
method that derives master equations[47], the ensemble average
$\mathcal{M}\left[P_{t}\bar{O}^{\dagger}\right]$ can be derived analytically
as,

$\displaystyle
R\left(t\right)=\mathcal{M}\left[P_{t}\bar{O}^{\dagger}\right]=\rho_{t}\bar{O}_{0}^{\dagger}+\iint_{0}^{t}ds_{1}ds_{2}\alpha\left(s_{1},s_{2}\right)O_{0}\left(t,s_{2}\right)\rho_{t}\bar{O}_{1}^{\dagger}\left(t,s_{1}\right),$

in the time-convolutionless form. As a result, the exact master equation can
be given by

$\displaystyle\partial_{t}\rho$
$\displaystyle=-i\left[H_{sys},\rho\right]+\left[L,R\left(t\right)\right]+\left[R^{\dagger}\left(t\right),L^{\dagger}\right].$

Figure 3: Evolution of the real and imaginary value of the complex coefficient
function $F_{1}(t)$. In (a) and (c), $\gamma=0.5$, non-Markovian regime. In
(b) and (d), $\gamma=5.0$, Markovian regime. Parameters are chosen as:
$\omega_{o}=\omega_{e}=\omega_{m}$, $g_{om}=g_{me}=0.1\omega_{m}$, $\Gamma=1$.

In Fig. 3, the dynamics of the real and imaginary parts of $F_{1}$ are similar
to that in the case I, when $\gamma=0.5$. They all approach to fixed values in
the long time limit, for $\Omega=0$ (Blue dashed line), $1.8$ (Green dash-
dotted line). But singularities in $F_{1}$ real part are observed when
$\Omega\approx 1$ (Red solid line). But the imaginary part of $F_{1}$ stays at
zero. In (b) and (d), When $\gamma=5$, the dynamics of $F_{1}$ is Markovian
and no correlated with the central frequency $\Omega$.

## Result

In this section, we look in the two classes of couplings that we have
considered, and illustrate the entanglement dynamics between the electric-,
optical- and mechanical- component in the system. However, for different
parameters, we observe that the non-Markovian dynamics transfers from one type
to another. In order to do so, we plot the entanglement dynamics versus the
central frequency of environment $\Omega$ for different types of coupling. We
provide a summary of all the results obtained in this section and connect the
results with the singularities observed in the coefficient functions of the
$O$ operator.

### Case I

Figure 4: Dynamics of entanglement between Om and me components for different
decay factor $\gamma$. Other parameters are chosen as
$\omega_{o}=\omega_{e}=\omega_{m}$, $g_{om}=g_{me}=0.1\omega_{m}$, $\Gamma=1$.

First of all, for this exciton hopping model, we choose the initial state that
there are two excitons existed in the optical cavity and the LC resonator,
$\psi_{sys}(t=0)=|1\rangle_{O}\otimes|0\rangle_{M}\otimes|1\rangle_{E}$, to
mimic the process that a photon and a phonon couple to each other via the
interactions with the mechanical oscillator. Then, we investigate how does the
memory time modify the entanglement dynamics. From the numerical analysis of
the coefficient functions in the $\bar{O}$ operator (Fig. 2), the memory
factor $1/\gamma$ is dominant when it is close to the Markovian limit. We plot
the dynamics of entanglement between optical cavity and mechanical oscillator
($\mathcal{N}_{om}$) in Fig. 4 (a), and the entanglement between mechanical
oscillator and LC circuit ($\mathcal{N}_{me}$) is displayed in Fig. 4 (b)
respectively. For the dissipative coupling described in case I, both
$\mathcal{N}_{om}$ and $\mathcal{N}_{me}$ decay exponentially to zero when
$\gamma=5$ and in the Markov limit. When $\gamma=0.1$ and $0.2$, the
entanglement dynamics shows the periodic behavior. And we notice that the
maximal $\mathcal{N}_{om}$ happens when the $\mathcal{N}_{me}$ achieves its
minimum value. This phenomena clearly display how the entanglement transfers
between different components in the system and offers a potential application
in realizing quantum non-demolition measurements.

Figure 5: Dynamics of entanglement between OM for different central
frequencies of the environment $\Omega$. Other parameters are chosen as
$\omega_{o}=\omega_{e}=\omega_{m}$, $g_{om}=g_{me}=0.1\omega_{m}$, $\Gamma=1$.
Figure 6: Dynamics of entanglement between ME for different central
frequencies of the environment $\Omega$. Other parameters are chosen as
$\omega_{o}=\omega_{e}=\omega_{m}$, $g_{om}=g_{me}=0.1\omega_{m}$, $\Gamma=1$.

In Fig. 5, we focus on the impact of the central frequency of the environment
on the entanglement dynamics. We choose $\gamma=1$, and adjust $\Omega$ from
$0$ to $5$. In the short time interval, the generation of entanglement
$\mathcal{N}_{om}$ is caused by the linear coupling, so it is independent of
the value of $\Omega$. After that, two diametrically different classes of non-
Markovian dynamics are observed. The boundary of the two classes of non-
Markovian dynamics locates at $\Omega\approx 1$. When $0\leq\Omega<1$,
$\mathcal{N}_{om}$ decay non-monotonically, and after sufficient time, it
decays to zero. When $\Omega>1$, the behavior of entanglement is periodic and
approaching to a non-zero value. From the above discussion, the two types of
non-Markovian dynamics can be categorized by the scale of the long-time
entanglement, due to the different range of central frequency of the
environment $\Omega$. In Fig. 6, the entanglement $\mathcal{N}_{me}$ achieve
the maximal value and decay to zero quickly after for different $\Omega$
values. However, we noticed that the $\omega$ value can modify the maximal
value of the negativity $\mathcal{N}_{me}$.

### Case II

Figure 7: Dynamics of entanglement between (a) OM and (b) ME, for different
$\gamma$ values. Other parameters are chosen as
$\omega_{o}=\omega_{e}=\omega_{m}$, $g_{om}=g_{me}=0.1\omega_{m}$, $\Gamma=1$.

In this subsection, we consider the strong coupling model. Since the number of
exciton is no longer conserved as in case I, we choose the initial state of
the system as $|\psi_{sys}(t=0)=|000\rangle$, as we will focus on the
generation of entanglement happened between the three components in the
system. Interestingly, we find that the dynamics of $\mathcal{N}_{om}$ and
$\mathcal{N}_{me}$ are very different, which can be seen in Fig. 7. For
$\mathcal{N}_{om}$ in Fig. 7 (a), the maximal generated entanglement is
positively correlated to the memory factor $1/\gamma$. While
$\mathcal{N}_{me}$, in Fig. 7 (b), entanglement sudden death (ESD) is observed
for all $\gamma$ values. Moreover, the time that the ESD happens is shorter
when the environment is further from the Markovian limit. At $\omega_{m}t>60$,
$\mathcal{N}_{om}$ and $\mathcal{N}_{me}$ decay to zero and indicate that the
three components of the system are disentangled with each other at that
instant. This feature can be used to synchronize the remote quantum network by
tuning the coupling strength with the mechanical oscillator.

Figure 8: Dynamics of entanglement between OM , for different $\omega$ values.
Other parameters are chosen as $\omega_{o}=\omega_{e}=\omega_{m}$,
$g_{om}=g_{me}=0.1\omega_{m}$, $\Gamma=1$. Figure 9: Dynamics of entanglement
between ME, for different $\Omega$ values. The dashed line marks the time that
entanglement sudden death happens. It shows that the time to observe ESD has a
linear relationship with the central frequency of the environment $\Omega$.
Other parameters are chosen as $\omega_{o}=\omega_{e}=\omega_{m}$,
$g_{om}=g_{me}=0.1\omega_{m}$, $\Gamma=1$.

In Fig. 8, impact of the central frequency of the environment on the
entanglement dynamics is that $\Omega$ can restrain the maximal generated
entanglement between OM. But the long-time limit entanglement approaches to
zero. In Fig. 9, the generation of entanglement $\mathcal{N}_{me}$ achieves
its maximal value and then decay to zero in the long-time limit. But we notice
that the non-Markovian behavior can be categorized in two types again: (1),
$\Omega<1$, $\mathcal{N}_{me}$ oscillates around a fixed value in the short-
time regime and quickly decay to zero right after; (2), $\Omega>1$,
$\mathcal{N}_{me}$ decays to zero linearly. In addition, we study the time
when ESD happens ($t_{ESD}$) and find out that the time and the central
frequency of the environment $\Omega$ have a quantitative relationship,

$\displaystyle t_{ESD}\approx\beta\Omega,$

which supplies a quantitatively method to estimate ESD and design the control
and measurement schemes, where $\beta$ is the proportional factor and its
value depends on the parameters chosen in the numerical simulations.

## Conclusions

In this paper, we study the entanglement dynamics and transfer between OM and
ME components in the piezoelectric optomechanical system. The non-Markovian
effects from the environment are discussed in weak coupling and strong
coupling respectively. In particular, the effects of the central frequency of
the environment on the entanglement dynamics are discussed. By using the
quantum-state diffusion approach and the corresponding master equation, We
analyze the coefficient functions evolution and discover that singularities
emerge when the central frequency of the environment is around a special
value. As a result, the non-Markovian dynamics of the system can be
categorized in two types: (1) the entanglement gradually decays to zero; (2)
the entanglement shows periodic behavior and approaches to a fixed value. This
categorization can help prepare the quantum system and its coupled environment
toward steady states and realize quantum teleportation in the quantum network.
Moreover, the entanglement generation is studied. An quantitative conclusion
is that the time to observe ESD $t_{ESD}$ is proportional to $\Omega$, which
offers a robust and simple method to estimate the entanglement in the quantum
network. Our work can be extended to study the non-equilibrium dynamics in a
full quantum framework when the environment temperature is close to zero
Kelvin.

## References

  * [1] Aspelmeyer, M., Kippenberg, T. J. & Marquardt, F.  Cavity optomechanics.  _Rev. Mod. Phys._ 86, 1391–1452, DOI: 10.1103/RevModPhys.86.1391 (2014).
  * [2] Brooks, D. W. C. _et al._ Non-classical light generated by quantum-noise-driven cavity optomechanics.  _Nature_ 488, 476–480, DOI: 10.1038/nature11325 (2012).
  * [3] Mancini, S., Vitali, D. & Tombesi, P.  Optomechanical cooling of a macroscopic oscillator by homodyne feedback.  _Phys. Rev. Lett._ 80, 688–691, DOI: 10.1103/PhysRevLett.80.688 (1998).
  * [4] Nimmrichter, S., Hornberger, K. & Hammerer, K.  Optomechanical sensing of spontaneous wave-function collapse.  _Phys. Rev. Lett._ 113, 020405, DOI: 10.1103/PhysRevLett.113.020405 (2014).
  * [5] Vitali, D. _et al._ Optomechanical entanglement between a movable mirror and a cavity field.  _Phys. Rev. Lett._ 98, 030405, DOI: 10.1103/PhysRevLett.98.030405 (2007).
  * [6] Otterstrom, N. T., Behunin, R. O., Kittlaus, E. A. & Rakich, P. T.  Optomechanical cooling in a continuous system.  _Phys. Rev. X_ 8, 041034, DOI: 10.1103/PhysRevX.8.041034 (2018).
  * [7] Kippenberg, T. J. & Vahala, K. J.  Cavity optomechanics: Back-action at the mesoscale.  _Science_ 321, 1172–1176, DOI: 10.1126/science.1156032 (2008).  https://science.sciencemag.org/content/321/5893/1172.full.pdf.
  * [8] Weis, S. _et al._ Optomechanically induced transparency.  _Science_ 330, 1520–1523, DOI: 10.1126/science.1195596 (2010).  https://science.sciencemag.org/content/330/6010/1520.full.pdf.
  * [9] Ma, Y. _et al._ Proposal for gravitational-wave detection beyond the standard quantum limit through epr entanglement.  _Nature Physics_ 13, 776–780, DOI: 10.1038/nphys4118 (2017).
  * [10] Chen, Y.  Macroscopic quantum mechanics: theory and experimental concepts of optomechanics.  _Journal of Physics B: Atomic, Molecular and Optical Physics_ 46, 104001, DOI: 10.1088/0953-4075/46/10/104001 (2013).
  * [11] Baldacci, L. _et al._ Thermal noise and optomechanical features in the emission of a membrane-coupled compound cavity laser diode.  _Scientific Reports_ 6, 31489, DOI: 10.1038/srep31489 (2016).
  * [12] Villafañe, V. _et al._ Optoelectronic forces with quantum wells for cavity optomechanics in gaas/alas semiconductor microcavities.  _Phys. Rev. B_ 97, 195306, DOI: 10.1103/PhysRevB.97.195306 (2018).
  * [13] Boales, J. A., Mateen, F. & Mohanty, P.  Micromechanical resonator driven by radiation pressure force.  _Scientific Reports_ 7, 16056, DOI: 10.1038/s41598-017-16063-4 (2017).
  * [14] Fang, K., Matheny, M. H., Luan, X. & Painter, O.  Optical transduction and routing of microwave phonons in cavity-optomechanical circuits.  _Nature Photonics_ 10, 489–496, DOI: 10.1038/nphoton.2016.107 (2016).
  * [15] Renninger, W. H., Kharel, P., Behunin, R. O. & Rakich, P. T.  Bulk crystalline optomechanics.  _Nature Physics_ 14, 601–607, DOI: 10.1038/s41567-018-0090-3 (2018).
  * [16] Santos, J. T., Li, J., Ilves, J., Ockeloen-Korppi, C. F. & Sillanpää, M.  Optomechanical measurement of a millimeter-sized mechanical oscillator approaching the quantum ground state.  _New Journal of Physics_ 19, 103014, DOI: 10.1088/1367-2630/aa83a5 (2017).
  * [17] Cripe, J. _et al._ Radiation-pressure-mediated control of an optomechanical cavity.  _Phys. Rev. A_ 97, 013827, DOI: 10.1103/PhysRevA.97.013827 (2018).
  * [18] Han, X., Zou, C.-L. & Tang, H. X.  Multimode strong coupling in superconducting cavity piezoelectromechanics.  _Phys. Rev. Lett._ 117, 123603, DOI: 10.1103/PhysRevLett.117.123603 (2016).
  * [19] Wu, S.-C. _et al._ Microwave-controlled optical double optomechanically induced transparency in a hybrid piezo-optomechanical cavity system.  _Phys. Rev. A_ 98, 013807, DOI: 10.1103/PhysRevA.98.013807 (2018).
  * [20] Stockill, R. _et al._ Gallium phosphide as a piezoelectric platform for quantum optomechanics.  _Phys. Rev. Lett._ 123, 163602, DOI: 10.1103/PhysRevLett.123.163602 (2019).
  * [21] Barg, A. _et al._ Carrier-mediated optomechanical forces in semiconductor nanomembranes with coupled quantum wells.  _Phys. Rev. B_ 98, 155316, DOI: 10.1103/PhysRevB.98.155316 (2018).
  * [22] Dahmani, Y. D., Sarabalis, C. J., Jiang, W., Mayor, F. M. & Safavi-Naeini, A. H.  Piezoelectric transduction of a wavelength-scale mechanical waveguide.  _Phys. Rev. Applied_ 13, 024069, DOI: 10.1103/PhysRevApplied.13.024069 (2020).
  * [23] Han, X. _et al._ Cavity piezo-mechanics for superconducting-nanophotonic quantum interface.  _Nature Communications_ 11, 3237, DOI: 10.1038/s41467-020-17053-3 (2020).
  * [24] Zou, C.-L., Han, X., Jiang, L. & Tang, H. X.  Cavity piezomechanical strong coupling and frequency conversion on an aluminum nitride chip.  _Phys. Rev. A_ 94, 013812, DOI: 10.1103/PhysRevA.94.013812 (2016).
  * [25] Shumeiko, V. S.  Quantum acousto-optic transducer for superconducting qubits.  _Phys. Rev. A_ 93, 023838, DOI: 10.1103/PhysRevA.93.023838 (2016).
  * [26] Arrangoiz-Arriola, P. & Safavi-Naeini, A. H.  Engineering interactions between superconducting qubits and phononic nanostructures.  _Phys. Rev. A_ 94, 063864, DOI: 10.1103/PhysRevA.94.063864 (2016).
  * [27] Wang, X., Qin, W., Miranowicz, A., Savasta, S. & Nori, F.  Unconventional cavity optomechanics: Nonlinear control of phonons in the acoustic quantum vacuum.  _Phys. Rev. A_ 100, 063827, DOI: 10.1103/PhysRevA.100.063827 (2019).
  * [28] Zhong, C., Han, X., Tang, H. X. & Jiang, L.  Entanglement of microwave-optical modes in a strongly coupled electro-optomechanical system.  _Phys. Rev. A_ 101, 032345, DOI: 10.1103/PhysRevA.101.032345 (2020).
  * [29] Weiß, M. _et al._ Multiharmonic frequency-chirped transducers for surface-acoustic-wave optomechanics.  _Phys. Rev. Applied_ 9, 014004, DOI: 10.1103/PhysRevApplied.9.014004 (2018).
  * [30] Yamazaki, R. _et al._ Radio-frequency-to-optical conversion using acoustic and optical whispering-gallery modes.  _Phys. Rev. A_ 101, 053839, DOI: 10.1103/PhysRevA.101.053839 (2020).
  * [31] Javerzac-Galy, C. _et al._ On-chip microwave-to-optical quantum coherent converter based on a superconducting resonator coupled to an electro-optic microresonator.  _Phys. Rev. A_ 94, 053815, DOI: 10.1103/PhysRevA.94.053815 (2016).
  * [32] Balram, K. C., Davanço, M. I., Song, J. D. & Srinivasan, K.  Coherent coupling between radiofrequency, optical and acoustic waves in piezo-optomechanical circuits.  _Nature Photonics_ 10, 346–352, DOI: 10.1038/nphoton.2016.46 (2016).
  * [33] Jiang, W. _et al._ Efficient bidirectional piezo-optomechanical transduction between microwave and optical frequency.  _Nature Communications_ 11, 1166, DOI: 10.1038/s41467-020-14863-3 (2020).
  * [34] Midolo, L., Schliesser, A. & Fiore, A.  Nano-opto-electro-mechanical systems.  _Nature Nanotechnology_ 13, 11–18, DOI: 10.1038/s41565-017-0039-1 (2018).
  * [35] Rueda, A., Hease, W., Barzanjeh, S. & Fink, J. M.  Electro-optic entanglement source for microwave to telecom quantum state transfer.  _npj Quantum Information_ 5, 108, DOI: 10.1038/s41534-019-0220-5 (2019).
  * [36] Sychev, D. V. _et al._ Entanglement and teleportation between polarization and wave-like encodings of an optical qubit.  _Nature Communications_ 9, 3672, DOI: 10.1038/s41467-018-06055-x (2018).
  * [37] Tavernarakis, A. _et al._ Optomechanics with a hybrid carbon nanotube resonator.  _Nature Communications_ 9, 662, DOI: 10.1038/s41467-018-03097-z (2018).
  * [38] Wang, C. _et al._ Integrated lithium niobate electro-optic modulators operating at cmos-compatible voltages.  _Nature_ 562, 101–104, DOI: 10.1038/s41586-018-0551-y (2018).
  * [39] Witmer, J. D. _et al._ High-q photonic resonators and electro-optic coupling using silicon-on-lithium-niobate.  _Scientific Reports_ 7, 46313, DOI: 10.1038/srep46313 (2017).
  * [40] Nielsen, M. A. & Chuang, I.  Quantum computation and quantum information.  _American Journal of Physics_ 70, 558–559, DOI: 10.1119/1.1463744 (2002).  https://doi.org/10.1119/1.1463744.
  * [41] Breuer, H.-P. & Petruccione, F.  _The theory of open quantum systems_ (Oxford University Press, Oxford [England]; New York, 2010).
  * [42] Chen, Y., You, J. Q. & Yu, T.  Non-markovian quantum interference in multilevel quantum systems: exact master equation approach.  _Quantum Inf. Comput._ 18, 1261–1271 (2018).
  * [43] Qiao, Y., Zhang, J., Chen, Y., Jing, J. & Zhu, S.  Quantumness protection for open systems in a double-layer environment.  _Science China Physics, Mechanics & Astronomy_ 63, 250312, DOI: 10.1007/s11433-019-1459-2 (2019).
  * [44] Chen, Y., Ding, Q., Shi, W., Jun, J. & Yu, T.  Exact entanglement dynamics mediated by leaky optical cavities.  _Journal of Physics B: Atomic, Molecular and Optical Physics_ 53, 125501, DOI: 10.1088/1361-6455/ab707c (2020).
  * [45] Diósi, L., Gisin, N. & Strunz, W. T.  Non-markovian quantum state diffusion.  _Phys. Rev. A_ 58, 1699–1712, DOI: 10.1103/PhysRevA.58.1699 (1998).
  * [46] Chen, Y. & Yu, T.  _Non-Markovian Dynamics of Qubit Systems: Quantum-State Diffusion Equations Versus Master Equations_ , 609–634 (Springer International Publishing, Cham, 2016).
  * [47] Chen, Y., You, J. Q. & Yu, T.  Exact non-markovian master equations for multiple qubit systems: Quantum-trajectory approach.  _Physical Review A_ 90, 052104, DOI: 10.1103/PhysRevA.90.052104 (2014).
  * [48] Ma, T., Chen, Y., Chen, T., Hedemann, S. R. & Yu, T.  Crossover between non-markovian and markovian dynamics induced by a hierarchical environment.  _Phys. Rev. A_ 90, 042108, DOI: 10.1103/PhysRevA.90.042108 (2014).

## Acknowledgements

Y.C. gratefully acknowledges financial supports by the Institutional Support
of Research and Creativity (ISRC) grant provided by New York Institute of
Technology.

## Author contributions statement

Y.C. designed the study and derived the mathematical model. Q.D. and P.Z.
performed numerical simulations and prepared the manuscript under the guidance
of Y.C. and Y. M. All authors have read and agreed to the published version of
the manuscript.

## Additional information

Competing Interests: The authors declare no competing interests.


#  mlr3summary: Concise and interpretable summaries for machine learning
models

Susanne Dandl
LMU Munich
MCML
&Marc Becker
LMU Munich
&Bernd Bischl
LMU Munich
MCML
0000-0003-4324-4163 0000-0002-8115-0400 0000-0001-6002-6980 Giuseppe
Casalicchio
LMU Munich
MCML &Ludwig Bothmann
LMU Munich
MCML  0000-0001-5324-5966 0000-0002-1471-6582<EMAIL_ADDRESS>muenchen.de

Susanne Dandl, Marc Becker, Bernd Bischl, Giuseppe Casalicchio, Ludwig
Bothmann mlr3summary: Concise and interpretable summaries for machine learning
models mlr3summary: Concise and interpretable summaries for machine learning
models This work introduces a novel R package for concise, informative
summaries of machine learning models. We take inspiration from the summary
function for (generalized) linear models in R, but extend it in several
directions: First, our summary function is model-agnostic and provides a
unified summary output also for non-parametric machine learning models;
Second, the summary output is more extensive and customizable – it comprises
information on the dataset, model performance, model complexity, model’s
estimated feature importances, feature effects, and fairness metrics; Third,
models are evaluated based on resampling strategies for unbiased estimates of
model performances, feature importances, etc. Overall, the clear, structured
output should help to enhance and expedite the model selection process, making
it a helpful tool for practitioners and researchers alike. Model summary,
interpretable machine learning, resampling-based evaluation Model summary,
interpretable machine learning, resampling-based evaluation Ludwig Bothmann
Institut für Statistik
Ludwig-Maximilians-Universität München, Germany
Ludwigstr. 33, 80539 Munich, Germany
Munich Center for Machine Learning (MCML), Germany
E-mail:

##  1 Introduction

Machine learning (ML) increasingly supports decision-making processes in
various domains. A data scientist has a wide range of models available,
ranging from intrinsically interpretable models such as linear models to
highly complex models such as random forests or gradient boosted trees.
Intrinsically interpretable models can come at the expense of generalization
performance, i.e., the model’s capability to predict accurately on future
data. Being able to interpret predictive models is either often a strict
requirement for scientific inference or at least a very desirable property to
audit models in other (more technical) contexts. Many methods have been
proposed for interpreting black-box ML models in the field of interpretable ML
(IML).

For comparing (generalized) linear models (GLMs), the stats package in R
offers a summary function, which only requires the model (fitted with lm or
glm) as input. As an example, glm is applied to a preprocessed version of the
German credit dataset (Hofmann, 1994) (available in the package via
data("credit", package = "mlr3summary")):

⬇

> logreg = glm(risk ~., data = credit,

+ family = binomial(link = "logit"))

> summary(logreg)

    ¯Call:
    ¯glm(formula = risk ~., data = credit, family = binomial(link = "logit"))
    ¯Coefficients:
    ¯Estimate Std. Error z value Pr(>|z|)
    ¯(Intercept)        1.057e+00  3.646e-01   2.900  0.00373 **
    ¯age                9.103e-03  8.239e-03   1.105  0.26925
    ¯...
    ¯Residual deviance: 656.19  on 515  degrees of freedom
    ¯AIC: 670.19
    ¯...

This (shortened) summary informs about the significance of variables
(Pr(>|z|)), their respective effect size and direction (Estimate), as well as
the goodness-of-fit of the model (Residual deviance and AIC). Unfortunately,
many other non-parametric ML models currently cannot be analyzed similarly:
either targeted implementations exist for specific model classes, or an array
of different model-agnostic interpretability techniques (e.g., to derive
feature importance) scattered across multiple packages (Molnar _et al._ ,
2018; Biecek, 2018; Zargari Marandi, 2023) must be employed. However,
especially in applied data science, a user often performs model selection or
model comparison across an often diverse pool of candidate models, so a
standardized diagnostic output becomes highly desirable.

Another issue is that in the glm-based summary, the goodness-of-fit is only
evaluated on the training data, but not on hold-out/test data. While this
might be appropriate for GLM-type models – provided proper model diagnosis has
been performed – this is not advisable for non-parametric and non-linear
models, which can overfit the training data.111For completeness’ sake:
Overfitting can happen for GLMs, e.g., in high-dimensional spaces with limited
sample size. Here, hold-out test data or in general resampling techniques like
cross-validation should be used for proper estimation of the generalization
performance Simon (2007). Such resampling-based performance estimation should
also be used for loss-based IML methods. For interpretability methods that
only rely on predictions, this might also be advisable but might not lead to
huge differences in results (Molnar _et al._ , 2022, 2023).

#### Contributions

With the mlr3summary package, we provide a novel model-agnostic summary
function for ML models and learning algorithms in R. This is facilitated by
building upon mlr3 (Lang _et al._ , 2019; Bischl _et al._ , 2024) – a package
ecosystem for applied ML, including resampling-based performance assessment.
The summary function returns a structured overview that gives information on
the underlying dataset and model, generalization performances, complexity of
the model, fairness metrics, and feature importances and effects. For the
latter two, the function relies on model-agnostic methods from the field of
IML. The output is customizable via a flexible control argument to allow
adaptation to different application scenarios. The mlr3summary package is
released under LGPL-3 on GitHub (https://github.com/mlr-org/mlr3summary) and
CRAN (https://cran.r-project.org/package=mlr3summary). Documentations in the
form of help pages are available as well as unit tests. The example code of
this manuscript is available via demo("credit", package = "mlr3summary").

##  2 Related work

Most R packages that offer model summaries are restricted to parametric models
and extend the stats summary method (e.g., modelsummary (Arel-Bundock, 2022),
broom (Robinson _et al._ , 2023)). Performance is only assessed based on
training data – generalization errors are not provided. Packages that can
handle diverse ML models focus primarily on performance assessment (e.g., mlr3
(Lang _et al._ , 2019), caret (Kuhn and Max, 2008)). Packages that primarily
consider feature importances and effects do not provide overviews in a
concise, decluttered format but provide extensive reports (e.g., modelDown
(Romaszko _et al._ , 2019) and modelStudio (Baniecki and Biecek, 2019) based
on DALEX (Biecek, 2018), or explainer (Zargari Marandi, 2023)). While it is
possible to base the assessment on hold-out/test data, assessment based on
resampling is not automatically supported by these packages. Overall, to the
best of our knowledge, there is no R package yet that allows for a concise yet
informative overview based on resampling-based performance assessment, model
complexity, feature importance and effect directions, and fairness metrics.

##  3 Design, functionality, and example

The core function of the mlr3summary package is the S3-based summary function
for mlr3 Learner objects. It has three arguments: object reflects a trained
model – a model of class Learner fitted with mlr3; resample_result reflects
the results of resampling – a ResampleResult object fitted with mlr3; control
reflects some control arguments – a list created with summary_control (details
in Section 3.2).

The mlr3 package is the basis of mlr3summary because it provides a unified
interface to diverse ML models and resampling strategies. A general overview
of the mlr3 ecosystem is given in Bischl _et al._ (Bischl _et al._ , 2024).
With mlr3, the modelling process involves the following steps: (1) initialize
a regression or classification task, (2) choose a regression or classification
learner, (3) train a model with the specified learner on the initialized task,
(4) apply a resampling strategy. The last step is necessary to receive valid
estimates for performances, importances, etc., as mentioned in Section 1. The
following lines of code illustrate steps (1)-(4) on the (preprocessed) credit
dataset from Section 1 using a ranger random forest. As a resampling strategy,
we conduct $3$-fold cross-validation.

⬇

> task = TaskClassif$new(id = "credit", backend = credit,

+ target = "risk")

> rf = lrn("classif.ranger", predict_type = "prob")

> rf$train(task)

> cv3 = rsmp("cv", folds = 3L)

> rr = resample(task = task, learner = rf, resampling = cv3,

+ store_models = TRUE)

Internally, the resample function fits, in each iteration, the model on the
respective training data, uses the model to predict the held-out test data,
and stores the predictions in the result object. To compute performances,
complexities, importances, and other metrics, the summary function iteratively
accesses the models and datasets within the resulting resample object, which
requires setting the parameter store_models = TRUE within the resample
function. For the final summary output, the results of each iteration are
aggregated (e.g., averages and standard deviations (sds)).

###  3.1 Summary function and output

This section shows the summary call and output for the random forest of the
previous credit example and provides some details on each displayed paragraph.

⬇

> summary(object = rf, resample_result = rr)

General provides an overview of the task, the learner (including its
hyperparameters), and the resampling strategy.222Currently, this is the only
paragraph that is based on object, all other paragraphs are based on
resample_result.

Residuals display the distribution of residuals of hold-out data over the
resampling iterations. For regression models, the residuals display the
difference between true and predicted outcome. For classifiers that return
class probabilities, the residuals are defined as the difference between
predicted probabilities and a one-hot-encoding of the true class. For
classifiers that return classes, a confusion matrix is shown.

Performance displays averages and sds (in $[\,]$) of performance measures over
the iterations.333Please note that there is no unbiased estimator of the
variance, see (Nadeau and Bengio, 1999) and Section 5 for a discussion. The
shown performance values are the area-under-the-curve (auc), the F-score
(fbeta), the binary Brier score (bbrier), and Mathew’s correlation coefficient
(mcc).

The arrows display whether lower or higher values refer to a better
performance. “(macro)” indicates a macro aggregation, i.e., measures are
computed for each iteration separately before averaging. “(micro)” would
indicate that measures are computed across all iterations (see Bischl _et al._
(2024) for details).

Complexity displays averages and sds of two model complexity measures proposed
by Molnar _et al._ (Molnar _et al._ , 2020): sparsity shows the number of used
features that have a non-zero effect on the prediction (evaluated by
accumulated local effects (ale) (Apley and Zhu, 2020)); interaction_strength
shows the scaled approximation error between a main effect model (based on
ale) and the prediction function.444The interaction strength has a value in
$[0,1]$, 0 means no interactions, 1 means no main effects but interactions.

Importance shows the averages and sds of feature importances over the
iterations. The first column (pdp) displays importances based on the sds of
partial dependence curves (Friedman, 2001; Greenwell _et al._ , 2018), the
second column (pfi.ce) shows the results for permutation feature importance
Breiman (2001); Fisher _et al._ (2019).

Effects shows average effect plots over the iterations – partial dependence
plots (pdp) and ale plots (Friedman, 2001; Apley and Zhu, 2020). For binary
classifiers, the effect plots are only shown for the positively-labeled class
(here, task$positive = "good"). For multi-class classifiers, the effect plots
are given for each outcome class separately (one vs. all). For categorical
features, the bars are ordered according to the factor levels of the feature.

The learner can also be a complete pipeline from mlr3pipelines (Binder _et
al._ , 2021), where the most common case would be an ML model with associated
pre-processing steps. Then, the summary output also shows some basic
information about the pipeline.555Linear pipelines can be displayed in the
console, non-linear parts are suppressed in the output. Since preprocessing
steps are treated as being part of the learner, the summary output is
displayed on the original data (e.g., despite one-hot encoding of categorical
features, importance results are not shown for each encoding level
separately). The learner can also be an AutoTuner from mlr3tuning, where
automatic processes for tuning the hyperparameters are conducted. Examples on
pipelining and tuning are given in the demo of the package.

###  3.2 Customizations

The output of the summary function can be customized via a control argument
which requires a list created with the function summary_control as an input.
If no control is specified, the following default setting is used:

⬇

> summary_control(measures = NULL,

+ complexity_measures = c("sparsity", "interaction_strength"),

+ importance_measures = NULL, n_important = 15L,

+ effect_measures = c("pdp", "ale"),

+ fairness_measures = NULL, protected_attribute = NULL,

+ hide = NULL, digits = max(3L, getOption("digits") - 3L))

Performances are adaptable via measures, complexities via complexity_measures,
importances via importance_measures and effects via effect_measures within
summary_control. Examples are given in the demo of the package. The default
for measures and importance_ measures is NULL, which results in a collection
of commonly reported measures being chosen, based on the task type – for
concrete measures see the help page (?summary_control). n_important reflects
that, by default, only the $15$ most important features are displayed in the
output. This is especially handy for high-dimensional data. With hide,
paragraphs of the summary output can be omitted (e.g., "performance") and with
digits, the number of printed digits is specified.

Fairness assessment for classification and regression models is also available
in mlr3summary based on the mlr3fairness package (Pfisterer _et al._ , 2023).
Therefore, a protected attribute must be specified. This can be done either
within the task by updating the feature roles or by specifying a
protected_attribute in summary_control. The following shows the code and
output when specifying sex as a protected attribute. The shown default
fairness measures are demographic parity (dp), conditional use accuracy
equality (cuae) and equalized odds (eod), other measures are possible via
fairness_measures in summary_control.

⬇

> summary(object = rf, resample_result = rr,

+ control = summary_control(protected_attribute = "sex"))

##  4 Runtime assessment

To assess how the runtime scales with differing numbers of features
$p\in\\{5,10,25,50,100\\}$ and numbers of observations
$n\in\\{50,100,500,1000,2000\\}$, we conducted a simulation study. Given
$X_{1},X_{2},X_{3}\sim U(0,1)$, $X_{4}\sim Bern(0.75)$, the data generating
process is $y=f(\mathbf{x})+\epsilon$ with
$f(\mathbf{x})=4x_{1}+4x_{2}+4x_{4}x_{3}^{2}$ and $\epsilon\sim N(0,0.1\cdot
f(\mathbf{x}))$. As noise variables, $X_{5}$ as a categorical feature with
five classes, and $X_{6},...,X_{p}\sim N(0,1)$ were added to the data. We
trained random forests and linear main effect models on the datasets and
conducted 3-fold cross-validation. The first two figures in Figure 1 show that
runtimes of the linear model were lower compared to the random forest. To
improve runtimes, we added parallelization over the resampling iterations (via
the future package (Bengtsson, 2021)) as another feature to mlr3summary –
results for the random forest (with $3$ cores) are on the right. Overall,
scaling of runtimes is worse in $p$ than in $n$.

Figure 1: Runtimes of the summary function for linear models (left), and
random forests without (middle) and with (right) parallelization, for
differing numbers of features $p$ and observations $n$.

##  5 Outlook and discussion

In conclusion, this paper introduces a novel R package for concise model
summaries. The summary output is highly adaptable due to a control argument
and might be extended in the future. We also plan to offer a report function
for detailed visualizations and model comparisons. To assess importance and
effects of single features, mlr3summary builds upon the iml and fastshap
packages. These packages only offer a limited set of interpretation methods.
Recommended alternatives to permutation feature importances like conditional
feature importance Molnar _et al._ (2022), are currently not available in a
proper R package (published on CRAN). Our summary also currently lacks proper
statistical tests for importances or confidence intervals for performances.
This is because unbiased estimates of the variance are required which is a
challenge for resampling strategies and the available methods that propose
unbiased estimates are computationally infeasible (e.g., due to many model
refits) (Molnar _et al._ , 2023; Stephen Bates and Tibshirani, 2023).
Addressing this issue requires some concerted efforts from the research
community. If methods are readily available in R, we are happy to integrate
them in mlr3summary.

## Acknowledgments

This work has been partially supported by the Federal Statistical Office of
Germany.

## References

  * Apley and Zhu (2020) Apley DW, Zhu J (2020).  “Visualizing the Effects of Predictor Variables in Black Box Supervised Learning Models.”  _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 82(4), 1059–1086.  10.1111/rssb.12377.
  * Arel-Bundock (2022) Arel-Bundock V (2022).  “modelsummary: Data and Model Summaries in R.”  _Journal of Statistical Software_ , 103(1), 1–23.  10.18637/jss.v103.i01.
  * Baniecki and Biecek (2019) Baniecki H, Biecek P (2019).  “modelStudio: Interactive Studio with Explanations for ML Predictive Models.”  _Journal of Open Source Software_ , 4(43), 1798.  10.21105/joss.01798.
  * Bengtsson (2021) Bengtsson H (2021).  “A Unifying Framework for Parallel and Distributed Processing in R using Futures.”  _The R Journal_ , 13(2), 208–227.  10.32614/RJ-2021-048.
  * Biecek (2018) Biecek P (2018).  “DALEX: Explainers for Complex Predictive Models in R.”  _Journal of Machine Learning Research_ , 19(84), 1–5.
  * Binder _et al._ (2021) Binder M, Pfisterer F, Lang M, Schneider L, Kotthoff L, Bischl B (2021).  “mlr3pipelines - Flexible Machine Learning Pipelines in R.”  _Journal of Machine Learning Research_ , 22(184), 1–7.
  * Bischl _et al._ (2024) Bischl B, Sonabend R, Kotthoff L, Lang M (2024).  _Applied machine learning using mlr3 in R_.  Chapman and Hall/CRC.  ISBN 9781003402848.  10.1201/9781003402848.
  * Breiman (2001) Breiman L (2001).  “Random Forests.”  _Machine Learning_ , 45(1), 5–32.  ISSN 0885-6125.  10.1023/a:1010933404324.
  * Fisher _et al._ (2019) Fisher A, Rudin C, Dominici F (2019).  “All Models are Wrong, but Many are Useful: Learning a Variable’s Importance by Studying an Entire Class of Prediction Models Simultaneously.”  _Journal of Machine Learning Research_ , 20(177).
  * Friedman (2001) Friedman JH (2001).  “Greedy Function Approximation: A Gradient Boosting Machine.”  _The Annals of Statistics_ , 29(5).  ISSN 0090-5364.  10.1214/aos/1013203451.
  * Greenwell _et al._ (2018) Greenwell BM, Boehmke BC, McCarthy AJ (2018).  “A Simple and Effective Model-Based Variable Importance Measure.”  _arXiv preprint arXiv:1805.04755_.  10.48550/arXiv.1805.04755.
  * Hofmann (1994) Hofmann H (1994).  “Statlog (German Credit Data).”  _UCI Machine Learning Repository_.  10.24432/C5NC77.
  * Kuhn and Max (2008) Kuhn, Max (2008).  “Building Predictive Models in R Using the caret Package.”  _Journal of Statistical Software_ , 28(5), 1–26.  10.18637/jss.v028.i05.
  * Lang _et al._ (2019) Lang M, Binder M, Richter J, Schratz P, Pfisterer F, Coors S, Au Q, Casalicchio G, Kotthoff L, Bischl B (2019).  “mlr3: A Modern Object-oriented Machine Learning Framework in R.”  _Journal of Open Source Software_.  10.21105/joss.01903.
  * Molnar _et al._ (2018) Molnar C, Bischl B, Casalicchio G (2018).  “iml: An R Package for Interpretable Machine Learning.”  _JOSS_ , 3(26), 786.  10.21105/joss.00786.
  * Molnar _et al._ (2020) Molnar C, Casalicchio G, Bischl B (2020).  _Quantifying model complexity via functional decomposition for better post-hoc interpretability_ , p. 193–204.  Springer International Publishing.  10.1007/978-3-030-43823-4_17.
  * Molnar _et al._ (2023) Molnar C, Freiesleben T, König G, Herbinger J, Reisinger T, Casalicchio G, Wright MN, Bischl B (2023).  “Relating the Partial Dependence Plot and Permutation Feature Importance to the Data Generating Process.”  In L Longo (ed.), _Explainable Artificial Intelligence_ , pp. 456–479. Springer Nature Switzerland, Cham.  10.1007/978-3-031-44064-9_24.
  * Molnar _et al._ (2022) Molnar C, König G, Herbinger J, Freiesleben T, Dandl S, Scholbeck CA, Casalicchio G, Grosse-Wentrup M, Bischl B (2022).  “General Pitfalls of Model-Agnostic Interpretation Methods for Machine Learning Models.”  In A Holzinger, et al (eds.), _xxAI - Beyond Explainable AI: International Workshop_ , pp. 39–68. Springer International Publishing, Cham.  ISBN 978-3-031-04083-2.  10.1007/978-3-031-04083-2_4.
  * Nadeau and Bengio (1999) Nadeau C, Bengio Y (1999).  “Inference for the Generalization Error.”  In S Solla, et al (eds.), _Advances in Neural Information Processing Systems_ , volume 12, pp. 1–7. MIT Press.
  * Pfisterer _et al._ (2023) Pfisterer F, Siyi W, Lang M (2023).  _mlr3fairness: Fairness auditing and debiasing for ’mlr3’_.  R package version 0.3.2, URL https://CRAN.R-project.org/package=mlr3fairness.
  * Robinson _et al._ (2023) Robinson D, Hayes A, Couch S (2023).  _broom: Convert statistical objects into tidy tibbles_.  R package version 1.0.5, URL https://CRAN.R-project.org/package=broom.
  * Romaszko _et al._ (2019) Romaszko K, Tatarynowicz M, Urbański M, Biecek P (2019).  “modelDown: Automated Website Generator with Interpretable Documentation for Predictive Machine Learning Models.”  _Journal of Open Source Software_ , 4(38).  10.21105/joss.01444.
  * Simon (2007) Simon R (2007).  _Resampling strategies for model assessment and selection_ , pp. 173–186.  Springer US, Boston, MA.  10.1007/978-0-387-47509-7_8.
  * Stephen Bates and Tibshirani (2023) Stephen Bates TH, Tibshirani R (2023).  “Cross-Validation: What Does It Estimate and How Well Does It Do It?”  _Journal of the American Statistical Association_ , pp. 1–12.  10.1080/01621459.2023.2197686.
  * Zargari Marandi (2023) Zargari Marandi R (2023).  _explainer: Machine learning model explainer_.  R package version 1.0.0, URL https://CRAN.R-project.org/package=explainer.