# Chiral current in Floquet cavity-magnonics

Shi-fan Qi  School of Physics, Zhejiang University, Hangzhou 310027, Zhejiang,
China  Jun Jing<EMAIL_ADDRESS>School of Physics, Zhejiang University,
Hangzhou 310027, Zhejiang, China

###### Abstract

Floquet engineering can induce complex collective behaviour and interesting
synthetic gauge-field in quantum systems through temporal modulation of system
parameters by periodic drives. Using a Floquet drive on frequencies of the
magnon modes, we realize a chiral state-transfer in a cavity-magnonic system.
The time-reversal symmetry is broken in such a promising platform for coherent
information processing. In particular, the photon mode is adiabatically
eliminated in the large-detuning regime and the magnon modes under conditional
longitudinal drives can be indirectly coupled to each other with a phase-
modulated interaction. The effective Hamiltonian is then used to generate
chiral currents in a circular loop, whose dynamics is evaluated to measure the
symmetry of the system Hamiltonian. Beyond the dynamics limited in the
manifold with a fixed number of excitations, our protocol applies to the
continuous-variable systems with arbitrary states. Also it is found to be
robust against the systematic errors in the photon-magnon coupling strength
and Kerr nonlinearity.

##  I Introduction

With unique properties observed in experiments such as large tunability and
long coherent-time, the magnon can be used as information carrier in a broad
variety of hybrid systems Rameshti  _et al._ (2021); Lachance-Quirion  _et
al._ (2019); Li  _et al._ (2020); Yuan  _et al._ (2022). Found to be directly
or effectively coupled to microwave photon Soykal and Flatté (2010a, b);
Tabuchi  _et al._ (2014); Zhang  _et al._ (2014); Wang  _et al._ (2018); Shen
_et al._ (2021), superconducting qubit Lachance-Quirion  _et al._ (2020);
Tabuchi  _et al._ (2015); Qi and Jing (2022a); Yan and Jing (2021), spin-
emitter Neuman  _et al._ (2020); Qi and Jing (2021a), and phonon Zhang  _et
al._ (2016); Li  _et al._ (2018a); Qi and Jing (2021b); Li  _et al._ (2021);
Hatanaka  _et al._ (2022); Potts  _et al._ (2021); Shen  _et al._ (2022), the
magnon modes provide novel applications in quantum computing Ladd  _et al._
(2010), quantum communication Reiserer and Rempe (2015), and quantum sensing
Degen  _et al._ (2017). Fundamentally, the coupling between magnon and other
types of mode motivates a number of insights, e.g., the indirect two-body
flip-flop interaction of magnons for Rabi oscillation, which is symmetrical to
time and parity. Moreover, Floquet engineering has been applied to a cavity
electromagnonic system Xu  _et al._ (2020). It could support more flexible and
interesting applications based on the complex-valued hopping, such as
generating chiral transfers of the quasiparticles of spin waves as a highly
promising bosonic platform.

Chirality has played an important role in the photon and phonon systems Lodahl
_et al._ (2017); Zhu  _et al._ (2018). Optical lattice could be used to
simulate the chirality in magnetic materials Sørensen  _et al._ (2005). In the
platform of superconducting qubits Liu  _et al._ (2020), spin chirality arises
from the three-body interactions. The chirality operator of three spins is
$\hat{O}=\vec{\sigma}_{1}\cdot(\vec{\sigma}_{2}\times\vec{\sigma}_{3})$, where
$\vec{\sigma}_{j}\equiv(\sigma^{x}_{j},\sigma^{y}_{j},\sigma^{z}_{j})$ is the
Pauli vector for the $j$th spin particle. It is straightforward to verify that
$\hat{O}$ breaks both the time reversal symmetry $\mathcal{T}$ (replacing
$\vec{\sigma}_{j}$ by $-\vec{\sigma}_{j}$) and the parity symmetry
$\mathcal{P}$ (exchanging $\vec{\sigma}_{j}$ with $\vec{\sigma}_{k}$) but
conserves the $\mathcal{PT}$ symmetry Liu  _et al._ (2020). The dynamics of
the spins driven by $\hat{O}$ features a chiral evolution, i.e.,
$|s_{1}s_{2}s_{3}\rangle\to$ $|s_{2}s_{3}s_{1}\rangle\to$
$|s_{3}s_{1}s_{2}\rangle$, where $s_{j}=0$ (spin down) or $1$ (spin up) Cai
and Wang (2021). In parallel to the protocol in the superconducting qubits, we
propose to generate a chiral transfer of an arbitrary bosonic state in the
Floquet cavity-magnonic system Xu  _et al._ (2020). Here the system-
Hamiltonian symmetry is described by a nonvanishing chiral current operator.

According to the Floquet theory Shirley (1965); Goldman and Dalibard (2014);
Bukov  _et al._ (2015); Petiziol  _et al._ (2021), the propagator $U(t)$
induced by a time-dependent Hamiltonian $H(t)=H(t+T)$ with a period $T$ can be
expressed by $U(t)=U_{F}(t)\exp(-iH^{\prime}t)$, where $U_{F}(t)$ is also
periodic with $T$. As a time-independent Hermitian operator,
$H^{\prime}=U^{\dagger}_{F}(t)H(t)U_{F}(t)-iU^{\dagger}_{F}(t)\dot{U}_{F}(t)$
can be understood as the system Hamiltonian in the frame defined by
$U_{F}(t)$. Floquet engineering by fast periodic modulation over the
characteristic frequencies of a quantum system is a major control approach to
the long-time dynamics of the system Goldman and Dalibard (2014); Shao  _et
al._ (2017); Wang  _et al._ (2016a); Li  _et al._ (2019), where the effective
Hamiltonian can be synthesized through designing proper drives. Floquet
engineering has been applied to quantum switch Wu  _et al._ (2018), chiral
ground state current Roushan1  _et al._ (2017), quantum simulation Kyriienko
and Sørensen (2018), and perfect state transfer Wang  _et al._ (2019); Liu
_et al._ (2020); Li  _et al._ (2018b) in superconducting qubits. Also it found
a very wide application in quantum gas systems Eckardt (2017); Weitenberg and
Simonet (2021) to reveal many interesting physics.

In this work, we apply the Floquet engineering to a hybrid photon-magnon
system in which three single-crystal yttrium iron garnet (YIG) spheres are
placed inside a microwave cavity. The uniform bias magnetic field excites the
Kittel mode in the YIG spheres and establishes a strong photon-magnon
coupling. In the large detuning regime, the common photon mode can be
adiabatically eliminated with the standard perturbation theory Garziano  _et
al._ (2016); Ma and Law (2015); Kaufman  _et al._ (2020) or the high-order
Fermi golden rule Combescot (2001). By periodically modulating over the three
magnon modes with well-controlled intensities, frequencies, and phases, we can
obtain an effective time-reversal symmetry broken Hamiltonian Koch  _et al._
(2010), that ensures chiral magnon currents Roushan1  _et al._ (2017) of
arbitrary states with a state-independent period. In contrast to the protocols
based on the superconducting qubits Roushan1  _et al._ (2017); Wang  _et al._
(2019); Liu  _et al._ (2020), our protocol is not constrained in a manifold of
a fixed number of excitations and is capable to transfer versatile states such
as Fock state, coherent state, Schrödinger cat state, and even two-body
entangled states in a chiral way. The magnon Fock state can be generated by
the magnon blockade Xie  _et al._ (2020); Liu  _et al._ (2019). A magnonic cat
state was proposed in Ref. Kounalakis  _et al._ (2022). And more magnonic-
state preparation protocols have been collected in a recent review Yuan  _et
al._ (2022).

The rest part of this work is structured as follows. In Sec. II, we introduce
the quantum model for a Floquet cavity-magnonic system and then derive an
effective Hamiltonian to generate a perfect chiral state transfer amongst
magnons. Clockwise and anticlockwise transfers of various states and chiral
currents to measure the symmetry of the system Hamiltonian are presented in
Sec. III.1 and Sec. III.2, respectively. State-transfer fidelity under the
magnon damping and the systematic errors arise from the photon-magnon coupling
strength and the Kerr effects are discussed in Sec. IV.1 and Sec. IV.2,
respectively. The whole work is summarized in Sec. V.

##  II Model

Consider a hybrid quantum model consists of a microwave cavity coupled to $N$
YIG spheres in their Kittle modes Xu  _et al._ (2020); Wang  _et al._ (2018,
2016b). The model Hamiltonian reads

$H=\omega_{a}a^{\dagger}a+\omega_{m}\sum_{k=1}^{N}m^{\dagger}_{k}m_{k}+g_{am}\sum_{k=1}^{N}\left(am^{\dagger}_{k}+a^{\dagger}m_{k}\right),$
(1)

where $a$ $(a^{\dagger})$ and $m_{k}$ $(m^{\dagger}_{k})$ are the annihilation
(creation) operators of the photon and the $k$th magnon modes, respectively.
$\omega_{a}$ and $\omega_{m}$ are their respective transition frequencies.
$g_{am}$ is the single-excitation coupling strength between photon and magnon
modes, which is much larger than their decay rates and much smaller than their
detuning $|\omega_{a}-\omega_{m}|$ in the dispersive regime. Counter-rotating
interaction between photon and magnons is omitted under the rotating-wave
approximation (RWA), as justified in Appendix A.

When one focuses on the state exchange within any pair of magnon modes, it is
instructive to show that the photon mode could be adiabatically eliminated by
the second-order perturbation Garziano  _et al._ (2016). The subsystem
Hamiltonian associated with the $k$th and $j$th magnon modes and the photon
mode can be written as $H^{kj}_{0}+H^{kj}_{I}$, where

$\displaystyle H^{kj}_{0}$
$\displaystyle=\omega_{a}a^{\dagger}a+\omega_{m}\left(m_{k}^{\dagger}m_{k}+m_{j}^{\dagger}m_{j}\right),$
(2) $\displaystyle H^{kj}_{I}$
$\displaystyle=g_{am}\left[a\left(m^{\dagger}_{k}+m^{\dagger}_{j}\right)+a^{\dagger}\left(m_{k}+m_{j}\right)\right]$

indicating the unperturbed and perturbation Hamiltonians, respectively. To the
second order in the perturbation Hamiltonian, the effective coupling strength
$g_{pq}$ between any pair of eigenstates $|p\rangle$ and $|q\rangle$ of
$H^{kj}_{0}$ in Eq. (2) is given by

$g_{pq}=\sum_{w\neq p,q}\frac{\langle q|H^{kj}_{I}|w\rangle\langle
w|H^{kj}_{I}|p\rangle}{E_{p}-E_{w}},$ (3)

where $E_{w}$ is the eigenenergy of an intermediate eigenstate $|w\rangle$,
namely, $H^{kj}_{0}|w\rangle=E_{w}|w\rangle$. Note that the large detuning
condition $\omega_{m}\neq\omega_{a}$ lifts the degeneracy of $H^{kj}_{0}$. Due
to the excitation conservation in both Eqs. (1) and (2), $g_{pq}$ can be
probed as the transition rate between
$|p\rangle=|n_{k}n_{a}n_{j}\rangle\equiv|n_{k}\rangle|n_{a}\rangle|n_{j}\rangle$
and $|q\rangle=|(n_{k}-1)n_{a}(n_{j}+1)\rangle$ in a subspace with a fixed
number of excitations $M=n_{k}+n_{a}+n_{j}$, where the subscript indicates the
mode and $n_{k}$, $n_{a}$, and $n_{j}$ are arbitrary nonnegative integers. By
Eqs. (2) and (3), the leading contribution for
$|n_{k}n_{a}n_{j}\rangle\to|(n_{k}-1)n_{a}(n_{j}+1)\rangle$ comes from two
physical paths, i.e.,
$|n_{k}n_{a}n_{j}\rangle\to|(n_{k}-1)(n_{a}+1)n_{j}\rangle\to|(n_{k}-1)n_{a}(n_{j}+1)\rangle$
and
$|n_{k}n_{a}n_{j}\rangle\to|n_{k}(n_{a}-1)(n_{j}+1)\rangle\to|(n_{k}-1)n_{a}(n_{j}+1)\rangle$.
Then in this subspace we have

$g_{n_{k}n_{j}}=\sqrt{n_{k}(n_{j}+1)}\frac{g^{2}_{am}}{\omega_{m}-\omega_{a}}=g\sqrt{n_{k}(n_{j}+1)}.$
(4)

Note $m_{k}|n_{k}\rangle=\sqrt{n_{k}}|n_{k}-1\rangle$ and
$m_{j}^{\dagger}|n_{j}\rangle=\sqrt{n_{j}+1}|n_{j}+1\rangle$. The effective
Hamiltonian for this pair of magnon modes in their full space Qi and Jing
(2021b) can thus be written as

$H_{kj}=g\left(m_{k}m^{\dagger}_{j}+m^{\dagger}_{k}m_{j}\right),\quad
g\equiv\frac{g^{2}_{am}}{\omega_{m}-\omega_{a}}.$ (5)

The preceding derivation as well as the construction of the effective
Hamiltonian applies to any indirect magnon-magnon interaction mediated by the
common photon mode. In the dispersive regime, the Hamiltonian for the whole
model in Eq. (1) can therefore be written as

$H=g\sum_{k<j}^{N}\left(m_{k}m^{\dagger}_{j}+m^{\dagger}_{k}m_{j}\right).$ (6)

We apply the periodic signals by Floquet driving Xu  _et al._ (2020); Liu  _et
al._ (2020) featured with the intensity $\Delta$, the frequency $\omega$, and
the local phases $\phi_{k}$, to the model Hamiltonian and obtain

$\displaystyle H(t)$ $\displaystyle=\Delta\sum_{k=1}^{N}\cos(\omega
t+\phi_{k})m^{\dagger}_{k}m_{k}$ (7)
$\displaystyle+g\sum_{k<j}^{N}\left(m_{k}m^{\dagger}_{j}+m^{\dagger}_{k}m_{j}\right).$

In the rotating frame with respect to

$\displaystyle U_{0}(t)$
$\displaystyle=\exp\left[i\int^{t}_{0}ds\Delta\sum_{k=1}^{N}\cos(\omega
s+\phi_{k})m^{\dagger}_{k}m_{k}\right]$ (8)
$\displaystyle=\exp\left\\{i\sum^{N}_{k=1}\frac{\Delta}{\omega}\left[\sin(\omega
t+\phi_{k})-\sin\phi_{k}\right]m^{\dagger}_{k}m_{k}\right\\},$

we have

$\displaystyle H_{I}(t)$
$\displaystyle=U_{0}(t)H(t)U^{\dagger}_{0}(t)-iU_{0}(t)\dot{U}^{\dagger}_{0}(t)$
(9) $\displaystyle=$
$\displaystyle\sum_{k<j}^{N}g_{kj}\exp\left[if_{kj}\sin(\omega
t+\alpha_{kj})\right]m_{k}^{\dagger}m_{j}+{\rm H.c.},$

where

$\displaystyle
g_{kj}=ge^{-i\beta_{kj}},\quad\beta_{kj}\equiv\frac{\Delta}{\omega}(\sin\phi_{k}-\sin\phi_{j}),$
$\displaystyle
f_{kj}=\frac{2\Delta}{\omega}\sin\left(\frac{\phi_{j}-\phi_{k}}{2}\right),$
$\displaystyle\alpha_{kj}\equiv\tan^{-1}\left(\frac{\sin\phi_{k}-\sin\phi_{j}}{\cos\phi_{k}-\cos\phi_{j}}\right).$

According to the Jacobi-Anger expansion $e^{ix\sin
y}=\sum_{n=-\infty}^{n=+\infty}J_{n}(x)e^{iny}$, where $J_{n}$ is the $n$th
Bessel function of the first kind, one can obtain

$H_{I}(t)=H_{0}+\sum_{n=1}^{\infty}\left(H_{n}e^{in\omega
t}+H_{-n}e^{-in\omega t}\right),$ (10)

with

$\displaystyle H_{0}$
$\displaystyle=\sum_{k<j}^{N}J_{0}(f_{kj})\left(g_{kj}m^{\dagger}_{k}m_{j}+g^{*}_{kj}m_{k}m_{j}^{\dagger}\right),$
(11) $\displaystyle H_{n}$
$\displaystyle=\sum_{k<j}^{N}J_{n}(f_{kj})e^{in\alpha_{kj}}\left[g_{kj}m^{\dagger}_{k}m_{j}+(-1)^{n}g_{kj}^{*}m_{k}m_{j}^{\dagger}\right],$

and $H_{-n}=H_{n}^{*}$. Therefore up to the order of $\mathcal{O}(1/\omega)$,
the Floquet-driving Hamiltonian (7) can be written as

$H_{\rm eff}=H_{0}+H_{\rm
eff}^{(2)}=H_{0}+\sum^{\infty}_{n=1}\frac{1}{n\omega}\left[H_{n},H_{-n}\right].$
(12)

We have several remarks for $H_{\rm eff}$. (a) The zeroth-order coupling-
strength between the $k$th and $j$th magnons is $g_{kj}J_{0}(f_{kj})$, which
can be tuned by the ratio of the driving intensity and frequency
$\Delta/\omega$ and the phase difference between the local driving signals
$\phi_{k}-\phi_{j}$. (b) When $N=2$, the leading-order correction by the
commutator $[H_{n},H_{-n}]$ vanishes and then the dynamics of magnons becomes
a Rabi oscillation with a frequency fully determined by $f_{kj}$. (c) When
$N>2$, the leading-order correction can be written as

$H_{\rm
eff}^{(2)}=\sum_{k<j}\left(g^{(2)}_{kj}e^{i\alpha_{kj}}m_{k}^{\dagger}m_{j}+h.c.\right),$
(13)

where the second-order coupling strength is

$g^{(2)}_{kj}\sim\frac{g^{2}}{\omega}\sum_{k_{1}\neq
k,j}\left[\sum^{\infty}_{n=1}\frac{J_{n}(f_{kk_{1}})J_{n}(f_{k_{1}j})}{n}\right],$
(14)

with $k_{1}$ indicating the other magnon modes rather than $k$ and $j$ in the
whole system.

The third and even higher orders of terms can be analyzed with the James’
effective-Hamiltonian method Shao  _et al._ (2017). The nonvanishing third-
order terms in the effective Hamiltonian can be obtained by

$\displaystyle H^{(3)}_{kj}$
$\displaystyle=\sum_{l=m+n}\frac{1}{(n+m)\omega^{2}}\left[H_{n},\left[H_{m},H_{-l}\right]\right]$
(15)
$\displaystyle=\sum_{k<j}\left(g^{(3)}_{kj}e^{i\alpha_{kj}}m_{k}^{\dagger}m_{j}+H.c.\right),$

where the third-order coupling strength $g^{(3)}_{kj}$ is about

$\sim\frac{g^{3}}{\omega^{2}}\sum_{k_{1}\neq k,k_{2}\neq
j}\left[\sum_{l=m+n}\frac{J_{n}(f_{kk_{1}})J_{m}(f_{k_{1}k_{2}})J_{l}(f_{k_{2}j})}{n+m}\right].$
(16)

Therefore the third and even higher orders of corrections can be safely
omitted in the dispersive regime $g\ll\omega$.

Chiral state-transfer (as described in Fig. 1) could manifest when $N=3$ by
setting $\phi_{j}=2\pi j/3$, $j=1,2,3$, which are uniformly distributed in the
range of $[0,2\pi]$. In this case, the magnitude of the ratio factor $f_{kj}$
becomes independent of the magnon pair, i.e.,
$f_{12}=f_{13}=f_{23}=\sqrt{3}\Delta/\omega\equiv f$. The relevant quantum
phases are found to be $\beta_{12}=f$, $\beta_{23}=-f/2$, $\beta_{13}=f/2$;
and $\alpha_{12}=\pi/2$, $\alpha_{13}=5\pi/6$, $\alpha_{23}=7\pi/6$, i.e.,
$\alpha_{kj}=(k+j)\pi/3-\pi/2$. Then using $H_{n}$ in Eq. (11), the second-
order term in the effective Hamiltonian (12) turns out to be

$H^{(2)}_{\rm eff}=-ig_{\rm
eff}\sum_{j=k+1}\left(e^{-i\beta_{kj}}m_{k}^{\dagger}m_{j}-e^{i\beta_{kj}}m_{k}m^{\dagger}_{j}\right).$
(17)

Note if $k=3$, then $j=1$, and

$g_{\rm eff}=g_{\rm
eff}(f)=\frac{2g^{2}}{\omega}\sum_{n=1}^{\infty}\frac{J^{2}_{n}(f)}{n}\sin\left(\frac{n\pi}{3}\right).$
(18)

Consequently, the effective Hamiltonian in Eq. (12) can be expressed by an
coefficient matrix,

$H_{\rm
eff}=[m^{\dagger}_{1},m^{\dagger}_{2},m^{\dagger}_{3}]\begin{bmatrix}0&G^{*}e^{-if}&Ge^{-if/2}\\\
Ge^{if}&0&G^{*}e^{if/2}\\\
G^{*}e^{if/2}&Ge^{-if/2}&0\end{bmatrix}\begin{bmatrix}m_{1}\\\ m_{2}\\\
m_{3}\end{bmatrix},$ (19)

where $G=gJ_{0}(f)+ig_{\rm eff}$ is the effective coupling strength. Using
$\tan\phi\equiv g_{\rm eff}/[gJ_{0}(f)]$, the Hamiltonian has a more compact
form,

$H_{\rm
eff}=|G|[m^{\dagger}_{1},m^{\dagger}_{2},m^{\dagger}_{3}]\begin{bmatrix}0&e^{i\phi_{12}}&e^{-i\phi_{31}}\\\
e^{-i\phi_{12}}&0&e^{i\phi_{23}}\\\
e^{i\phi_{31}}&e^{-i\phi_{23}}&0\end{bmatrix}\begin{bmatrix}m_{1}\\\ m_{2}\\\
m_{3}\end{bmatrix},$ (20)

where $\phi_{12}=-\phi-f$, $\phi_{23}=-\phi+f/2$, and $\phi_{31}=-\phi+f/2$.
The eigenvalues of the coefficient matrix satisfy $E^{3}-3E-2\cos\Phi=0$,
where $\Phi\equiv\phi_{12}+\phi_{23}+\phi_{31}=-3\phi$ is defined as the
closed-loop phase up to $2n\pi$ with $n$ an integer. The system dynamics is
then found to be relevant to the driving intensity and frequency rather than
solely determined by the ratio factor $f$. $\Phi$ is thus regarded as the
synthetic magnetic flux. It is gauge-invariant by noting that in our model the
three YIGs form a closed loop and the accumulated phase of the state-current
(either chiral or not) has to be single-valued when going around this loop.
$\Phi=0$ or $\pi$ corresponds to the absence of the Floquet engineering. Then
the evolution of the magnons is symmetrical when one of the magnons (say
mode-$1$) is prepared at the target state and the rest two are in the same
states. The quantum state propagates from mode-$1$ to mode-$2$ and mode-$3$
simultaneously and then back to mode-$1$. This pattern repeats itself with no
indication of any preferred circulation direction.

##  III Numerical simulation of Chiral currents

###  III.1 Chiral state transfer

Figure 1: (a) and (b): Diagrams of the chiral state transfers in the magnons
along the clockwise and anticlockwise circles, respectively. (c): Diagram of
the chiral entanglement transfer along the clockwise direction, where two
yellow nodes connected by a solid line represent the two entangled magnons and
the blue node represents a separable magnon.

With a Hamiltonian similar to Eq. (20), that is obtained by the tunneling-term
modulation Roushan1  _et al._ (2017) rather than the longitudinal Floquet
driving, the time-reversal symmetry broken was distinguished when the loop
phase $\Phi=\pm\pi/2$. This idea can be generalized to realize a chiral state-
transfer in our continuous-variable system for arbitrary target states.

Due to Eq. (19), $\Phi=\pi/2$ yields $G=ig_{\rm eff}$ under the condition of
$J_{0}(f)=0$. It renders that $f=\sqrt{3}\Delta/\omega=2.4048$. Consequently,
$J_{1}(f)=0.5192$, $J_{2}(f)=0.4318$, $J_{3}(f)=0.1990$, $J_{4}(f)=0.0647$,
$J_{5}(f)=0.0164$, $J_{6}(f)=0.0034$. The Hamiltonian becomes

$H_{\rm eff}=ig_{\rm
eff}[m^{\dagger}_{1},m^{\dagger}_{2},m^{\dagger}_{3}]\begin{bmatrix}0&-e^{-if}&e^{-if/2}\\\
e^{if}&0&-e^{if/2}\\\
-e^{if/2}&e^{-if/2}&0\end{bmatrix}\begin{bmatrix}m_{1}\\\ m_{2}\\\
m_{3}\end{bmatrix}.$ (21)

Note when $N\geq 4$, $f_{kj}$’s are not the same in magnitude. Consequently,
the zeroth-order time-reversal symmetry Hamiltonian $H_{0}$ in Eq. (11)
remains in any case. It is then hard to observe a perfect chiral state-
transfer. The canonical transformation in the Heisenberg picture induced by
Eq. (21) can be written as

$\begin{bmatrix}m_{1}(t)\\\ m_{2}(t)\\\
m_{3}(t)\end{bmatrix}=T(t)\begin{bmatrix}m_{1}(0)\\\ m_{2}(0)\\\
m_{3}(0)\end{bmatrix}$ (22)

with

$T(t)=\frac{1}{3}\begin{bmatrix}x(t)&e^{-if}y(t)&e^{-\frac{if}{2}}z(t)\\\
e^{if}z(t)&x(t)&e^{\frac{if}{2}}y(t)\\\
e^{\frac{if}{2}}y(t)&e^{-\frac{if}{2}}z(t)&x(t)\end{bmatrix},$ (23)

where

$\displaystyle x(t)$ $\displaystyle=1+2\cos\left(\sqrt{3}g_{\rm eff}t\right),$
(24) $\displaystyle y(t)$ $\displaystyle=1-2\cos\left(\sqrt{3}g_{\rm
eff}t+\frac{\pi}{3}\right),$ $\displaystyle z(t)$
$\displaystyle=1-2\cos\left(\sqrt{3}g_{\rm eff}t-\frac{\pi}{3}\right).$

With respect to the time-evolved operator, it is interesting to find that
$m_{3}(t)=m_{1}(0)$ and $m_{2}(t)=m_{3}(0)$ when
$t=(2\pi/3+2n\pi)/(\sqrt{3}g_{\rm eff})$; $m_{2}(t)=m_{1}(0)$ and
$m_{3}(t)=m_{2}(0)$ when $t=(4\pi/3+2n\pi)/(\sqrt{3}g_{\rm eff})$; and
$m_{j}(t)=m_{j}(0)$ when $t=2n\pi/(\sqrt{3}g_{\rm eff})$ with $n$ an integer.
This transfer is exactly described by the clockwise rotation $m_{1}\to
m_{3}\to m_{2}\to m_{1}$ in Fig. 1(a), corresponding to a chiral evolution in
the Schrödinger picture, i.e., $|\varphi_{1}\varphi_{2}\varphi_{3}\rangle\to$
$|\varphi_{2}\varphi_{3}\varphi_{1}\rangle\to$
$|\varphi_{3}\varphi_{1}\varphi_{2}\rangle$, where $\varphi_{j}$ is arbitrary
for the $j$th mode.

The state-transfer fidelity in the Schrödinger picture can be measured by the
time-dependent state population $P_{j}(t)$, $j=1,2,3$, which can be defined by

$P_{j}(t)=\sum_{C_{n}\neq
0}\left|\left\langle\varphi(t)\Big{|}\frac{(m^{\dagger}_{j})^{n}}{\sqrt{n!}}\Big{|}000\right\rangle\right|^{2}.$
(25)

It describes the overlap between the target state and the $j$th mode’s state
without regarding the influence from local dynamical phases Qi and Jing
(2022b). Here the initial state of the magnons is assumed to be
$|\varphi(0)\rangle=\sum_{n}C_{n}|n00\rangle\equiv\sum_{n}C_{n}|n\rangle_{1}|0\rangle_{2}|0\rangle_{3}$
with the normalized coefficients $C_{n}$, i.e., the first magnon is prepared
as an arbitrary target state and the other two are in their ground states. The
dynamics of the magnons is analytically determined by Eq. (22) or numerically
calculated by the Floquet-driving Hamiltonian in Eq. (7).

For $|\varphi(0)\rangle=|100\rangle$, we have

$|\varphi(t)\rangle=\frac{1}{3}\left[x(t)|100\rangle+e^{-if}z(t)|010\rangle+e^{\frac{-if}{2}}y(t)|001\rangle\right].$
(26)

While for an arbitrary superposed state

$|\varphi(0)\rangle=\sum_{n}C_{n}|n00\rangle=\sum_{n}\frac{C_{n}}{\sqrt{n!}}\left[m^{\dagger}_{1}(0)\right]^{n}|000\rangle,$
(27)

the time-evolved state can be written as

$|\varphi(t)\rangle=\sum_{n}\frac{C_{n}}{\sqrt{n!}}\left[\sum_{j}m^{\dagger}_{1}(0)T_{j1}^{\dagger}(t)\right]^{n}|000\rangle.$
(28)

Then due to Eqs. (22) and (23), we have
$\varphi(t_{3})=\sum_{n}C_{n}|00n\rangle$ at the desired time
$t_{3}=2\pi/(3\sqrt{3}g_{\rm eff})$ and
$\varphi(t_{2})=\sum_{n}C_{n}|0n0\rangle$ at $t_{2}=4\pi/(3\sqrt{3}g_{\rm
eff})$, up to certain local phases. By virtue of Eq. (25), we have
$P_{1}=\sum_{n}x^{2n}(t)|C_{n}|^{2}/9^{n}$,
$P_{2}=\sum_{n}z^{2n}(t)|C_{n}|^{2}/9^{n}$, and
$P_{3}=\sum_{n}y^{2n}(t)|C_{n}|^{2}/9^{n}$. During the state transfer of
$m_{1}\to m_{3}$, the population of $m_{2}$ achieve the side peak value at
$t=(\pi/3+2n\pi)/(\sqrt{3}g_{\rm eff})$, which is undesired for $m_{1}\to
m_{3}$ yet cannot be under control. In particular, $P_{2}^{\rm
max}=\sum_{n}|C_{n}|^{2}/9^{n}$ at those moments. Similarly, when
$t=(\pi+2n\pi)/(\sqrt{3}g_{\rm eff})$ and $t=(5\pi/3+2n\pi)/(\sqrt{3}g_{\rm
eff})$, $P_{1}$ and $P_{3}$ are expected to achieve their side peak values in
the middle of $m_{3}\to m_{2}$ and $m_{2}\to m_{1}$, respectively. In
addition, one can even transfer a mixed state by
$\rho(t)=\sum_{j}p_{j}|\varphi(t)\rangle_{jj}\langle\varphi(t)|$ based on Eq.
(28).

Figure 2: Time evolutions of the state populations $P_{j}(t)$, $j=1,2,3$, for
the three magnons by the numerical simulation with the Floquet-driving
Hamiltonian (7) and the analytical evaluation with the canonical
transformation (22). The initial state is
$|\varphi(0)\rangle=1/\sqrt{3}\sum_{n=1}^{3}|n00\rangle$. $\omega=20g$ and
$\phi_{j}=2\pi j/3$.

Figure 2 is used to verify the time evolutions of the state populations
$P_{j}(t)$ by the Floquet-driving Hamiltonian (7) under the chiral condition
of $J_{0}(f)=0$ in comparison with those by the effective Hamiltonian (21) or
the canonical transformation (22). It is found that the analytical results
(lines with markers) do match with their numerical counterparts (lines without
markers) for the normalized time evolution. The periodic dynamics is
determined once $g$ is fixed.

Figure 3: Clockwise state-transfer of the three magnons described by $P_{j}$’s
under the Floquet-driving Hamiltonian (7). The initial states for $m_{1}$ are
(a) the Fock state $|n=1\rangle$, (b) the coherent state with $\beta=1$, (c)
the cat state with $\zeta=1$, and (d) the thermal state with $\bar{n}=1$.
$\omega=20g$ and $\phi_{j}=2\pi j/3$.

In Fig. 3, the chiral transfer protocol is applied to Fock state, Glauber
coherent state $|\beta\rangle$, cat state
$(|\zeta\rangle+|-\zeta\rangle)/\sqrt{2+2e^{-2|\zeta|^{2}}}$ with
$|\zeta\rangle$ a coherent state, and thermal state
$\rho=\sum_{n}p_{n}|n\rangle\langle n|$ with the Fock-state occupation
$p_{n}=(\bar{n})^{n}/(1+\bar{n})^{n+1}$ and the average excitation number
$\bar{n}$. Periodically all the state populations approach unit in sequence
indicating a perfect clockwise transfer amongst $m_{1}$, $m_{2}$, and $m_{3}$.
It is interesting to find that the side peak for the Fock state is more
significant than the other states that occupy the vacuum state. It is caused
by the distinction between the side peak and initial values $P_{2}^{\rm
max}-P_{2}(0)=\sum_{n}(1/9)^{n}|C_{n}|^{2}-|C_{0}|^{2}$. In Fig. 3(a), (b) and
(c), the distinctions are about $0.11$, $0.05$ and $0.01$, respectively.
Figures 2 and 3 demonstrate that the Floquet-driving Hamiltonian in Eq. (21)
under the chiral conditions $J_{0}(f)=0$ and $\phi_{j}=2\pi j/3$ would ensure
a high-fidelity chiral state-transfer for arbitrary target states. And the
transfer period is irrespective to the chosen initial state.

Figure 4: Anticlockwise state-transfer of the magnons described by $P_{j}$’s
under the Floquet-driving Hamiltonian (7). The initial states of the magnon-
mode $m_{1}$ are (a) the Fock state $|n=1\rangle$, (b) the coherent state with
$\beta=1$, (c) the cat state with $\zeta=1$, and (d) the thermal state with
$\bar{n}=1$. $\omega=20g$, $\phi_{1}=2\pi/3$, $\phi_{2}=2\pi$, and
$\phi_{3}=4\pi/3$.

By setting $\phi_{1}=2\pi/3$, $\phi_{2}=2\pi$, and $\phi_{3}=4\pi/3$, which is
an alternative uniform distribution in the range of $[0,2\pi]$, the
coefficient matrix in Eq. (19) becomes

$\begin{bmatrix}0&Ge^{-if/2}&G^{*}e^{-if}\\\ G^{*}e^{if/2}&0&G^{*}e^{if/2}\\\
Ge^{if}&Ge^{if/2}&0\end{bmatrix}$ (29)

where $G=gJ_{0}(f)-ig_{\rm eff}$. That induces the synthetic magnetic flux
with $\Phi=-\pi/2$. It amounts to exchange the magnon mode-$2$ and $3$, giving
rise to the anticlockwise rotation shown in Fig 1(b). Then the effective
Hamiltonian in Eq. (21) is modified to

$H_{\rm eff}=ig_{\rm
eff}[m^{\dagger}_{1},m^{\dagger}_{2},m^{\dagger}_{3}]\begin{bmatrix}0&e^{-if/2}&-e^{-if}\\\
-e^{if/2}&0&-e^{-if/2}\\\
e^{if}&e^{if/2}&0\end{bmatrix}\begin{bmatrix}m_{1}\\\ m_{2}\\\
m_{3}\end{bmatrix}$ (30)

Through a canonical transformation similar to Eq. (22), one can easily analyse
the anticlockwise state-transfer. For example, the time-dependent state
$|\varphi(t)\rangle$ in Eq. (26) for the single-excitation state becomes

$|\varphi(t)\rangle=\frac{1}{3}\left[x(t)|100\rangle+e^{\frac{-if}{2}}y(t)|010\rangle+e^{-if}z(t)|001\rangle\right],$
(31)

where $x(t)$, $y(t)$, and $z(t)$ are given by Eq. (24). The numerical results
for various states are provided in Fig. 4, exhibiting the opposite chirality
to Fig. 3. During the anticlockwise state-transfer $m_{1}\to m_{2}\to m_{3}$,
all the magnon modes achieve a nearly unit fidelity in a periodical sequence.
Anticlockwise chiral state-transfer also applies to arbitrary states.

Figure 5: Clockwise transfer of the pairwise concurrence in the three-magnon
system under the Floquet-driving Hamiltonian (7). Magnonic mode-$1$ and
mode-$3$ are initially set as (a) the Bell state
$(|00\rangle+|11\rangle)/\sqrt{2}$ and (b) the NOON state
$(|0N\rangle+|N0\rangle)/\sqrt{2}$ with $N=5$. The other parameters are the
same as Fig. 2.

If any two magnons in our model could be prepared as an entangled state and
the third one is separably initialized as the ground state, we can demonstrate
a high-fidelity chiral transfer for the pairwise quantum entanglement [see
Fig. 1(c)]. In particular, we consider that mode-$1$ and mode-$3$ are prepared
as a double-excitation Bell state $(|00\rangle+|11\rangle)/\sqrt{2}$ or a NOON
state $(|0N\rangle+|N0\rangle)/\sqrt{2}$ under the clockwise-transfer
condition, i.e., $\Phi=\pi/2$.

It is well known that the concurrence could be used as a sufficient and
necessary entanglement criterion for any bipartite system in pure state and
any bipartite two-state system. According to its definition Wootters (1998)
for a two-qubit system $AB$,
$C=\max\\{0,\lambda_{1}-\lambda_{2}-\lambda_{3}-\lambda_{4}\\}$, where
$\lambda_{j}$ are the square root of the eigenvalues of
$\rho_{AB}(\sigma^{A}_{y}\otimes\sigma^{B}_{y})\rho^{*}_{AB}(\sigma^{A}_{y}\otimes\sigma^{B}_{y})$
in the decreasing order. For a pure state
$a|11\rangle+b|10\rangle+c|01\rangle+d|00\rangle$ with normalized
coefficients, its concurrence can be easily derived as $C=2|bc-ad|$.

Inspired by the expression for the pure state, we here employ a modified
concurrence to measure the chiral-transfer of the entangled state in our
model. If the initial state is the Bell state
$(|00\rangle+|11\rangle)/\sqrt{2}$, then the concurrence between mode-$j$ and
mode-$k$ is defined as
$C_{jk}(t)=2|\langle\varphi(t)|00\rangle_{jk}||\langle\varphi(t)|11\rangle_{jk}|$;
and if the target state is initially chosen as the NOON state
$(|0N\rangle+|N0\rangle)/\sqrt{2}$, then
$C_{jk}(t)=2|\langle\varphi(t)|0N\rangle_{jk}||\langle\varphi(t)|N0\rangle_{jk}|$.
These definitions amount to the original occurrence under the perfect state
transfer and become sufficient but not necessary for the entangled states
during the time evolution, which is calculated with the Floquet-driving
Hamiltonian (7).

In Fig. 5, the blue solid lines, the red dot-dashed lines, and the yellow
dotted lines represent the concurrences for mode-$1$ and mode-$3$, mode-$3$
and mode-$2$, and mode-$2$ and mode-$1$, respectively. It is found that both
the double-excitation Bell state and the NOON state could be transferred in a
chiral way along the magnon triangle. At the same desired moments for the
faithful state transfer, one can observe a nearly perfect entanglement
transfer. Comparing Fig. 5(a) and Fig. 5(b), it is found that the undesired
magnon pair obtains a lower probability of entanglement generation in
transferring NOON state than Bell state. The side peaks of the Bell state
$(|00\rangle+|11\rangle)/\sqrt{2}$ are dominant because its constituent
$|00\rangle$ does not participate in the time evolution. Our model then widens
the application range of the chiral state-transfer.

###  III.2 Chiral current

Chiral current is an effective measure for the symmetry of the system
Hamiltonian in condensed matter systems, which can be probed by measuring the
energy spectrum under varying synthetic field $\Phi$ Roushan1  _et al._
(2017). Using the definition in Ref. Roushan1  _et al._ (2017), we can
consider the continuity equation for each magnon mode, e.g.,
$\partial_{t}\hat{n}_{k}=I_{\rm in}-I_{\rm out}$, where $\hat{n}_{k}\equiv
m^{\dagger}_{k}m_{k}$ is the occupation operator of the $k$th magnon mode and
$I_{\rm in}$ ($I_{\rm out}$) represents the input (output) current with
respect to mode-$k$. Under the necessary condition for the clockwise transfer
in Fig. 1(a), i.e., $\phi_{j}=2\pi j/3$, $j=1,2,3$, we have three equations of
motion:

$\displaystyle\partial_{t}\hat{n}_{1}$
$\displaystyle=I_{12}-I_{31}=i\left[H_{\rm eff},\hat{n}_{1}\right],$ (32)
$\displaystyle\partial_{t}\hat{n}_{2}$
$\displaystyle=I_{23}-I_{12}=i\left[H_{\rm eff},\hat{n}_{2}\right],$
$\displaystyle\partial_{t}\hat{n}_{3}$
$\displaystyle=I_{31}-I_{23}=i\left[H_{\rm eff},\hat{n}_{3}\right],$

where $H_{\rm eff}$ is the effective Hamiltonian in Eq. (19) or Eq. (20) and
$I_{jk}$ represents the current from mode-$j$ to mode-$k$. Then we have

$\displaystyle
I_{12}=iGe^{if}m_{1}m^{\dagger}_{2}-iG^{*}e^{-if}m^{\dagger}_{1}m_{2},$ (33)
$\displaystyle
I_{31}=-iG^{*}e^{-if/2}m_{1}m^{\dagger}_{3}+iGe^{if/2}m^{\dagger}_{1}m_{3},$
$\displaystyle
I_{23}=iGe^{-if/2}m_{2}m^{\dagger}_{3}-iG^{*}e^{if/2}m^{\dagger}_{2}m_{3}.$

Normalized by the coupling-strength $|G|$ of $H_{\rm eff}$, the effective
current operator from mode-$j$ to mode-$k$ can be defined as

$\hat{I}_{kj}=ie^{-i\phi_{kj}}m_{k}m^{\dagger}_{j}-ie^{i\phi_{kj}}m^{\dagger}_{k}m_{j}.$
(34)

Accordingly, the circle current operator is defined as
$\hat{I}=\hat{I}_{12}+\hat{I}_{23}+\hat{I}_{31}$. It is straightforward to
find that the circle current operator is an integral of motion because it
commutates with the effective Hamiltonian, i.e., $[H_{\rm eff},\hat{I}]=0$. It
means that the expectation value of $\langle\hat{I}\rangle$ can not be used to
measure the time-reversal symmetry of our model.

We have known that the magnon system has a symmetrical evolution for a
synthetic magnetic flux $\Phi=0$ or $\Phi=\pi$ under the effective Hamiltonian
in Eq. (20). In this case, arbitrary state of mode-$1$ simultaneously
propagates to mode-$2$ and mode-$3$ and periodically goes back to mode-$1$.
Consequently, the current $\langle\hat{I}_{23}\rangle$ from mode-$2$ to
mode-$3$ vanishes. A nonzero current $\langle\hat{I}_{23}\rangle$ appears when
the time-reversal symmetry is broken when $\Phi\neq 0,\pi$ and becomes maximum
when $\Phi=\pm\pi/2$. We check the condition of $\Phi=\pi/2$ under
$\phi_{j}=2\pi j/3$. Using $|\varphi(t)\rangle$ in Eq. (26) with the initial
Fock state $|100\rangle$, it is found that

$\langle\hat{I}_{23}(t)\rangle=\frac{4}{9}\left[2\cos^{2}\left(\sqrt{3}g_{\rm
eff}t\right)-\cos\left(\sqrt{3}g_{\rm eff}t\right)-1\right].$ (35)

The maximal value is achieved as $\langle\hat{I}_{23}^{\rm max}\rangle=8/9$
when $\sqrt{3}g_{\rm eff}t=\pi+2n\pi$. Given an arbitrary state in Eq. (28),
we have

$\langle\hat{I}_{23}^{\rm
max}\rangle=\frac{8}{9}\sum_{n}nC^{2}_{n}=\frac{8}{9}\bar{n},$ (36)

where $\bar{n}$ is the average excitation number of the state. For a coherent
state $|\beta\rangle$, $\bar{n}=|\beta|^{2}$ and for a cat state
$(|\zeta\rangle+|-\zeta\rangle)/\sqrt{2+2e^{-2|\zeta|^{2}}}$,
$\bar{n}=|\zeta|^{2}[1-\exp(-2|\zeta|^{2})]/[1+\exp(-2|\zeta|^{2})]$. Then the
expectation value of $\langle\hat{I}_{23}(t)\rangle$ is effective to measure
the time-reversal or parity-reversal symmetry of the system Hamiltonian.

Figure 6: (a) Maximal value of the current from $m_{2}$ to $m_{3}$ as a
function of the closed-loop phase $\Phi$ under the effective Hamiltonian (20).
The initial states for $m_{1}$ are the Fock state $|n=1\rangle$ (the blue
solid line), the coherent state with $\beta=1$ (the red dot-dashed line), and
the cat state with $\zeta=1$ (the yellow dotted line). (b) Maximal value of
the current from $m_{2}$ to $m_{3}$ for various Fock states as a function of
$\Phi$.

The preceding analysis is consistent with the numerical simulation in Fig. 6
about $\langle\hat{I}_{23}^{\rm max}\rangle$ as a function of the close-loop
phase or the synthetic magnetic field $\Phi$. It is found that the maximal
current vanishes when $\Phi=0$ and $\Phi=\pi$. It can be understood that at
these points, the system Hamiltonian becomes real that cannot break the time-
reversal symmetry. Otherwise, any finite $\Phi$ breaks the symmetry and leads
to the chiral currents. It is interesting to find in Fig. 6(a) that the
maximal currents for the Fock state $|n=1\rangle$ and the coherent state with
$\beta=1$ are equivalent to each other. Both of them equal to $0.89\approx
8/9$ at the peak point $\Phi=\pi/2$. They are larger than the cat state with
$\zeta=1$, whose peak value is $0.68$. A larger average excitation number
yields a larger $\langle\hat{I}_{23}^{\rm max}\rangle$ with the same $\Phi$.
This result is also supported by Fig. 6(b), in which a monotonic pattern
appears with increasing $n$.

##  IV Fidelity analysis

###  IV.1 Fidelity analysis under magnon damping

Taking account the magnon damping into consideration, we can study the chiral
state-transfer fidelity in a standard open-quantum-system framework. In
general, the full Hamiltonian including the environments can be written as

$H_{\rm
tot}=H(t)+\sum_{k=1}^{3}\left(m_{k}B_{k}^{\dagger}+m_{k}^{\dagger}B_{k}\right),$
(37)

where $H(t)$ is the Floquet-driving Hamiltonian (7) with $N=3$ and
$B_{k}=\sum_{j}g_{kj}a^{\dagger}_{kj}e^{-i\omega_{kj}t}$ is a collective
environmental operator. $\omega_{kj}$ is the eigenfrequency for the $j$th
field mode in the $k$th environment and $g_{kj}$ is its coupling strength with
the $k$th magnon mode. Under the Markovian approximation and tracing out the
degrees of freedom of the external environments (assumed to be at the vacuum
state), we arrive at the master equation Breuer and Petruccione (2022) for the
density-matrix operator $\rho(t)$ for the magnon system:

$\dot{\rho}(t)=-i\left[H(t),\rho(t)\right]+\kappa\sum_{k=1}^{3}\mathcal{L}[m_{k}]\rho(t),\\\
$ (38)

where $\kappa$ is the damping rate. The dissipative superoperator
$\mathcal{L}$ is in a Lindblad form,

$\mathcal{L}[o]\rho\equiv\frac{1}{2}\left(2o\rho
o^{\dagger}-o^{\dagger}o\rho-\rho o^{\dagger}o\right),$ (39)

where $o=m_{1},m_{2},m_{3}$ are the decay operators.

Figure 7: Clockwise state-transfer of the three magnons described by $P_{j}$’s
under various damping rates: (a) $\kappa/g=0.001$ and (b) $\kappa/g=0.01$. The
initial state is $|n=1\rangle|0\rangle|0\rangle$. $\omega=20g$ and
$\phi_{j}=2\pi j/3$.

In Fig. 7(a), one can observe that the open-system dynamics under
$\kappa/g=0.001$ is very close to the closed-system result in Fig. 3(a). And
in Fig. 7(b), it is found that at least one round of chiral transfer presents
when $\kappa/g=0.01$. It is thus required that the effective coupling strength
$g$ is about two orders higher than the damping rate in magnitude. For YIG
spheres of a typical size about $250$ $\mu$m in recent experiments Tabuchi
_et al._ (2014, 2015); Lachance-Quirion  _et al._ (2020), the coupling
strength $g$ was found to be in the order of $1\sim 10$ MHz and the damping
rate $\kappa\sim 1$ MHz. Then our protocol could become practical in a near
future by increasing the strength of the magnetic dipole interaction to the
order of $100$ MHz Lachance-Quirion  _et al._ (2019).

Alternatively, we can apply the external quantum control, such as the
dynamical decoupling method Viola and Lloyd (1998), to our system to
neutralize the magnon damping effects. The system coherence can then be saved
by the $\pi$-pulses exerted on the magnons. It is known that

$e^{i\int_{\tau_{n}}^{\tau_{n+1}}dt\Omega(t)m_{k}^{\dagger}m_{k}}m_{k}e^{-i\int_{\tau_{n}}^{\tau_{n+1}}dt\Omega(t)m_{k}^{\dagger}m_{k}}=-m_{k}$
(40)

as long as $\int_{\tau_{n}}^{\tau_{n+1}}dt\Omega(t)=\pi$. In the ideal
situation for dynamical decoupling,
$\Omega(t)=\pi/(\tau_{n+1}-\tau_{n})\delta(t-s)$ with
$s\in(\tau_{n},\tau_{n+1})$. Thus it could reverse the decoherence dynamics
induced by the unwanted terms $m_{k}B_{k}^{\dagger}+H.c.$ in Eq. (37), but
does not change the effective Hamiltonian in our protocol.

###  IV.2 Fidelity under systematic errors

Our protocol based on Eq. (1) has been discussed under the ideal condition
that all the coupling strengths $g_{am}$ between photon and magnon modes are
the same in magnitude. While in practice Zhang  _et al._ (2014); Lachance-
Quirion  _et al._ (2019), they are associated with the individual locations of
the YIG spheres in the cavity. We first consider the systematic errors raised
by the nonequal coupling strengths. In particular, it is assumed that
$g_{am}^{(1)}=g_{am}$, $g_{am}^{(2)}=g_{am}(1+\delta)$, and
$g_{am}^{(3)}=g_{am}(1-\delta)$, where the superscript $j\in\\{1,2,3\\}$ marks
the coupling strength between the common photon and the $j$th magnon and
$\delta$ is used to estimate the magnitude of the relative error. With $N=3$,
the system Hamiltonian in Eq. (1) is then rewritten as

$\displaystyle H$
$\displaystyle=\omega_{a}a^{\dagger}a+\omega_{m}\sum_{k=1}^{3}m^{\dagger}_{k}m_{k}+g_{am}\left(am^{\dagger}_{1}+a^{\dagger}m_{1}\right)$
(41)
$\displaystyle+g_{am}(1+\delta)\left(am^{\dagger}_{2}+a^{\dagger}m_{2}\right)$
$\displaystyle+g_{am}(1-\delta)\left(am^{\dagger}_{3}+a^{\dagger}m_{3}\right).$

Under the conditions that $\phi_{j}=2\pi j/3$ and $J_{0}(f)=0$ for the
clockwise chiral transfer, the effective Hamiltonian in Eq. (21) is modified
by changing the coefficient matrix to be

$ig_{\rm eff}\begin{bmatrix}0&-(1-\delta^{2})e^{if}&(1+\delta)e^{if/2}\\\
(1-\delta^{2})e^{-if}&0&-(1-\delta)e^{-if/2}\\\
-(1+\delta)e^{-if/2}&(1-\delta)e^{if/2}&0\end{bmatrix}.$ (42)

Figure 8: (a) State population of the $3$rd magnon at the desired moment
$t_{3}=2\pi/(3\sqrt{3}g_{\rm eff})$ and (b) State population of the $2$nd
magnon at the desired moment $t_{2}=4\pi/(3\sqrt{3}g_{\rm eff})$ as functions
of the magnitude of the relative error $\delta$ under the nonideal Hamiltonian
(42). The initial states for $m_{1}$ are the Fock state $|n=1\rangle$ (the
blue solid line), the coherent state with $\beta=1$ (the red dot-dashed line),
and the cat state with $\zeta=1$ (the yellow dotted line).

In Fig. 8, we present the sensitivity of the state populations [see Eq. (25)]
to the error $\delta$ for $g_{ma}$ at the desired moments for the chiral
state-transfer obtained by Eq. (24). It is found that the robustness of our
protocol against the systematic errors becomes weaker with a larger average
population. The cat state with $\zeta=1$ is found to be more robust in
comparison with the Fock state $|n=1\rangle$ and the coherent state with
$\beta=1$. For all of them, both transferred populations of mode-$2$ and
mode-$3$ could be maintained above $0.98$ in the presence of about $10\%$
errors or fluctuations in the photon-magnon coupling strength.

Another systematic error in our model is associated with the nonlinearity in
the expansion of the YIG crystal or the dipolar anisotropy energy quadratic in
the magnon numbers, which is parameterized with the Kerr coefficient $K$. When
the excitation number in the spin wave is much smaller than the total number
of spins in YIG sphere, the system Hamiltonian in Eq. (6) can be modified to
Wang  _et al._ (2018, 2016b)

$\tilde{H}_{\rm
Kerr}=g\sum_{k<j}\left(m_{k}m^{\dagger}_{j}+m^{\dagger}_{k}m_{j}\right)+K\sum_{k}\left(m^{\dagger}_{k}m_{k}\right)^{2}.$
(43)

The Kerr coefficient $K/2\pi\approx 6.4\times 10^{-9}$ Hz, which is much
smaller than the effective coupling strength $g$. Given the timescale we
considered in the preceding chiral dynamics, we can safely apply our model in
such a hybrid photon-magnon system when the number of the excited spins
$\langle m^{\dagger}_{k}m_{k}\rangle$ is less than $10^{6}$.

##  V Conclusion

In summary, we have proposed a chiral state-transfer proposal in a Floquet
cavity-magnonic system, which consists of a microwave cavity coupled to three
YIG spheres in their Kittle mode. The three magnon modes constitute a triangle
loop of mutual interaction through eliminating the photon mode in a state-
resolved way. The Floquet driving with desired parameters forms complex-valued
hopping of magnons, that generates a tunable synthetic magnetic flux. Our
results open a path toward simulation of a time-reversal-broken continuous-
variable system, yielding the chiral transfer for arbitrary states and even
quantum entanglement. It is beyond the limitation of fixed excitation numbers.
We evaluate the chiral current to measure the symmetry of the system
Hamiltonian. We discuss the feasibility of the chiral transfer in the presence
of magnon damping. The stability of our proposal is also tested with respect
to the systematic errors in the coupling strength and the nonlinear terms in
the system. Our work in pursuit of the quantum chiral state transfer provides
an interesting application in cavity magnonics as a promising hybrid platform
for control over the emerging phases.

Our theory can be straightforwardly generalized to arbitrary continuous-
variable systems that could be indirectly coupled to each other by a common
mediator. For example, one can replace our magnon modes with microwave photon
modes in superconducting circuit with a much reduced ratio of decoherence and
coupling strength Niemczyk  _et al._ (2010). In addition, although our
protocol has no direct relation with topological variables, the anomalous edge
states, and the bulk-edge correspondence in the momentum space, we expect that
the edge state Cai and Wang (2021); Roushan1  _et al._ (2017); Rudner  _et
al._ (2013) can arise from the cavity-magnonic system of a larger size under
Floquet engineering. It will bring extra interesting physics.

## Acknowledgments

We acknowledge financial support from the National Science Foundation of China
(Grants No. 11974311 and No. U1801661).

##  Appendix A RWA in Hamiltonian (1)

This appendix contributes to verifying that it is valid to apply the rotating-
wave approximation in the model Hamiltonian (1). Without loss of generality,
it is instructive to focus on the state exchange (Rabi oscillation) between
two magnons mediated by a common photon mode as what we concerned in Eq. (2).
Including the counter-rotating terms, the Hamiltonian reads

$H=\omega_{a}a^{\dagger}a+\omega_{m}\sum_{k=1}^{2}m^{\dagger}_{k}m_{k}+g_{am}\left(a+a^{\dagger}\right)\sum_{k=1}^{2}\left(m^{\dagger}_{k}+m_{k}\right).$
(44)

And it becomes

$H_{\rm
rwa}=\omega_{a}a^{\dagger}a+\omega_{m}\sum_{k=1}^{2}m^{\dagger}_{k}m_{k}+g_{am}\sum_{k=1}^{2}\left(am^{\dagger}_{k}+a^{\dagger}m_{k}\right)$
(45)

under RWA.

Figure 9: Time evolution of the state populations $P(t)$ for the two magnon
modes evaluated by the Hamiltonian in Eq. (44) and that under RWA in Eq. (45).
The initial state is $|\varphi(0)\rangle=|1\rangle_{m_{1}}|0\rangle_{m_{2}}$
and the detuning is chosen as $\Delta=10g_{am}$. In (a),
$\omega_{m}=100g_{am}$ and in (b), $\omega_{m}=200g_{am}$.

To demonstrate the effect from the counter-rotating terms in a practical
cavity-magnon system Wang  _et al._ (2018, 2016b) on Rabi oscillation of the
two modes, we plot their population dynamics using the total Hamiltonian in
Eq. (44) and that under RWA in Eq. (45). One can find that in Fig. 9(a), the
difference between $H$ and $H_{\rm rwa}$ is not significant when
$\omega_{m}=100g_{am}$; and in Fig. 9(b) it becomes ignorable when
$\omega_{m}=200g_{am}$. Thus we are sure that it is not necessary to include
the counter-rotating terms since $g_{am}\ll\Delta=\omega_{a}-\omega_{m}$ is in
the strong instead of ultrastrong regimes.

Note the counter-rotating terms does not constitute an obstacle for our state-
resolved method Qi and Jing (2021b) to eliminate the mediator [see from Eq.
(2) to Eq. (5)]. We can still obtain an effective Hamiltonian for the magnons
as in Eq. (5), and it is little different from the current result in the
dispersive regime of $g_{am}$.

## References

  * Rameshti  _et al._ (2021) B. Z. Rameshti, S. V. Kusminskiy, J. A. Haigh, K. Usami, D. Lachance-Quirion, Y. Nakamura, C.-M. Hu, H. X. Tang, G. E. Bauer, and Y. M. Blanter,  _Cavity magnonics,_ arxiv: 2106. 09312v1 (2021).
  * Lachance-Quirion  _et al._ (2019) D. Lachance-Quirion, Y. Tabuchi, A. Gloppe, K. Usami, and Y. Nakamura,  _Hybrid quantum systems based on magnonics,_Appl. Phys. Express 12, 070101 (2019).
  * Li  _et al._ (2020) Y. Li, W. Zhang, V. Tyberkevych, W. K. Kwok, and V. Novosad,  _Hybrid magnonics: physics, circuits, and applications for coherent information processing,_J. Appl. Phys. 128, 130902 (2020).
  * Yuan  _et al._ (2022) H. Yuan, Y. Cao, A. Kamra, R. A. Duine, and P. Yan,  _Quantum magnonics: When magnon spintronics meets quantum information science,_Phys. Rep. 965, 1 (2022).
  * Soykal and Flatté (2010a) O. O. Soykal and M. E. Flatté,  _Strong field interactions between a nanomagnet and a photonic cavity,_Phys. Rev. Lett. 104, 077202 (2010a).
  * Soykal and Flatté (2010b) O. O. Soykal and M. E. Flatté,  _Size dependence of strong coupling between nanomagnets and photonic cavities,_Phys. Rev. B 82, 104413 (2010b).
  * Tabuchi  _et al._ (2014) Y. Tabuchi, S. Ishino, T. Ishikawa, R. Yamazaki, K. Usami, and Y. Nakamura,  _Hybridizing ferromagnetic magnons and microwave photons in the quantum limit,_Phys. Rev. Lett. 113, 083603 (2014).
  * Zhang  _et al._ (2014) X. Zhang, C.-L. Zou, L. Jiang, and H. X. Tang,  _Strongly coupled magnons and cavity microwave photons,_Phys. Rev. Lett. 113, 156401 (2014).
  * Wang  _et al._ (2018) Y.-P. Wang, G.-Q. Zhang, D. Zhang, T.-F. Li, C.-M. Hu, and J. Q. You,  _Bistability of cavity magnon polaritons,_Phys. Rev. Lett. 120, 057202 (2018).
  * Shen  _et al._ (2021) R.-C. Shen, Y.-P. Wang, J. Li, S.-Y. Zhu, G. S. Agarwal, and J. Q. You,  _Long-time memory and ternary logic gate using a multistable cavity magnonic system,_Phys. Rev. Lett. 127, 183202 (2021).
  * Lachance-Quirion  _et al._ (2020) D. Lachance-Quirion, S. Piotr Wolski, Y. Tabuchi, S. Kono, K. Usami, and Y. Nakamura,  _Entanglement-based single-shot detection of a single magnon with a superconducting qubit,_Science 367, 425 (2020).
  * Tabuchi  _et al._ (2015) Y. Tabuchi, S. Ishino, A. Noguchi, T. Ishikawa, R. Yamazaki, K. Usami, and Y. Nakamura,  _Coherent coupling between a ferromagnetic magnon and a superconducting qubit,_Science 349, 405 (2015).
  * Qi and Jing (2022a) S.-f. Qi and J. Jing,  _Generation of bell and greenberger-horne-zeilinger states from a hybrid qubit-photon-magnon system,_Phys. Rev. A 105, 022624 (2022a).
  * Yan and Jing (2021) J.-s. Yan and J. Jing,  _External-level assisted cooling by measurement,_Phys. Rev. A 104, 063105 (2021).
  * Neuman  _et al._ (2020) T. c. v. Neuman, D. S. Wang, and P. Narang,  _Nanomagnonic cavities for strong spin-magnon coupling and magnon-mediated spin-spin interactions,_Phys. Rev. Lett. 125, 247702 (2020).
  * Qi and Jing (2021a) S.-f. Qi and J. Jing,  _Magnon-mediated quantum battery under systematic errors,_Phys. Rev. A 104, 032606 (2021a).
  * Zhang  _et al._ (2016) X. Zhang, C.-L. Zou, L. Jiang, and H. Tang,  _Cavity magnonmechanics,_Sci. Adv. 2, e1501286 (2016).
  * Li  _et al._ (2018a) J. Li, S.-Y. Zhu, and G. S. Agarwal,  _Magnon-photon-phonon entanglement in cavity magnomechanics,_Phys. Rev. Lett. 121, 203601 (2018a).
  * Qi and Jing (2021b) S.-f. Qi and J. Jing,  _Magnon-assisted photon-phonon conversion in the presence of structured environments,_Phys. Rev. A 103, 043704 (2021b).
  * Li  _et al._ (2021) J. Li, Y.-P. Wang, W.-J. Wu, S.-Y. Zhu, and J. You,  _Quantum network with magnonic and mechanical nodes,_PRX Quant. 2, 040344 (2021).
  * Hatanaka  _et al._ (2022) D. Hatanaka, M. Asano, H. Okamoto, Y. Kunihashi, H. Sanada, and H. Yamaguchi,  _On-chip coherent transduction between magnons and acoustic phonons in cavity magnomechanics,_Phys. Rev. Appl. 17, 034024 (2022).
  * Potts  _et al._ (2021) C. A. Potts, E. Varga, V. A. S. V. Bittencourt, S. V. Kusminskiy, and J. P. Davis,  _Dynamical backaction magnomechanics,_Phys. Rev. X 11, 031053 (2021).
  * Shen  _et al._ (2022) R.-C. Shen, J. Li, Z.-y. Fan, Y.-P. Wang, and J. You,  _Mechanical bistability in kerr-modeified cavity magnomechanics,_ arxiv: 2206. 14588 (2022).
  * Ladd  _et al._ (2010) T. D. Ladd, F. Jelezko, R. Laflamme, Y. Nakamura, C. Monroe, and J. L. O’Brien,  _Quantum computers,_Nature (London) 464, 45 (2010).
  * Reiserer and Rempe (2015) A. Reiserer and G. Rempe,  _Cavity-based quantum networks with single atoms and optical photons,_Rev. Mod. Phys. 87, 1379 (2015).
  * Degen  _et al._ (2017) C. L. Degen, F. Reinhard, and P. Cappellaro,  _Quantum sensing,_Rev. Mod. Phys. 89, 035002 (2017).
  * Xu  _et al._ (2020) J. Xu, C. Zhong, X. Han, D. Jin, L. Jiang, and X. Zhang,  _Floquet cavity electromagnonics,_Phys. Rev. Lett. 125, 237201 (2020).
  * Lodahl  _et al._ (2017) P. Lodahl, S. Mahmoodian, S. Stobbe, A. Rauschenbeutel, P. Schneeweiss, J. Volz, H. Pichler, and P. Zoller,  _Chiral quantum optics,_Nature 541, 473 (2017).
  * Zhu  _et al._ (2018) H. Zhu, J. Yi, M.-Y. Li, J. Xiao, L. Zhang, C.-W. Yang, R. A. Kaindl, L.-J. Li, Y. Wang, and X. Zhang,  _Observation of chiral phonons,_Science 359, 579 (2018).
  * Sørensen  _et al._ (2005) A. S. Sørensen, E. Demler, and M. D. Lukin,  _Fractional quantum hall states of atoms in optical lattices,_Phys. Rev. Lett. 94, 086803 (2005).
  * Liu  _et al._ (2020) W. Liu, W. Feng, W. Ren, D.-W. Wang, and H. Wang,  _Synthesizing three-body interaction of spin chirality with superconducting qubits,_Appl. Phys. Lett. 116, 114001 (2020).
  * Cai and Wang (2021) H. Cai and D.-W. Wang,  _Topological phases of quantized light,_Nat. Sci. Rev. 8, nwaa196 (2021).
  * Shirley (1965) J. H. Shirley,  _Solution of the schrödinger equation with a hamiltonian periodic in time,_Phys. Rev. 138, B979 (1965).
  * Goldman and Dalibard (2014) N. Goldman and J. Dalibard,  _Periodically driven quantum systems: Effective hamiltonians and engineered gauge fields,_Phys. Rev. X 4, 031027 (2014).
  * Bukov  _et al._ (2015) M. Bukov, L. D’Alessio, and A. Polkovnikov,  _Universal high-frequency behavior of periodically driven system: from dynamical stabilization to floquet engineering,_Adv. in Phys. 64, 139 (2015).
  * Petiziol  _et al._ (2021) F. Petiziol, M. Sameti, S. Carretta, S. Wimberger, and F. Mintert,  _Quantum simulation of three-body interactions in weakly driven quantum systems,_Phys. Rev. Lett. 126, 250504 (2021).
  * Shao  _et al._ (2017) W. Shao, C. Wu, and X.-L. Feng,  _Generalized james’ effective hamiltonian method,_Phys. Rev. A 95, 032124 (2017).
  * Wang  _et al._ (2016a) D.-W. Wang, H. Cai, R.-B. Liu, and M. O. Scully,  _Mesoscopic superposition states generated by synthetic spin-orbit interaction in fock-state lattices,_Phys. Rev. Lett. 116, 220502 (2016a).
  * Li  _et al._ (2019) H. Li, H. Cai, J. Xu, V. V. Yakovlev, Y. Yang, and D.-W. Wang,  _Quantum photonic transistor controlled by an atom in a floquet cavity-qed system,_Opt. Express 27, 6946 (2019).
  * Wu  _et al._ (2018) Y. Wu, L.-P. Yang, M. Gong, Y. Zheng, H. Deng, Z. Yan, Y. Zhao, K. Huang, A. D. Castellano, W. J. Munro, K. Nemoto, D.-N. Zheng, C. Sun, Y.-x. Liu, X. Zhu, and L. Lu,  _An efficient and compact switch for quantum circuits,_npj. Quantum Inf. 4, 50 (2018).
  * Roushan1  _et al._ (2017) P. Roushan1, C. Neill, A. Megrant, Y. Chen, R. Babbush, R. Barends, B. Campbell, Z. Chen, B. Chiaro, A. Dunsworth, A. Fowler, E. Jeffrey, J. Kelly, E. Lucero, J. Mutus, P. ÓMalley, M. Neeley, C. Quintana, D. Sank, A. Vainsencher, J. Wenner, T. White, E. Kapit, H. Neven, and J. Martinis,  _Chiral ground-state currents of interacting photons in a synthetic magnetic field,_Nat. Phys. 13, 146 (2017).
  * Kyriienko and Sørensen (2018) O. Kyriienko and A. S. Sørensen,  _Floquet quantum simulation with superconducting qubits,_Phys. Rev. Appl. 9, 064029 (2018).
  * Wang  _et al._ (2019) D.-W. Wang, C. Song, W. Feng, and et al.,  _Synthesis of antisymmetric spin exchange interaction and chiral spin clusters in superconducting circuits,_Nat. Phys. 15, 382 (2019).
  * Li  _et al._ (2018b) X. Li, Y. Ma, J. Han, T. Chen, Y. Xu, W. Cai, H. Wang, Y. Song, Z.-Y. Xue, Z.-q. Yin, and L. Sun,  _Perfect quantum state transfer in a superconducting qubit chain with parametrically tunable couplings,_Phys. Rev. Appl. 10, 054009 (2018b).
  * Eckardt (2017) A. Eckardt,  _Colloquium: Atomic quantum gases in periodically driven optical lattices,_Rev. Mod. Phys. 89, 011004 (2017).
  * Weitenberg and Simonet (2021) C. Weitenberg and J. Simonet,  _Tailoring quantum gases by floquet engineering,_Nat. Phys. 17, 1342 (2021).
  * Garziano  _et al._ (2016) L. Garziano, V. Macrì, R. Stassi, O. Di Stefano, F. Nori, and S. Savasta,  _One photon can simultaneously excite two or more atoms,_Phys. Rev. Lett. 117, 043601 (2016).
  * Ma and Law (2015) K. K. W. Ma and C. K. Law,  _Three-photon resonance and adiabatic passage in the large-detuning rabi model,_Phys. Rev. A 92, 023842 (2015).
  * Kaufman  _et al._ (2020) B. Kaufman, T. Rozgonyi, P. Marquetand, and T. Weinacht,  _Adiabatic elimination in strong-field light-matter coupling,_Phys. Rev. A 102, 063117 (2020).
  * Combescot (2001) M. Combescot,  _On the generalized golden rule for transition probabilities,_J. Phys. A: Math. Gen. 34, 6087 (2001).
  * Koch  _et al._ (2010) J. Koch, A. A. Houck, K. L. Hur, and S. M. Girvin,  _Time-reversal-symmetry breaking in circuit-qed-based photon lattices,_Phys. Rev. A 82, 043811 (2010).
  * Xie  _et al._ (2020) J.-k. Xie, S.-l. Ma, and F.-l. Li,  _Quantum-interference-enhanced magnon blockade in an yttrium-iron-garnet sphere coupled to superconducting circuits,_Phys. Rev. A 101, 042331 (2020).
  * Liu  _et al._ (2019) Z.-X. Liu, H. Xiong, and Y. Wu,  _Magnon blockade in a hybrid ferromagnet-superconductor quantum system,_Phys. Rev. B 100, 134421 (2019).
  * Kounalakis  _et al._ (2022) M. Kounalakis, G. E. W. Bauer, and Y. M. Blanter,  _Analog quantum control of magnonic cat states on a chip by a superconducting qubit,_Phys. Rev. Lett. 129, 037205 (2022).
  * Wang  _et al._ (2016b) Y.-P. Wang, G.-Q. Zhang, D. Zhang, X.-Q. Luo, W. Xiong, S.-P. Wang, T.-F. Li, C.-M. Hu, and J. Q. You,  _Magnon kerr effect in a strongly coupled cavity-magnon system,_Phys. Rev. B 94, 224410 (2016b).
  * Qi and Jing (2022b) S.-f. Qi and J. Jing,  _Accelerated adiabatic passage in cavity magnomechanics,_Phys. Rev. A 105, 053710 (2022b).
  * Wootters (1998) W. K. Wootters,  _Entanglement of formation of an arbitrary state of two qubits,_Phys. Rev. Lett. 80, 2245 (1998).
  * Breuer and Petruccione (2022) H.-P. Breuer and F. Petruccione,  _The theory of open quantum systems,_ Oxford University Press, Oxford (2022).
  * Viola and Lloyd (1998) L. Viola and S. Lloyd,  _Dynamical suppression of decoherence in two-state quantum systems,_Phys. Rev. A 58, 2733 (1998).
  * Niemczyk  _et al._ (2010) T. Niemczyk, F. Deppe, H. Huebl, E. P. Menzel, F. Hocke, M. J. Schwarz, J. J. Garciaripoll, D. Zueco, T. Hömmer, and E. Solano,  _Circuit quantum electrodynamics in the ultrastrong-coupling regime,_Nat. Phys. 6, 772 (2010).
  * Rudner  _et al._ (2013) M. S. Rudner, N. H. Lindner, E. Berg, and M. Levin,  _Anomalous edge states and the bulk-edge correspondence for periodically driven two-dimensional systems,_Phys. Rev. X 3, 031005 (2013).


# Robust Optimal Control With Binary Adjustable Uncertainties

Yun Li1, Neil Yorke-Smith2 and Tamas Keviczky1 *The work was supported by the
Brains4Buildings project under the Dutch grant programme for Mission-Driven
Research, Development and Innovation (MOOI).1Yun Li and Tamas Keviczky are
with Delft Center for Systems and Control, Delft University of Technology,
Delft, the Netherlands<EMAIL_ADDRESS><EMAIL_ADDRESS>Yorke-
Smith is with STAR Lab, Delft University of Technology, the Netherlands.
<EMAIL_ADDRESS>
###### Abstract

Robust Optimal Control (ROC) with adjustable uncertainties has proven to be
effective in addressing critical challenges within modern energy networks,
especially the reserve and provision problem. However, prior research on ROC
with adjustable uncertainties has predominantly focused on the scenario of
uncertainties modeled as continuous variables. In this paper, we explore ROC
with binary adjustable uncertainties, where the uncertainties are modeled by
binary decision variables, marking the first investigation of its kind. To
tackle this new challenge, firstly we introduce a metric designed to
quantitatively measure the extent of binary adjustable uncertainties. Then, to
balance computational tractability and adaptability, we restrict control
policies to be affine functions with respect to uncertainties, and propose a
general design framework for ROC with binary adjustable uncertainties. To
address the inherent computational demands of the original ROC problem,
especially in large-scale applications, we employ strong duality (SD) and big-
M-based reformulations to create a scalable and computationally efficient
Mixed-Integer Linear Programming (MILP) formulation. Numerical simulations are
conducted to showcase the performance of our proposed approach, demonstrating
its applicability and effectiveness in handling binary adjustable
uncertainties within the context of modern energy networks.

##  I INTRODUCTION

In recent years, the so-called reserve and provision problem has gained
increasing attention owing to its ability to characterize the demand side
management issues in modern energy networks. As an efficient approach for
modeling the reserve and provision problem, a new optimal control framework
called robust optimal control (ROC) with adjustable uncertainties was proposed
[1, 2, 3, 4, 5, 6]. Unlike the conventional robust optimal control problem,
where the uncertainty sets are fixed and are determined by exogenous factors,
ROC with adjustable uncertainties considers the case where the scope of
uncertainties is adjustable and is determined by decision makers.

By allowing the decision maker to actively decide the admissible scope of
uncertainty, the corresponding robust optimal control framework can be more
versatile to include more practical cases. However, this versatility comes
with more computational challenges. It is worth pointing out that existing
research about ROC with adjustable uncertainties all focus on the case of
continuous adjustable uncertainties. In the continuous case, uncertainty sets
are typically represented as polyhedra or norm balls. As continuous variables
yield an infinite number of uncertain scenarios within a given set, this leads
to a semi-infinite optimization problem in the ROC formulation with adjustable
uncertainties. To create a computationally tractable reformulation, strong
duality theory is employed, reducing the need to satisfy universal constraints
for all uncertain scenarios and instead focusing only on the worst-case
scenario [1, 2, 3, 4, 5, 6].

Despite the success of ROC with continuous adjustable uncertainties in
supporting modern energy networks, it should be noted that it fails to deal
with the case where uncertainties can be binary. For example, in building
climate control, some heating, ventilation and air conditioning (HVAC) devices
only operate with two modes: on and off. If such devices are considered for
providing demand response services, the existing ROC with continuous
adjustable uncertainties fails to characterize this scenario. In addition, the
existing design framework for ROC with continuous adjustable uncertainties
cannot be directly generalized to deal with the binary case. On the one hand,
polyhedra or norm balls used in the continuous case to describe the scope of
uncertainties is not suitable for binary uncertainties. On the other hand,
duality-based reformulation, which is built upon the assumption of continuous
uncertainties, utilized in the continuous case for finding the worst
uncertainty scenario cannot be directly applied to deal with the binary case.
Motivated by this fact, in this paper, we investigate the problem of ROC with
binary adjustable uncertainties. Our main contributions can be summarized as
follows:

  * •

The problem of ROC with binary adjustable uncertainties is investigated for
the first time. By introducing a new metric to evaluate the scope of binary
uncertainties, we propose a general design formulation for ROC with binary
adjustable uncertainties.

  * •

Due to the high computational demand and poor scalability of the orignal
formulation of ROC with binary adjustable uncertainties, we provide a
computationally efficient and scalable alternative together with detailed
derivation procedures. In addition, the probabilistic robustness of the
optimal solution for our proposed approach is analyzed via Markov inequality.

  * •

We consider a practical case of reserve and provision problem: utilizing the
energy flexibility of buildings to provide demand response services and show
how our proposed approach can be utilized to explore the energy flexibility
potential. Simulation results are provided to demonstrate the effectiveness of
the proposed approach.

The remaining parts of this paper are organized as follows. Section II
introduces the control problem and proposes a novel formulation of the ROC
with binary adjustable uncertainties. For this original formulation, Section
III gives a computationally efficient alternative with detailed derivation
procedures and also analyzes the probabilistic robustness of the optimal
solution. Simulation results are presented in Section IV, and the conclusion
follows in Section V.

Notation: Bold-face letters are used to denote the stacked time sequences of
related signals; the operator $|\cdot|$ for a given set of indices denotes the
cardinality of the set; subscript $i$ of a given vector/matrix represents the
$i$-th element/row of the corresponding vector/matrix; $\mathbb{N}_{[0,a]}$
denotes the set of nonnegative integers less than equal to $a$; $\max$ or
$\min$ for a given vector implies row-wise maximization/minimization.
$\sum_{k}x_{k}$ denotes the summation of all terms $x_{k}$, $\mathbf{1}_{N}$
denotes the $N$-dimensional all 1 vector.

##  II Problem Formulation

In this section, we will introduce the problem of robust optimal control with
binary adjustable uncertainties. The considered system is

$x(t+1)=Ax(t)+Br(t)+Du(t)+Ev(t)$ (1)

where the $x(t)\in\mathbb{R}^{n}$ are system states, $r(t)\in\mathbb{B}^{m}$
are reference control inputs, $u(t)\in\mathbb{R}^{p}$ and
$v(t)\in\mathbb{B}^{q}$ are continuous and binary recourse control variables,
respectively. For system (1), we assume that the following linear state and
input constraints should be satisfied

$\displaystyle G_{x}x(t)\leq g_{x},$ (2a) $\displaystyle
G_{r}r(t)+G_{u}u(t)+G_{v}v(t)\leq g_{r}.$ (2b)

In our design framework, the reference control input is assumed to be
predetermined but suffers from adjustable uncertainties, which is equivalent
to the case that the system is subject to additive uncertainties.

Within the context of the reserve and provision problem, system (1) can model
the dynamics of several utilities in energy systems, such as building thermal
dynamics or energy storage devices. For example, for building climate control,
system (1) can model the thermal dynamics of buildings that draw power from
both main power grid $r(t)$, which is usually determined according to the day-
ahead electricity market, and local renewable energy sources (RES) ($u(t)$ and
$v(t)$), which can be flexibly scheduled in real time.

Based on (1), the predicted state evolution over $N$ time steps is

$\mathbf{x}=F_{x}x(0)+F_{r}\mathbf{r}+F_{u}\mathbf{u}+F_{v}\mathbf{v},$ (3)

where
$\mathbf{x}=[x(1)^{\mathrm{T}},\cdots,x(N)^{\mathrm{T}}]^{\mathrm{T}}\in\mathbb{R}^{nN}$,
$\mathbf{r}=[r(0)^{\mathrm{T}},\cdots,r(N-1)^{\mathrm{T}}]^{\mathrm{T}}\in\mathbb{B}^{mN}$,
$\mathbf{u}=[u(0)^{\mathrm{T}},\cdots,u(N-1)^{\mathrm{T}}]\in\mathbb{R}^{pN}$,
and
$\mathbf{v}=[v(0)^{\mathrm{T}},\cdots,v(N-1)^{\mathrm{T}}]^{\mathrm{T}}\in\mathbb{B}^{qN}$,
and $x^{0}$ is the initial state vector. The detailed format of the matrixes
$F_{x}$, $F_{r}$, $F_{u}$ and $F_{v}$ can be found in [7].

Correspondingly, constraints (2) within the prediction horizon can be
compactly denoted as

$\displaystyle\mathbb{X}^{N}:=\\{\mathbf{x}|G_{\mathbf{x}}\mathbf{x}\leq
g_{\mathbf{x}}\\},$ (4a)
$\displaystyle\mathbb{U}^{N}:=\\{(\mathbf{u},\mathbf{v})|G_{\bf
r}\mathbf{r}+G_{\mathbf{u}}\mathbf{u}+G_{\mathbf{v}}\mathbf{v}\leq g_{\bf
r}\\},$ (4b)

where $G_{\bf x}=\text{diag}(G_{x},\cdots,G_{x})$, $G_{\bf
r}=\text{diag}(G_{r},\cdots,G_{r})$, $G_{\bf
u}=\text{diag}(G_{u},\cdots,G_{u})$, $G_{\bf
v}=\text{diag}(G_{v},\cdots,G_{v})$, $g_{\bf x}=\mathbf{1}_{N}\otimes g_{x}$,
and $g_{\bf r}=\mathbf{1}_{N}\otimes g_{r}$.

In the reserve and provision framework, the reference control signal suffers
from unknown external manipulations to provide services for external entities,
and the design objective is to optimally quantify how much uncertainty (also
called flexibility) can be reserved while guaranteeing system constraints. For
the case of continuous adjustable uncertainties, the scope of flexibility is
denoted as adjustable polyhedrons or norm balls, and the volume of uncertainty
sets are adopted as metrics to measure the scope of uncertainties[2]. However,
in the case of binary uncertainty, this definition of uncertainty set is not
applicable because polyhedron/norm balls cannot accurately characterize
uncertain binary variables, and their volumes also cannot truly reflect the
scope of uncertainties.

In the following, we introduce a new metric to deal with binary adjustable
uncertainties. The indices $k\in\mathbb{N}_{[0,mN-1]}$ for $\mathbf{r}_{k}$
within the prediction horizon are partitioned into two disjunctive groups
$\mathcal{C}$ and $\mathcal{U}$. For all $k\in\mathcal{U}$, the corresponding
reference inputs $\mathbf{r}_{k}$ are allowed to be adjusted flexibly. For all
$k\in\mathcal{C}$, the reference inputs $\mathbf{r}_{k}$ are not influenced by
uncertainties and are fixed to their nominal values $\mathbf{\bar{r}}_{k}$.

However, due to system constraints, not all reference inputs
$\mathbf{r}_{k}(k\in\mathcal{U})$ can be flexibly adjusted without violating
(2). Here, we define $\mathcal{S}\subseteq\mathcal{U}$ as the index set in
which the corresponding $\mathbf{r}_{k}$ ($k\in\mathcal{S}$) is adjusted from
$\bar{\bf r}_{k}$ to a new value, denoted as $\tilde{\bf r}_{k}$ with
$\tilde{\bf r}_{k}\neq\bar{\bf r}_{k}$, without violating system constraints.
Based on the above definition, the values of $\mathbf{r}_{k}$ within the
prediction horizon can be represented as

$\mathbf{r}_{k}=\begin{cases}\tilde{\mathbf{r}}_{k},&k\in\mathcal{S}\\\
\bar{\mathbf{r}}_{k},&k\in\mathcal{C}\cup\mathcal{U}\setminus\mathcal{S}\end{cases}$
(5)

A schematic diagram of flexible reference input signals with $m=1$ is given in
Fig. 1. Similarly to the case of continuous adjustable uncertainty, $\bar{\bf
r}_{k}=0$ and $\tilde{\bf r}_{k}=1$ means up-ward flexibility, and $\bar{\bf
r}_{k}=1$ and $\tilde{\bf r}_{k}=0$ means down-ward flexibility.

Figure 1: Schemetic diagram of flexible reference input signal $\bf r$.

In contrast to the continuous case, where the volume of the uncertainty set is
adopted to measure the scope of reserved flexibility, according to the
definition (5), a new variable $\Gamma$ is defined as the largest cardinality
of $\mathcal{S}$ that can be freely selected among $\mathcal{U}$, to measure
the scope of the reserved uncertainty/flexibility. An admissible $\Gamma$
implies that there exist feasible recourse inputs $(\mathbf{u,v})$ to ensure
system constraint satisfaction for any possible reference input adjustment as
long as $|\mathcal{S}|\leq\Gamma$. Clearly, a larger value of $\Gamma$ implies
that the controlled system can reserve more flexibility for providing
services.

For the problem of ROC with adjustable uncertainties, a general control
objective is to balance the minimization of operational cost and the
maximization of the scope of reserved flexibility, while guaranteeing
constraint satisfaction under all possible scenarios of
uncertainties/flexibilities [1, 2]. Accordingly, the design objective for our
considered problem takes the form

$\displaystyle\min_{\mathbf{u,v},\Gamma}\ $
$\displaystyle\big{\\{}\max_{\mathcal{S}}\\{J(\mathbf{x},\mathbf{u},\mathbf{v})\\}-\lambda\Gamma\big{\\}}$
(6a) s.t.
$\displaystyle\mathbf{x}=F_{x}x^{0}+F_{r}\mathbf{r}+F_{u}\mathbf{u}+F_{v}\mathbf{v}$
(6b) $\displaystyle\mathbf{x}\in\mathbb{X}^{N},\
(\mathbf{u},\mathbf{v})\in\mathbb{U}^{N},$ (6c)
$\displaystyle\mathbf{r}_{k}=\bar{\bf r}_{k},\ \forall
k\in\mathcal{C}\cup\mathcal{U}\setminus\mathcal{S},$ (6d)
$\displaystyle\mathbf{r}_{k}=\tilde{\mathbf{r}}_{k},\ \forall
k\in\mathcal{S},$ (6e)
$\displaystyle\forall\mathcal{S}:|\mathcal{S}|\leq\Gamma,$ (6f)

where
$J(\mathbf{x},\mathbf{u},\mathbf{v}):=\sum_{k=0}^{N}l_{k}(\mathbf{x},\mathbf{u},\mathbf{v})$
is the summation of linear stage cost functions
$l_{k}(\mathbf{x},\mathbf{u},\mathbf{v})$, $\lambda>0$ is a user-selected
weighting parameter.

The objective function (6a) aims at minimizing the worst-case operational cost
given $\Gamma$, namely $\max_{\mathcal{S}}J(\mathbf{x,u,v})$, w.r.t. the input
uncertainties determined by $\mathcal{S}$ and in the meanwhile maximizing the
scope of flexibility $\Gamma$. If the worst-case operational cost function
$\max_{\mathcal{S}}J(\mathbf{x,u,v})$ is removed, optimization problem (6) is
to find out the largest feasible $\Gamma$. Namely, the largest cardinality of
$\mathcal{S}$ that can be freely selected without violating system
constraints. For example, in Fig. 1, four inputs among
$\\{\mathbf{r}_{k}:\forall k\in\mathcal{U}\\}$ are adjusted, if the optimal
$\Gamma^{*}$ computed in (6) satisfies $\Gamma^{*}\geq 4$, we are still
capable of ensuring system constraints for this scenario of flexible inputs.

To ensure the feasibility of (6), it is assumed that for the predefined
reference input $\bar{\mathbf{r}}$, there exists at least one input sequence
$(\mathbf{u,v})$ such that the input constraints (6c) are satisfied.

Remark 1: We refer to the optimal control problem (6), similarly to its
continuous counterparts: ROC with continuous adjustable uncertainties [1, 2],
as ROC with binary adjustable uncertainties. The main difference between ROC
with adjustable uncertainties and conventional ROC problem lies in the
uncertainty set. For conventional ROC problems, the uncertainty sets are
predetermined and fixed. But for the ROC with adjustable uncertainties, the
shape/size of uncertainty sets are not fixed but will be optimally determined
by decision variables, such as $\Gamma$ in (6).

Remark 2: We highlight that the problem formulation proposed in (6) can be
applied to several important issues in modern energy systems. For example,
similarly to [3, 5], formulation (6) can model the reserve and provision
problem of buildings with on-off types of HVAC devices for providing demand
response services. Besides, as shown in [8], it is also possible to extend the
formulation (6) to analyze and design resilient control strategies for power
grid networks to mitigate faults or cyberattacks.

##  III Robust Optimal Control Design and Robustness Analysis

This section is devoted to solving the problem of ROC with binary adjustable
uncertainties formulated in (6), and to analyzing the robustness of the
derived solution.

###  III-A Robust Optimal Control Design

Considering the system dynamics (1), state and input constraints (4), and the
linear property of the objective function, the optimization problem in (6) can
be rewritten in the following compacted form

$\displaystyle\min_{\theta,\Gamma,\mathbf{u,v}}\ $
$\displaystyle\theta-\lambda\Gamma$ (7a) s.t.
$\displaystyle\mathbf{O}\mathbf{r}+\mathbf{P}\mathbf{u}+\mathbf{Q}\mathbf{v}\leq\mathbf{h},$
(7b) $\displaystyle\mathbf{r}_{k}=\bar{\bf r}_{k},\ \forall
k\in\mathcal{C}\cup\mathcal{U}\setminus\mathcal{S},$ (7c)
$\displaystyle\mathbf{r}_{k}=\tilde{\bf r}_{k},\ \forall k\in\mathcal{S},$
(7d) $\displaystyle\forall\ \mathcal{S}:\ |\mathcal{S}|\leq\Gamma,$ (7e)

where matrixes $\mathbf{O}$, $\mathbf{P}$, $\mathbf{Q}$ and $\mathbf{h}$ are
constructed according to system dynamics, constraints and problem data,
$\theta$ is introduced for the epigraph reformulation of the worst-case
operational cost in (6) as robust constraints $J(\mathbf{x,u,v})\leq\theta$
$\forall\mathcal{S}:|\mathcal{S}|\leq\Gamma$, which are compactly included in
(7b). Similar reformulation can be found in [2].

Unlike the continuous case, in which the original formulation of ROC with
adjustable uncertainties is semi-infinite and is computationally intractable,
the optimization problem (7) can be solved directly without further
reformulation since it only entails a finite number of uncertainties
$\tilde{\bf r}_{k}$ and $\mathcal{S}$ given a fixed $\Gamma$. By trying out
all possible $\Gamma$ and considering all corresponding uncertain scenarios,
the optimal solution can be found by exhaustive search. However, it should be
noted that given a fixed $\Gamma$, the number of constraints that need to be
considered in (7) is proportional to
$\sum_{i=0}^{\Gamma}\binom{|\mathcal{U}|}{i}\cdot N$, which can result in a
large scale MILP that is computationally demanding. In the following, we will
propose a novel reformulation for (7) such that the resulting problem is
scalable and computationally efficient.

For the recourse control inputs $\bf u$ and $\bf v$, we will design closed-
loop control policies instead of open-loop control actions since it has been
demonstrated in many works that closed-loop control policies can tolerate more
uncertainties than open-loop control actions [9, 10, 11]. Because considering
arbitrary control policies can make the problem computationally intractable,
we consider the following affine control policies to balance robustness,
optimality, and computational effort for our proposed scheme

$\displaystyle\mathbf{u}_{k}=\sum_{i\in\mathcal{U}}M_{ki}\mathbf{r}_{i}+\eta_{k},$
(8a)
$\displaystyle\mathbf{v}_{k}=\sum_{i\in\mathcal{U}}L_{ki}\mathbf{r}_{i}+\epsilon_{k},$
(8b)

where $M_{ki}\in\mathbb{R}$, $\eta_{k}\in\mathbb{R}$, $L_{ki}\in\mathbb{Z}$
and $\epsilon_{k}\in\mathbb{Z}$ are control policy parameters to be optimized.
Since the exact value of the flexible inputs $\tilde{\bf
r}_{k}(k\in\mathcal{U})$ is not determined when solving (7), the control
policy (8) will adapt inputs $\mathbf{u}_{k}$ and $\mathbf{v}_{k}$ once
${\mathbf{r}}_{i}(i\in\mathcal{U})$ are revealed, which can increase the
robustness of the controlled system so that more flexibility can be reserved.
Note that binary restriction for $\mathbf{v}_{k}$ can be imposed by a linear
constraint $0\leq\mathbf{v}_{k}\leq 1$ due to the fact that the control policy
$\eqref{eq:ploicy_v}$ will only generate integer outputs.

Theorem 1: Considering the system dynamics (1), the binary adjustable input
uncertainties (5) and the control policies (8), the robust optimal control
problem in (7) can be reformulated as the following MILP problem

$\displaystyle\min_{\begin{subarray}{c}\theta,\delta_{k},\mu_{ij},\phi_{ij},y_{ij}\beta_{ij},\\\
M_{kj},L_{kj},\eta_{k},\epsilon_{k}\end{subarray}}\ $
$\displaystyle\theta-\lambda\sum_{k\in\mathcal{U}}\delta_{k}$ (9a) s.t.
$\displaystyle\sum_{j\in\mathcal{U}}\mu_{ij}+\sum_{j\in\mathcal{U}}\beta_{ij}+\phi_{i}+\sum_{k}\mathbf{P}_{ik}\eta_{k}+$
$\displaystyle\qquad\sum_{k}\mathbf{Q}_{ik}\epsilon_{k}+\sum_{j\in\mathcal{C}}\mathbf{O}_{ij}\bar{\mathbf{r}}_{j}\leq\mathbf{h}_{i},$
(9b)
$\displaystyle\mu_{ij}+\pi_{i}\geq\frac{1}{2}\bigg{(}\mathbf{O}_{ij}+\sum_{k}\mathbf{P}_{ik}M_{kj}+$
$\displaystyle\qquad\sum_{k}\mathbf{Q}_{ik}L_{kj}+y_{ij}\bigg{)}+\bigg{(}\mathbf{O}_{ij}+$
$\displaystyle\qquad\sum_{k}\mathbf{P}_{ik}M_{kj}+\sum_{k}\mathbf{Q}_{ik}L_{kj}\bigg{)}\bar{\bf
r}_{j},$ (9c) $\displaystyle-y_{ij}\leq\mathbf{O}_{ij}+\sum_{k}{\bf
P}_{ik}M_{kj}+$ $\displaystyle\qquad\sum_{k}\mathbf{Q}_{ik}L_{kj}\leq y_{ij},$
(9d)
$\displaystyle\phi_{i}\geq\sum_{j\in\mathcal{U}}\bigg{(}\mathbf{O}_{ij}+\sum_{k}\mathbf{P}_{ik}M_{kj}+$
$\displaystyle\qquad\sum_{k}\mathbf{Q}_{ik}L_{kj}\bigg{)}\cdot\mathbf{\bar{r}}_{j},$
(9e) $\displaystyle 0\leq\beta_{ij}\leq M\delta_{j},$ (9f) $\displaystyle
0\leq\pi_{i}-\beta_{ij}\leq M(1-\delta_{j}),$ (9g)
$\displaystyle\delta_{j}\in\mathbb{B},$ (9h) $\displaystyle\mu_{ij}\geq 0,\
\pi_{i}\geq 0,\ y_{ij}\geq 0,$ (9i) $\displaystyle i\in\mathcal{I},\
j\in\mathcal{U},$ (9j)

where $\mathcal{I}$ is the set of row indices of constraints (7b), and the
optimal solution $\delta_{k}^{*}$ satisfies
$\sum_{j\in\mathcal{U}}\delta_{j}^{*}=\Gamma^{*}$.

Proof: Given a fixed $\Gamma$, system constraint (7b) has to be satisfied for
all possible $\tilde{\mathbf{r}}_{k}$ and $\mathcal{S}$ such that
$k\in\mathcal{S}$ and $|\mathcal{S}|\leq\Gamma$. To ensure constraint
satisfaction for all possible uncertainties, which entails the order of
$\sum_{i=0}^{\Gamma}\binom{|\mathcal{U}|}{i}$ uncertainty scenarios, we
alternatively design a control scheme to guarantee constraint satisfaction
under the worst case scenario of uncertainty, which can be represented as

$\max_{\mathcal{S}}\big{\\{}\mathbf{Or+Pu+Qv}:|\mathcal{S}|\leq\Gamma,\mathbf{r}_{j}\text{
satisfies \eqref{eq:uncertainty}}\big{\\}}\leq\mathbf{h},$ (10)

where the maximization and inequality in (10) is computed and hold row-wise,
respectively. The $i$-th row of constraints (10), combining the control
policies (8) and the reference input pattern (5) can be expressed as

$\begin{split}\max_{\mathcal{S}:|\mathcal{S}|\leq\Gamma}\
&\bigg{\\{}\mathbf{O}_{i}^{\mathrm{T}}\mathbf{r}+\sum_{k}\mathbf{P}_{ik}\bigg{(}\sum_{j\in\mathcal{U}}M_{kj}\mathbf{r}_{j}+\eta_{k}\bigg{)}+\\\
&\sum_{k}\mathbf{Q}_{ik}\bigg{(}\sum_{j\in\mathcal{U}}L_{kj}\mathbf{r}_{j}+\epsilon_{k}\bigg{)}\bigg{\\}}\leq\mathbf{h}_{i},\end{split}$

which can be further reformulated as

$\displaystyle\max_{\mathcal{S}:|\mathcal{S}|\leq\Gamma}\bigg{\\{}\sum_{j\in\mathcal{S}}\bigg{(}\mathbf{O}_{ij}+\sum_{k}\mathbf{P}_{ik}M_{kj}+\sum_{k}\mathbf{Q}_{ik}L_{kj}\bigg{)}\tilde{\bf
r}_{j}$
$\displaystyle\qquad+\sum_{j\in\mathcal{U}\setminus\mathcal{S}}\bigg{(}\mathbf{O}_{ij}+\sum_{k}\mathbf{P}_{ik}M_{kj}+\sum_{k}\mathbf{Q}_{ik}L_{kj}\bigg{)}\bar{\bf
r}_{j}\bigg{\\}}$
$\displaystyle\qquad+\sum_{k}\mathbf{P}_{ik}\eta_{k}+\sum_{k}\mathbf{Q}_{ik}\epsilon_{k}+\sum_{j\in\mathcal{C}}\mathbf{O}_{ij}\bar{\mathbf{r}}_{j}\leq\mathbf{h}_{i}.$
(11)

As explained in [12], it can be verified that

$\begin{split}&\bigg{(}\mathbf{O}_{ij}+\sum_{k}\mathbf{P}_{ik}M_{kj}+\sum_{k}\mathbf{Q}_{ik}L_{kj}\bigg{)}\tilde{\bf
r}_{j}\leq\\\
&\qquad\frac{1}{2}\bigg{(}\mathbf{O}_{ij}+\sum_{k}\mathbf{P}_{ik}M_{kj}+\sum_{k}\mathbf{Q}_{ik}L_{kj}+\big{|}\mathbf{O}_{ij}+\\\
&\qquad\sum_{k}\mathbf{P}_{ik}M_{kj}+\sum_{k}\mathbf{Q}_{ik}L_{kj}\big{|}\bigg{)}.\end{split}$
(12)

Accordingly, by following a similar line as in [13] and introducing additional
decision variables $z_{j}\in\mathbb{R}\ (j\in\mathcal{U})$, the maximization
problem in (III-A) is relaxed as

$\displaystyle\max_{z_{j}}\ $
$\displaystyle\sum_{j\in\mathcal{U}}\frac{z_{j}}{2}\cdot\bigg{(}\mathbf{O}_{ij}+\sum_{k}\mathbf{P}_{ik}M_{kj}+\sum_{k}\mathbf{Q}_{ik}L_{kj}+$
$\displaystyle\quad|\mathbf{O}_{ij}+\sum_{k}\mathbf{P}_{ik}M_{kj}+\sum_{k}\mathbf{Q}_{ik}L_{kj}|\bigg{)}+$
$\displaystyle\quad(1-z_{j})\cdot\bigg{(}\mathbf{O}_{ij}+\sum_{k}\mathbf{P}_{ik}M_{kj}+\sum_{k}\mathbf{Q}_{ik}L_{kj}\bigg{)}\bar{\bf
r}_{j}$ (13a) s.t. $\displaystyle 0\leq z_{j}\leq
1,\quad\sum_{j\in\mathcal{U}}z_{j}\leq\Gamma.$ (13b)

Since (13) is a linear program (LP) w.r.t. $z_{j}$, according to strong
duality theory of LP, see [14], the objective value of (13) coincides with
that of the following dual problem

$\displaystyle\min_{\mu_{ij},\pi_{i},\phi_{i}}\ $
$\displaystyle\sum_{j\in\mathcal{U}}\mu_{ij}+\Gamma\pi_{i}+\phi_{i}$ (14a)
s.t.
$\displaystyle\mu_{ij}+\pi_{i}\geq\frac{1}{2}\bigg{(}\mathbf{O}_{ij}+\sum_{k}\mathbf{P}_{ik}M_{kj}+\sum_{k}\mathbf{Q}_{ik}L_{kj}$
$\displaystyle\quad+|\mathbf{O}_{ij}+\sum_{k}\mathbf{P}_{ik}M_{kj}+\sum_{k}\mathbf{Q}_{ik}L_{kj}|\bigg{)}-$
$\displaystyle\quad\bigg{(}\mathbf{O}_{ij}+\sum_{k}\mathbf{P}_{ik}M_{kj}+\sum_{k}\mathbf{Q}_{ik}L_{kj}\bigg{)}\bar{\bf
r}_{j},$ (14b)
$\displaystyle\phi_{i}\geq\sum_{j\in\mathcal{U}}\bigg{(}\mathbf{O}_{ij}+\sum_{k}\mathbf{P}_{ik}M_{kj}+\sum_{k}\mathbf{Q}_{ik}L_{kj}\bigg{)}\mathbf{\bar{r}}_{j},$
(14c) $\displaystyle\mu_{ij}\geq 0,\ \pi_{i}\geq 0,\ \forall j\in\mathcal{U},$
(14d)

where $\mu_{ij}$ and $\pi_{i}$ are Lagrange multipliers.

Optimization problem (14) is nonlinear because it contains the absolute
function, which will incur increased computational effort. By introducing
auxiliary decision variables $y_{ij}\geq 0$, optimization problem (14) can be
further relaxed as the following linear optimization problem

$\displaystyle\min_{\mu_{ij},\pi_{i},\phi_{i}}\ $
$\displaystyle\sum_{j\in\mathcal{U}}\mu_{ij}+\Gamma\pi_{i}+\phi_{i}$ (15a)
s.t.
$\displaystyle\mu_{ij}+\pi_{i}\geq\frac{\mathbf{O}_{ij}+\sum_{k}\mathbf{P}_{ik}M_{kj}+\sum_{k}\mathbf{Q}_{ik}L_{kj}+y_{ij}}{2}$
$\displaystyle\quad-\bigg{(}\mathbf{O}_{ij}+\sum_{k}\mathbf{P}_{ik}M_{kj}+\sum_{k}\mathbf{Q}_{ik}L_{kj}\bigg{)}\cdot\bar{\bf
r}_{j},$ (15b) $\displaystyle-
y_{ij}\leq\mathbf{O}_{ij}+\sum_{k}\mathbf{P}_{ik}M_{kj}+\sum_{k}\mathbf{Q}_{ik}L_{kj}\leq
y_{ij},$ (15c)
$\displaystyle\phi_{i}\geq\sum_{j\in\mathcal{U}}\bigg{(}\mathbf{O}_{ij}+\sum_{k}\mathbf{P}_{ik}M_{kj}+\sum_{k}\mathbf{Q}_{ik}L_{kj}\bigg{)}\cdot\mathbf{\bar{r}}_{j},$
(15d) $\displaystyle\mu_{ij}\geq 0,\ \pi_{i}\geq 0,\ y_{ij}\geq 0,\ \forall
j\in\mathcal{U}.$ (15e)

For the inequality (III-A), the maximization problem can be relaxed by a
feasible solution of the above LP problem to yield the following alternative
constraints

$\displaystyle\sum_{j\in\mathcal{U}}\mu_{ij}+\Gamma\pi_{i}+\phi_{i}+\sum_{k}\mathbf{P}_{ik}\eta_{k}+\sum_{k}\mathbf{Q}_{ik}\epsilon_{k}+$
$\displaystyle\qquad\sum_{j\in\mathcal{C}}\mathbf{O}_{ij}\bar{\bf
r}_{j}\leq\mathbf{h}_{i},$ (16a)
$\displaystyle\mu_{ij}+\pi_{i}\geq\frac{1}{2}\bigg{(}\mathbf{O}_{ij}+\sum_{k}\mathbf{P}_{ik}M_{kj}+\sum_{k}\mathbf{Q}_{ik}L_{kj}+y_{ij}\bigg{)}$
$\displaystyle\quad-\bigg{(}\mathbf{O}_{ij}+\sum_{k}\mathbf{P}_{ik}M_{kj}+\sum_{k}\mathbf{Q}_{ik}L_{kj}\bigg{)}\cdot\bar{\bf
r}_{j},$ (16b) $\displaystyle-
y_{ij}\leq\mathbf{O}_{ij}+\sum_{k}\mathbf{P}_{ik}M_{kj}+\sum_{k}\mathbf{Q}_{ik}L_{kj}\leq
y_{ij},$ (16c)
$\displaystyle\phi_{i}\geq\sum_{j\in\mathcal{U}}\bigg{(}\mathbf{O}_{ij}+\sum_{k}\mathbf{P}_{ik}M_{kj}+\sum_{k}\mathbf{Q}_{ik}L_{kj}\bigg{)}\cdot\mathbf{\bar{r}}_{j},$
(16d) $\displaystyle\mu_{ij}\geq 0,\ \pi_{i}\geq 0,\ y_{ij}\geq 0,\ \forall
j\in\mathcal{U}.$ (16e)

Furthermore, introducing $\delta_{j}\in\mathbb{B}\ (j\in\mathcal{U})$ and
$\beta_{ij}\in\mathbb{R}\ (i\in\mathcal{I},j\in\mathcal{U})$ and defining
$\Gamma:=\sum_{j\in\mathcal{U}}\delta_{j}$, the term $\Gamma\pi_{i}$ in (16a)
can be further relaxed via big-M formulation as

$\displaystyle\Gamma\pi_{i}=\sum_{j\in\mathcal{U}}\beta_{ij},$ (17a)
$\displaystyle 0\leq\beta_{ij}\leq M\delta_{j},$ (17b) $\displaystyle
0\leq\pi_{i}-\beta_{ij}\leq M(1-\delta_{j}),$ (17c)
$\displaystyle\delta_{j}\in\mathbb{B},$ (17d)

where $M>0$ is a sufficiently large constant. Finally, constraints (7b) – (7e)
can be replaced by (16) and (17), which result in the optimization problem
(9). This completes the proof. $\square$

Remark 3: For the original formulation in (7), given a fixed $\Gamma$, in
order to guarantee constraint satisfaction robustly, the number of constraints
considered is proportional to
$\sum_{i=0}^{\Gamma}\binom{|\mathcal{U}|}{i}\cdot N$. In contrast, the number
of constraints in the reformulated optimization problem in (9) is only
proportional to $|\mathcal{U}|\cdot N$, which is more favourable when
$\sum_{i=0}^{\Gamma}\binom{|\mathcal{U}|}{i}\gg|\mathcal{U}|$. In addition,
the constraint (16a) has nonconvex term $\Gamma\pi_{i}$ since
$\Gamma\in\mathbb{N}_{[0,|\mathcal{U}|]}$ is a decision variable. While this
nonconvexity in (16a) can be dealt with by several solvers, e.g. Gurobi, with
the big-M based reformualation in (17) the final optimization problem (9),
which is a MILP, can be solved more efficiently by more off-the-shelf solvers,
such as Gurobi, CPLEX, and GLPK.

###  III-B Robustness of the Optimal Solution

The obtained optimal control strategy, namely
$(M^{*}_{kj},\eta^{*}_{k},L_{kj}^{*},\epsilon_{k}^{*})$, is able to guarantee
constraint satisfaction if no more than $\Gamma$ of the flexible reference
inputs ${\mathbf{r}}_{j}\ (j\in\mathcal{U})$ are changed. In practice,
however, it is possible that in some cases more than $\Gamma$ number of
reference inputs $\mathbf{r}_{k}$ are subjected to uncertainties, e.g., in
case of device faults. As a result, it is of interest and importance to
analyze the robustness of the derived solution for such cases. In the
following, we will investigate the possibility of constraint violation with
the derived optimal solution when more than $\Gamma$ number of
${\mathbf{r}}_{j}\ (j\in\mathcal{U})$ might be changed.

In our analysis, $\mathbf{r}_{j}\ (j\in\mathcal{U})$ are regarded as random
variables (r.v.). Consequently, the left-hand side of constraint (7b) are also
r.v. For the $i$-th row of the constraint (7b), its probability of violation
is

$\mathbb{P}_{\text{vio}}=\text{Pr}\big{(}\mathbf{O}_{i}^{\mathrm{T}}\mathbf{r}+\sum_{j\in\mathcal{U}}\mathbf{r}_{j}\sum_{k}\mathbf{P}_{ik}M^{*}_{kj}+\sum_{j\in\mathcal{U}}\mathbf{r}_{j}\sum_{k}\mathbf{Q}_{ik}L^{*}_{kj}+\\\
\sum_{k}\mathbf{P}_{ik}\eta^{*}_{k}+\sum_{k}\mathbf{Q}_{ik}\epsilon^{*}_{k}\geq\mathbf{h}_{i}\big{)}.$
(18)

For notational brevity, we define

$\displaystyle
a_{ij}:=\mathbf{O}_{ij}+\sum_{k}\mathbf{P}_{ik}M_{kj}^{*}+\sum_{k}\mathbf{Q}_{ik}L^{*}_{kj},$
(19a) $\displaystyle
b_{i}:=\mathbf{h}_{i}-\sum_{j\in\mathcal{C}}\mathbf{O}_{ij}{\bf
r}_{j}-\sum_{k}\mathbf{P}_{ik}\eta_{k}^{*}-\sum_{k}\mathbf{Q}_{ik}\epsilon_{k}^{*}.$
(19b)

Proposition 1: Assuming that $\mathbf{r}_{j}\in\mathbb{B}\ (j\in\mathcal{U})$
are independent random variables, then the probability of violation for the
$i$-th constraint (18) satisfies

$\mathbb{P}_{\text{vio}}\leq\frac{\prod_{j\in\mathcal{U}}\mathbb{E}\big{(}\exp(a_{ij}{\bf
r}_{j})\big{)}}{\exp(b_{i})}.$ (20)

Proof: Based on the definitions of $a_{ij}$ and $b_{j}$ in (19), the
probability of violation (18) can be rewritten as

$\displaystyle\mathbb{P}_{\text{vio}}$
$\displaystyle=\text{Pr}\bigg{[}\sum_{j\in\mathcal{U}}a_{ij}\mathbf{r}_{j}\geq
b_{i}\bigg{]}\ $ (21a)
$\displaystyle=\text{Pr}\left\\{\exp\big{(}\sum_{j\in\mathcal{U}}a_{ij}\mathbf{r}_{j}\big{)}\geq\exp(b_{i})\right\\}.$
(21b)

Then, it follows from the Markov inequality [15] and the independence property
of ${\bf r}_{j}$ that

$\mathbb{P}_{\text{vio}}\leq\frac{\prod_{j\in\mathcal{U}}\mathbb{E}\big{(}\exp(a_{ij}{\bf
r}_{j})\big{)}}{\exp(b_{i})}.$ (22)

This completes the proof. $\square$

It should be noted that the bound of the probability of violation in (20)
relies on system information in $a_{ij}$ and $b_{i}$. In the following, we
give a probability bound that does not rely on system information and is easy
to compute.

Proposition 2: Assuming that ${\bf r}_{j}\ (j\in\mathcal{U})$ are independent
random variables with $\text{Pr}({\bf r}_{j}\neq\bar{\bf
r}_{j})=\epsilon_{j}$, Then, we have
$\mathbb{P}_{\text{vio}}\leq\frac{\sum_{j\in\mathcal{U}}\epsilon_{j}}{\Gamma}$.

Proof: Define random variables $z_{j}$ as the indicator functions of the
events ${\bf r}_{j}\neq\bar{\bf r}_{j}\ (j\in\mathcal{U})$. Then, it is
readily concluded that
$\mathbb{P}_{\text{vio}}\leq\text{Pr}(\sum_{j\in\mathcal{U}}z_{j}\geq\Gamma)$.
By applying Markov inequality and the independence of ${\bf r}_{j}$, it yields

$\displaystyle\mathbb{P}_{\text{vio}}\ $
$\displaystyle\leq\frac{\mathbb{E}(\sum_{j\in\mathcal{U}}z_{j})}{\Gamma}$
(23a)
$\displaystyle=\frac{\sum_{j\in\mathcal{U}}\mathbb{E}(z_{j})}{\Gamma}=\frac{\sum_{j\in\mathcal{U}}\epsilon_{j}}{\Gamma}.$
(23b)

This completes the proof. $\square$

While the bound of probability violation in Proposition 2 is independent of
system information and is easy to compute, it is very conservative when
$\Gamma$ is small or $\sum_{j\in\mathcal{U}}\epsilon_{j}$ is large.

##  IV Simulation Results

Figure 2: Diagram of energy consumption options.

The effectiveness of our proposed approach is illustrated via a numerical case
study about the reserve and provision problem for a smart building system. A
schematic diagram of the considered building system is depicted in Fig. 2. In
our simulation, it is assumed that the thermal energy of the building is
provided by three heat pumps (HP) that draw electricity from two sectors:
power grid and local renewable energy sources (a solar panel). The HP powered
by the grid has only two operation modes: on and off, which is denoted as $\bf
r$, and the operation status is predetermined according to day-ahead
electricity price. For the remaining two HPs drawing power from local RES,
their control signals ($\bf u,v$) are available for adjusting in real-time.
One HP has only on-off operational modes, and the other can operate
continuously, e.g., PWM based operation. A solar panel is equipped to provide
RES.

As ubiquitous storage devices of energy flexibility, buildings are utilized to
provide demand response (DR) services to the power grid. During demand
response period, the building management system is asked to adjust its
electricity consumption, namely the operational pattern of $\bf r$, according
to the requests of the grid operator.

In our simulation, the sampling period is selected as 30 min and the
prediction horizon is set as $N=48$, i.e., one day. The starting time instant
of flexibility assessment period $\mathcal{U}$ is selected as $22$ (11am). A
solar panel with a maximal 2000W electricity output is included to provide
RES. The initial indoor temperature is selected as $x(0)=21^{\circ}$C, and the
indoor temperature bound is $[20^{\circ}\text{C},24^{\circ}\text{C}]$.

The objective for the optimization problem (7) is selected to maximize
$\Gamma$. Namely, we want to assess the largest flexibility potential of the
considered building system within the period of $\mathcal{U}$. Three schemes
are adopted in our simulation:

  * •

Scheme 1: our proposed formulation (9).

  * •

Scheme 2: original formulation (7) solved via exhaustive search.

  * •

Scheme 3: our proposed formulation (9) but only considering open-loop control
actions ($M_{ki}=0$, $L_{ki}=0$ in (8)).

All simulations are carried out using Gurobi 9.5.1 [16] on an Intel Xeon
W-2223 CPU at 3.60GHz with 16G RAM. Optimization problems are modeled via the
Python package Pyomo [17].

Figure 3: Computation time and optimal $\Gamma^{*}$ for all schemes. Figure 4:
Indoor temperature envelopes for all schemes.

Fig. 3 gives the computation time and the optimal $\Gamma^{*}$, i.e., the
largest number of flexible switches of $\mathbf{r}_{j}$ $(j\in\mathcal{U})$,
for all schemes with an increasing size of $|\mathcal{U}|$. Clearly, it can be
observed from Fig. 3a that the computation time of Scheme 2 (exhaustive
search) increases exponentially with $|\mathcal{U}|$. On the contrary, our
proposed approach is much more computationally efficient, and its computation
time increases linearly with $|\mathcal{U}|$. Compared with Scheme 1, Scheme
3, which only considers open-loop control actions instead of closed-loop
control policies as in Scheme 1, is computationally less demanding due to the
lower number of decision variables. Obviously, the downside of Scheme 3 is its
suboptimality, as shown in Fig. 3b. In addition, we can also conclude from
Fig. 3b that, for this case study, the affine control policy (8) restricted in
our proposed approach (Scheme 1) does not compromise the optimality of the
solution and can find the optimal $\Gamma^{*}$ as with Scheme 2.

Fig. 4 depicts the indoor temperature envelopes with all 3 schemes under all
possible uncertain scenarios of $\bf r$ given $|\mathcal{U}|=8$ and
$\Gamma=5$. Clearly, Scheme 3 is unable to robustly guarantee indoor comfort
constraints, which validates the conclusion that open-loop control actions are
less robust than closed-loop control policies.

##  V CONCLUSIONS

This paper investigated a class of ROC problems with binary adjustable
uncertainties, that offers a framework for addressing the reserve and
provision problem for energy systems with on-off type devices. In contrast to
the existing methods for ROC with continuous adjustable uncertainties, a novel
metric is introduced to measure the extent of binary uncertainties, and a
general design framework for ROC with binary adjustable uncertainties is
formulated. Subsequently, we reformulated and relaxed the ROC formulation by
introducing auxiliary variables and applying strong duality, transforming it
into a MILP problem amenable to efficient numerical solvers. This innovative
approach empowers us to quantitively evaluate the flexibility potential
inherent in modern energy systems with on-off type equipment. The efficacy of
our proposed methodology is demonstrated through numerical experiments.

Future extensions include considering nonlinear constraints and nonaffine
control policies, as well as hybrid types of adjustable uncertainties.

## References

  * [1] A. Bitlislioğlu, T. T. Gorecki, and C. N. Jones, “Robust tracking commitment,” _IEEE Transactions on Automatic Control_ , vol. 62, no. 9, pp. 4451–4466, 2017.
  * [2] X. Zhang, M. Kamgarpour, A. Georghiou, P. Goulart, and J. Lygeros, “Robust optimal control with adjustable uncertainty sets,” _Automatica_ , vol. 75, pp. 249–259, 2017.
  * [3] E. Vrettos, F. Oldewurtel, and G. Andersson, “Robust energy-constrained frequency reserves from aggregations of commercial buildings,” _IEEE Transactions on Power Systems_ , vol. 31, no. 6, pp. 4272–4285, 2016.
  * [4] O. Nohadani and K. Sharma, “Optimization under decision-dependent uncertainty,” _SIAM Journal on Optimization_ , vol. 28, no. 2, pp. 1773–1795, 2018.
  * [5] F. Bünning, J. Warrington, P. Heer, R. S. Smith, and J. Lygeros, “Robust MPC with data-driven demand forecasting for frequency regulation with heat pumps,” _Control Engineering Practice_ , vol. 122, p. 105101, 2022.
  * [6] N. H. Lappas and C. E. Gounaris, “Robust optimization for decision-making under endogenous uncertainty,” _Computers & Chemical Engineering_, vol. 111, pp. 252–266, 2018.
  * [7] Y. Li, N. Yorke-Smith, and T. Keviczky, “Robust optimal control with inexact state measurements and adjustable uncertainty sets,” _IFAC-PapersOnLine_ , 2023.
  * [8] L. Zhao and B. Zeng, “Vulnerability analysis of power grids with line switching,” _IEEE Transactions on Power Systems_ , vol. 28, no. 3, pp. 2727–2736, 2013.
  * [9] A. Ben-Tal, A. Goryashko, E. Guslitzer, and A. Nemirovski, “Adjustable robust solutions of uncertain linear programs,” _Mathematical programming_ , vol. 99, no. 2, pp. 351–376, 2004.
  * [10] D. Bertsimas and D. D. Hertog, _Robust and adaptive optimization_. Dynamic Ideas, 2022.
  * [11] D. Bertsimas and A. Georghiou, “Binary decision rules for multistage adaptive mixed-integer optimization,” _Mathematical Programming_ , vol. 167, pp. 395–433, 2018.
  * [12] J. Ramirez-Calderon, V. Jorge Leon, and B. Lawrence, “Robust binary linear programming under implementation uncertainty,” _Engineering Optimization_ , pp. 1–21, 2022.
  * [13] D. Bertsimas and M. Sim, “The price of robustness,” _Operations Research_ , vol. 52, no. 1, pp. 35–53, 2004.
  * [14] S. P. Boyd and L. Vandenberghe, _Convex optimization_. Cambridge University Press, 2004.
  * [15] J. Jacod and P. Protter, _Probability essentials_. Springer Science & Business Media, 2004.
  * [16] Gurobi Optimization, LLC, “Gurobi Optimizer Reference Manual,” 2022. [Online]. Available: https://www.gurobi.com
  * [17] M. L. Bynum, G. A. Hackebeil, W. E. Hart, C. D. Laird, B. L. Nicholson, J. D. Siirola, J.-P. Watson, and D. L. Woodruff, _Pyomo–optimization modeling in python_ , 3rd ed. Springer Science & Business Media, 2021, vol. 67.