# Nonreciprocal heat flux via synthetic fields in linear quantum systems

S.-A. Biehs<EMAIL_ADDRESS>Institut für Physik, Carl von Ossietzky
Universität, D-26111 Oldenburg, Germany  Center for Nanoscale Dynamics
(CeNaD), Carl von Ossietzky Universität, D-26129 Oldenburg, Germany
Laboratoire Charles Coulomb (L2C), UMR 5221 CNRS-University of Montpellier,
F-34095 Montpellier, France  P. Rodriguez-Lopez<EMAIL_ADDRESS>Área de
Electromagnetismo and Grupo Interdisciplinar de Sistemas Complejos (GISC),
Universidad Rey Juan Carlos, 28933 Móstoles, Madrid, Spain  Laboratoire
Charles Coulomb (L2C), UMR 5221 CNRS-University of Montpellier, F-34095
Montpellier, France  M. Antezza<EMAIL_ADDRESS>Laboratoire
Charles Coulomb (L2C), UMR 5221 CNRS-University of Montpellier, F-34095
Montpellier, France  Institut Universitaire de France, 1 rue Descartes, Paris
Cedex 05, F-75231, France  G. S. Agarwal<EMAIL_ADDRESS>Institute
for Quantum Science and Engineering and Department of Biological and
Agricultural Engineering Department of Physics and Astronomy, Texas A & M
University, College Station, Texas 77845, USA

###### Abstract

We study the heat transfer between N coupled quantum resonators with applied
synthetic electric and magnetic fields realized by changing the resonators
parameters by external drivings. To this end we develop two general methods,
based on the quantum optical master equation and on the Langevin equation for
$N$ coupled oscillators where all quantum oscillators can have their own heat
baths. The synthetic electric and magnetic fields are generated by a dynamical
modulation of the oscillator resonance with a given phase. Using Floquet
theory we solve the dynamical equations with both methods which allow us to
determine the heat flux spectra and the transferred power. With apply these
methods to study the specific case of a linear tight-binding chain of four
quantum coupled resonators. We find that in that case, in addition to a non-
reciprocal heat flux spectrum already predicted in previous investigations,
the synthetic fields induce here non-reciprocity in the total heat flux hence
realizing a net heat flux rectification.

##  I Introduction

In the last decade a great number of experiments have verified the near-field
enhancement of thermal radiation between two macroscopic objects down to
distances of a few nanometer HuEtAl2008 ; OttensEtAl2011 ; Kralik2012 ;
LimEtAl2015 ; WatjenEtAl2016 ; BernadiEtAl2016 ; SongEtAl2016 ;
FiorinoEtAl2018b ; Mittapally2023 . In particular, the theoretically proposed
effects of thermal rectification with a phase-change diode Ito ;
FiorinoEtAl2018a , a phase-change material based memory ItoMem  and active
heat flux switching or modulations Minnich ; Minnich2 ; ShiEtAl2022  have been
realized, experimentally. Also several proposals for non-reciprocal heat
transport have been made, but these effects have not been demonstrated
experimentally, yet. Typically, these proposals rely on the application of
magnetic fields to nanoscale setups involving magneto-optical materials or by
using Weyl semi-metals with intrinsic nonreciprocal optical properties. It can
be shown theoretically that by means of magnetic fields the magnitude of the
heat flux and its direction can be manipulated Latella2017 ; He2020 ; Cuevas ;
reviewMagneto ; MoncadaVilla1 ; MoncadaVilla2 ; GeEtAl2019 ; WuEtAl2021 . Due
to the broken time-reversal symmetry also non-reciprocal heat fluxes can exist
in such cases leading to persistent heat currents and fluxes zhufan ;
CircularHeatFlux , persistent angular momenta and spins Silveirinha ;
CircularHeatFlux ; Zubin2019 , normal and anomalous Hall effect for thermal
radiation hall ; ahall , and diode effects by coupling to non-reciprocal
surface modes Doyeux2017 ; NRdiode ; HuEtAl2023 , and spin-directional near-
and far-field thermal emission FanEtAl2020 ; DongEtAl2021 . A tradeoff of
using magneto-optical materials is that for having observable non-reciprocal
heat fluxes, experiments with large magnetic fields in a nanoscale setup are
necessary. On the other hand, using Weyl semi-metals with intrinsic
nonreciprocity does not allow for a dynamic tuning.

Figure 1:  Sketch of $N$ coupled quantum resonators each coupled to its own
heat bath.

Recently, the modulation of resonance frequencies of a system of resonators
with a single modulation frequency but different phases has been interpreted
as a mean to create synthetic electric and magnetic fields Lustig . For the
energy transmission in a setup of two resonators with applied synthetic
electric and magnetic fields, i.e. with a modulation of the resonance
frequencies and a phase shift, it could be shown experimentally and
theoretically that monochromatic waves are transmitted in a nonreciprocal
manner PetersonEtAl2019  if there is a non-zero phase shift, i.e. a synthetic
magnetic field. Now, if the two resonators with applied synthetic electric and
magnetic fields are coupled to two thermal reservoirs within a master equation
approach JOSA ; Barton ; Karthik ; Antezza  then the transmission coefficients
for the heat current in both directions are not the same which is a
manifestation of a broken detailed balance SABGSA2023 . However, in this case
the total transferred power between both resonances itself is reciprocal even
in presence of synthetic electric and magnetic fields SABGSA2023 .

That the transferred power is reciprocal might not be surprising for two
reasons. First of all, in the context of Rytov’s fluctuational electrodynamics
it can easily shown that the total radiative heat flux between two objects is
always reciprocal Latella2017 . Non-reciprocal effects necessitate at least a
third object and non-reciprocal material properties of the objects or
environment Herz ; RMP . Another argument is that within the quantum master
equation approach for linearly coupled oscillators typically non-linear
effects need to be included to have nonreciprocal heat flow LandiEtAl2022
even though it seems that nonreciprocal heat flow can also be generated by
specific choices of temperatures in a linear chain of oscillators
PurkayasthaEtAl2016 ; KalantarEtAl2021 . However, as we will show below the
application of synthetic electric and magnetic fields can indeed generate
nonreciprocal heat flow in a tight-binding configuration of four coupled
resonators without the need of non-linearity due to the presence of the
synthetic magnetic field.

Now, we distinguish our work from previous studies. Several kinds of
modulations have been proposed like the periodic modulation of the
permittivity LiraEtAl2012 ; ChamanaraEtAl2017 ; EstepEtAl2014 . Such
modulations have been shown to introduce synthetic magnetic fields for photons
TzuanEtAl2014  and consequently related effects like the Aharonov-Bohm effect
for photons FangEtAl2012 . In the context of thermal radiation it could be
demonstrated that permittivity modulations can introduce nonreciprocity which
manifests itself in a breakdown of the detailed balance in Kirchhoff’s law
GhanekarEtAl2022  and can be employed for photonic refrigeration
BuddhirajuEtAl2020 . In similar approaches a combined dynamical modulation of
the resonances of heat exchanging objects and their interaction strength were
applied resulting in a heat pumping effect and non-reciprocal heat fluxes in a
three resonator configurations LiEtAl2019 ; AlcazarEtAl2021 . Heat pumping
effects also exist when only the interaction strengths in three-body
configurations are dynamically modulated MessinaPBA2020 . It must be
emphasized that these effects are different from the heat shuttling effect
where temperatures or chemical potential of two reservoirs are periodically
modulated around their equilibrium values in order to have a heat transport
despite the fact that the system is in average in equilibrium LiEtAl2008 ;
LiEtAl2009 ; Latella . Indeed, in that case the modulation affects the baths
only and not resonator parameters. Finally, it could be demonstrated
theoretically that adiabatic dynamical modulation of resonators with non-
reciprocal conductances geometrical phases can increase or reduce the thermal
relaxation SABPBABerry  and rapid magnetic field modulations in magneto-
optical systems can substantially increase the cooling Messina .

In this work, we extend the quantum Langevin equation (qLE) and quantum master
equation (qME) approach used in Ref. SABGSA2023  to the case of $N$ coupled
arbitrary resonators with their own heat baths as sketched in Fig. 1 with
applied syntetic electric and magnetic synthetic fields. Both methods can be
used to calculate the heat flux between any two resonators which are coupled
to their own reservoirs. We show numerically that both methods give the same
values for the heat flux. The qLE approach naturally allows for calculating
the heat flux spectra, whereas the master-equation method is a better choice
for fast numerical calculations of the heat flux. We use both methods to show
that the heat flux itself is nonreciprocal in the presence of synthetic fields
in a linear tight-binding chain of four-resonators.

The manuscript is organized as follows: First, in Sec. II we introduce the
standard master-equation for $N$ coupled resonators with $N$ reservoirs. We
derive the dynamical equations for the mean values of products of the
resonator amplitudes and introduce the qLE for the coupled resonator system.
In Sec. III we introduce the synthetic fields in the qLE approach and provide
a formal solution in Fourier space. In Sec. IV we introduce the synthetic
fields in the master-equation approach and give a formal solution by making a
Fourier series ansatz. In Sec. V we show the occurrence of nonreciprocal heat
flux in presence of synthetic electric and magnetic fields in a four-resonator
chain. We conclude with a summary of the main results in Sec. VI.

##  II Langevin and master equations

We start with writing down the Hamiltonian of a coupled harmonic oscillator
system (each oscillator coupled to its own heat bath of oscillators) which is
given by AgarwalBook2012

$H=H_{S}+\sum_{i}H_{B,i}+\sum_{i}H_{SB,i}$ (1)

with the Hamiltonian of the system of coupled oscillators

$H_{S}=\sum_{i}\hbar\omega_{i}a^{\dagger}_{i}a_{i}+\sum_{i,j,i\neq j}\hbar
g_{ij}a^{\dagger}_{i}a_{j},$ (2)

with resonance frequencies $\omega_{i}$ and coupling constants
$g_{ij}=g_{ji}^{*}$ for hermitian system $H_{S}^{\dagger}=H_{S}$ and the
bosonic creation and annihilation operators $a^{\dagger}_{i}$ and $a_{i}$. The
bath oscillator Hamiltonians are given by($i=1,\ldots,N$)

$H_{B,i}=\sum_{j}\hbar\omega_{ij}b_{ij}^{\dagger}b_{ij}$ (3)

with bosonic creation and annihilation operators $b^{\dagger}_{ij}$ and
$b_{ij}$ and the Hamiltonians describing the linear coupling between the
system oscillators and their baths are given by

$H_{SBi}={\rm
i}\hbar\sum_{j}g_{B,ij}(a_{i}+a^{\dagger}_{i})(b_{ij}-b_{ij}^{\dagger})$ (4)

with corresponding coupling constants $g_{B,ij}$. By assuming the validity of
the Born-Markov approximation and tracing out the bath variables one can
arrive at the qME AgarwalBook2012

$\begin{split}\frac{\partial\rho_{S}}{\partial t}&=-{\rm
i}\sum_{i}\omega_{i}[a^{\dagger}_{i}a_{i},\rho_{S}]\\\ &\quad-{\rm
i}\sum_{i,j;i\neq j}g_{ij}[a^{\dagger}_{i}a_{j},\rho_{S}]\\\
&\quad-\sum_{i}\kappa_{i}(n_{i}+1)\bigl{(}a^{\dagger}_{i}a_{i}\rho_{S}-2a_{i}\rho_{S}a^{\dagger}_{i}+\rho_{S}a^{\dagger}_{i}a_{i}\bigr{)}\\\
&\quad-\sum_{i}\kappa_{i}n_{i}\bigl{(}a_{i}a^{\dagger}_{i}\rho_{S}-2a^{\dagger}_{i}\rho_{S}a_{i}+\rho_{S}a_{i}a^{\dagger}_{i}\bigr{)}\end{split}$
(5)

where the coupling to the bath oscillators is formally given in terms of the
coupling constants
$\kappa_{i}=\pi\sum_{j}g_{B,ij}^{2}\delta(\omega_{ij}-\omega_{i})$ and
$n_{i}=[\exp(\hbar\omega_{i}/k_{\rm B}T_{i})-1]^{-1}$ are the mean occupation
numbers at the bath temperatures $T_{i}$. For the sake of generality we assume
that $g_{ij}\neq g_{ji}$ which will allow to include systems without
“inversion symmetry”.

From the qME we can derive the dynamical equation for the mean values of any
observable. For example, for the mean values of products of raising and
lowering operators we obtain the set of equations ($k,l=1,\ldots,N;k\neq l$)

$\displaystyle\frac{{\rm d}}{{\rm d}t}\langle a^{\dagger}_{k}a_{k}\rangle$
$\displaystyle=-{\rm i}\sum_{j,j\neq k}\bigl{(}g_{kj}\langle
a^{\dagger}_{k}a_{j}\rangle-g_{jk}\langle a_{k}a^{\dagger}_{j}\rangle\bigr{)}$
$\displaystyle\quad-2\kappa_{k}\langle
a^{\dagger}_{k}a_{k}\rangle+2\kappa_{k}n_{k},$ (6) $\displaystyle\frac{{\rm
d}}{{\rm d}t}\langle a^{\dagger}_{k}a_{l}\rangle$
$\displaystyle={\Omega}_{kl}\langle a^{\dagger}_{k}a_{l}\rangle-{\rm
i}\sum_{j\neq k;j\neq l}\bigl{(}g_{lj}\langle a^{\dagger}_{k}a_{j}\rangle-
g_{jk}\langle a^{\dagger}_{j}a_{l}\rangle\bigr{)}$ $\displaystyle\quad-{\rm
i}g_{lk}\bigl{(}\langle a_{k}^{\dagger}a_{k}\rangle-\langle
a_{l}^{\dagger}a_{l}\rangle\bigr{)}$ (7)

with

${\Omega}_{kl}={\rm i}(\omega_{k}-\omega_{l})-\kappa_{k}-\kappa_{l}.$ (8)

In the following we will refer to this set of equations for the mean values of
operator products (6) and (7) as master-equation approach as they are derived
from the qME in Eq. (5).

Similarly, one obtains for the time evolution of the mean values of the
raising/lowering operators of each oscillator $a_{i}$ the set of equations
($k=1,\ldots,N$)

$\frac{{\rm d}}{{\rm d}t}\langle a_{k}\rangle=-\Omega_{k}\langle
a_{k}\rangle-{\rm i}\sum_{i;i\neq k}g_{ki}\langle a_{i}\rangle$ (9)

with $\Omega_{k}\equiv{\rm i}\omega_{k}+\kappa_{k}$. The set of equations for
the mean values of the lowering operators of the two oscillators in Eq. (9)
motivates the introduction of a set of qLE for the operators themselves
instead for their expectation values

$\dot{a}_{k}=-{\rm i}\omega_{k}a_{k}-\kappa_{k}a_{k}-{\rm i}\sum_{i,i\neq
k}g_{ki}a_{i}+F_{k},$ (10)

where the coupling to baths is taken into account by the bath operators
$F_{k}$ which obviously must fulfill $\langle F_{k}\rangle=0$ to retrieve Eqs.
(9). Furthermore, to be consistent with the set of Eqs. (6)-(7) the
correlation functions of the bath operators are given by

$\displaystyle\langle F^{\dagger}_{k}(t)F_{k}(t^{\prime})\rangle$
$\displaystyle=2\kappa_{k}n_{k}\delta(t-t^{\prime}),$ (11)
$\displaystyle\langle F_{k}(t)F^{\dagger}_{k}(t^{\prime})\rangle$
$\displaystyle=2\kappa_{k}(n_{k}+1)\delta(t-t^{\prime})$ (12)

and $\langle F_{k}F_{k}\rangle=\langle
F_{k}^{\dagger}F_{k}^{\dagger}\rangle=0$. Furthermore, the bath operators of
different baths are uncorrelated.

##  III Langevin equations with synthetic fields

We now use the set of qLEs as introduced above and include a frequency
modulation ($k=1,\ldots,N$)

$\displaystyle\omega_{k}$
$\displaystyle\rightarrow\omega_{k}+m_{k}\beta\cos(\Omega t+\theta_{k}),$ (13)

with phase shifts $\theta_{k}$ and $m_{k}=\\{0,1\\}$ ($m_{k}=0$ modulation of
oscillator $k$ turned of, $m_{k}=1$ modulation turned on). The set of coupled
qLE in frequency space is therefore ($k=1,\ldots,N$)

$X_{k}a_{k}(\omega)+{\rm i}\sum_{l\neq
k}g_{kl}a_{l}(\omega)=F_{k}+\frac{\beta}{2{\rm i}}\bigl{(}a_{k,-}{\rm
e}^{-{\rm i}\theta_{k}}+a_{k,+}{\rm e}^{+{\rm i}\theta_{k}}\bigr{)}$ (14)

introducing

$X_{k}={\rm i}(\omega_{k}-\omega)+\kappa_{k}$ (15)

and the short-hand notation

$a_{k,\pm}=a_{k}(\omega\pm\Omega).$ (16)

The coupled qLEs can now be brought in matrix form

$\boldsymbol{\psi}=\mathds{M}\mathbf{F}+\frac{\beta}{2{\rm
i}}\mathds{M}\mathds{Q}_{+}\boldsymbol{\psi}_{+}+\frac{\beta}{2{\rm
i}}\mathds{M}\mathds{Q}_{-}\boldsymbol{\psi}_{-}$ (17)

by introducing the vectors

$\boldsymbol{\psi}=\begin{pmatrix}a_{1}(\omega)\\\ \vdots\\\
a_{N}(\omega)\end{pmatrix},\boldsymbol{\psi}_{\pm}=\begin{pmatrix}a_{1}(\omega\pm\Omega)\\\
\vdots\\\
a_{N}(\omega\pm\Omega)\end{pmatrix},\mathbf{F}=\begin{pmatrix}F_{1}(\omega)\\\
\vdots\\\ F_{N}(\omega)\end{pmatrix},$ (18)

and the matrices

$\mathds{M}=\mathds{A}^{-1}\quad\text{with}\quad\mathds{A}=\begin{pmatrix}X_{1}&{\rm
i}g_{12}&\ldots&{\rm i}g_{1N}\\\ {\rm i}g_{21}&X_{2}&\ldots&{\rm i}g_{2N}\\\
\vdots&\dots&\vdots&\vdots\\\ {\rm i}g_{N1}&g_{N2}&\ldots&X_{N}\end{pmatrix}$
(19)

and

$\mathds{Q}_{\pm}={\rm diag}({\rm e}^{\pm{\rm i}\theta_{1}}m_{1},\ldots,{\rm
e}^{\pm{\rm i}\theta_{N}}m_{N}).$ (20)

In Eq. (17) it can be clearly seen that due to the modulation there are
couplings to the next sidebands $\omega\pm\Omega$ so that this set of
equations is recursive and infinitely large. These side bands can be
understood as being the consequent of a synthetic constant electric field.
Furthermore, the phase shift adds a phase $\pm\theta_{k}$ to this coupling
which can be understood as a consequence of a synthetic magnetic field.

The solution of the coupled qLEs (17) can formally be written down for all
orders. By introducing the block vectors

$\displaystyle\uline{\boldsymbol{\psi}}$
$\displaystyle=(\ldots,\boldsymbol{\psi}_{++},\boldsymbol{\psi}_{+},\boldsymbol{\psi},\boldsymbol{\psi}_{-},\boldsymbol{\psi}_{--}\ldots)^{t},$
(21) $\displaystyle\uline{\mathbf{F}}$
$\displaystyle=(\ldots,\mathbf{F}_{++},\mathbf{F}_{+},\mathbf{F},\mathbf{F}_{-},\mathbf{F}_{--},\ldots)^{t},$
(22)

the diagonal block matrix

$\uline{\mathds{M}}=\begin{pmatrix}\ldots&\ldots&\ldots&\ldots&\ldots\\\
\ldots&\mathds{M}_{+}&\mathds{O}&\mathds{O}&\ldots\\\
\ldots&\mathds{O}&\mathds{M}&\mathds{O}&\ldots\\\
\ldots&\mathds{O}&\mathds{O}&\mathds{M}_{-}&\ldots\\\
\ldots&\ldots&\ldots&\ldots&\ldots\end{pmatrix}$ (23)

and tridiagonal block matrix

$\uline{\mathds{L}}=\begin{pmatrix}\ldots&\ldots&\ldots&\ldots&\ldots\\\
\frac{{\rm i}\beta}{2}\mathds{M}_{+}\mathds{Q}_{+}&\mathds{1}&\frac{{\rm
i}\beta}{2}\mathds{M}_{+}\mathds{Q}_{-}&\mathds{O}&\ldots\\\ \ldots&\frac{{\rm
i}\beta}{2}\mathds{M}\mathds{Q}_{+}&\mathds{1}&\frac{{\rm
i}\beta}{2}\mathds{M}\mathds{Q}_{-}&\ldots\\\ \ldots&\mathds{O}&\frac{{\rm
i}\beta}{2}\mathds{M}_{-}\mathds{Q}_{+}&\mathds{1}&\frac{{\rm
i}\beta}{2}\mathds{M}_{-}\mathds{Q}_{-}\\\
\ldots&\ldots&\ldots&\ldots&\ldots\end{pmatrix}$ (24)

we can rewrite the coupled qLE (17) as a matrix equation

$\uline{\mathds{L}}\uline{\boldsymbol{\psi}}=\uline{\mathds{M}}\uline{\mathbf{F}}.$
(25)

Hence

$\uline{\boldsymbol{\psi}}=\uline{\mathds{L}}^{-1}\uline{\mathds{M}}\uline{\mathbf{F}}.$
(26)

By considering only block vectors $\uline{\boldsymbol{\psi}}$ of $2n+1$
vectors $\boldsymbol{\psi}$ with the corresponding block matrices of size
$(2n+1)\times(2n+1)$ submatrices we obtain the perturbation results up to
order $n$. Note, that the full size of the block vectors and matrices is
$N(2n+1)$ and $N^{2}(2n+1)^{2}$, resp.

To evaluate these spectra in our general formalism, we start with Eq. (26) and
introduce the block matrices $\uline{\mathds{Y}}_{1}={\rm
diag}(1,0,0,0,1,0,0,0,\ldots)$, $\uline{\mathds{Y}}_{2}={\rm
diag}(0,1,0,0,0,1,0,0,\ldots)$ etc. so that
$\sum_{k}\uline{\mathds{Y}}_{k}=\uline{\mathds{1}}$. These matrices allow us
to split the contributions from all baths $k$ so that

$\uline{\boldsymbol{\psi}}=\sum_{k=1}^{N}\uline{\mathds{L}}^{-1}\uline{\mathds{M}}\mathds{Y}_{k}\uline{\mathbf{F}}.$
(27)

To evaluate products we use the fluctuation-dissipation theorem in the form

$\langle
F_{k}^{\dagger}(\omega+l\Omega)F_{k^{\prime}}(\omega^{\prime}+l^{\prime}\Omega)\rangle=\delta_{k,k^{\prime}}\delta_{l,l^{\prime}}2\pi\delta(\omega-\omega^{\prime})\langle
F_{k}^{\dagger}F_{k}\rangle_{\omega},$ (28)

where $\langle F_{k}^{\dagger}F_{k}\rangle_{\omega}=2\kappa_{k}n_{k}$. Here in
agreement with the treatment in the qME approach we are assuming that $n_{k}$
is constant as demanded by the assumption of white noise. This assumption is
justified for $\beta\ll\omega_{k}$ and $\Omega\ll k_{\rm B}T/\hbar$. Then we
have

$\langle\uline{\boldsymbol{\psi}}_{\alpha}^{\dagger}\uline{\boldsymbol{\psi}}_{\epsilon}\rangle_{\omega}=\sum_{k=1}^{N}2\kappa_{k}n_{k}\bigl{(}\uline{\mathds{L}}^{-1}\uline{\mathds{M}}\uline{\mathds{Y}}_{k}\uline{\mathds{M}}^{\dagger}{\uline{\mathds{L}}^{-1}}^{\dagger}\bigr{)}_{\epsilon,\alpha}$
(29)

using the properties $\uline{\mathds{Y}}_{k}^{\dagger}=\uline{\mathds{Y}}_{k}$
and $\uline{\mathds{Y}}_{k}\uline{\mathds{Y}}_{k}=\uline{\mathds{Y}}_{k}$.
From this expression we can numerically calculate all spectral correlation
function.

The mean heat flux (transferred power over one oscillation period) from
oscillator $k$ at temperature $T_{k}\neq 0\,{\rm K}$ to an oscillator $l$ at
temperature $T_{l}=0\,{\rm K}$ is defined by LandiEtAl2022

$P_{k\rightarrow l}=\int\frac{{\rm
d}\omega}{2\pi}\hbar\omega_{k}2\kappa_{l}\langle
a^{\dagger}_{l}a_{l}\rangle_{\omega}$ (30)

where $\langle a^{\dagger}_{l}a_{l}\rangle_{\omega}$ is given by
$\langle\uline{\boldsymbol{\psi}}_{\alpha}^{\dagger}\uline{\boldsymbol{\psi}}_{\epsilon}\rangle_{\omega}$
from Eq. (29) with $\epsilon=\alpha=Nn+l$ coming from the term involving
$n_{k}$ due to bath $k$. The total power emitted by the hot oscillator $k$ is
given by

$P_{k}^{\rm em}=\int\frac{{\rm
d}\omega}{2\pi}\hbar\omega_{k}2\kappa_{k}(n_{k}-\langle
a^{\dagger}_{k}a_{k}\rangle_{\omega}).$ (31)

Since energy conservation is fulfilled the sum of the power received by all
resonators $l$ equals that emitted by resonator $k$ so that we have

$P_{k}^{\rm em}=\sum_{l\neq k}P_{k\rightarrow l}.$ (32)

This relation is very useful for validating the numerical results.

##  IV Master equations with synthetic fields

Now, instead of the qLEs we use the qME (6)-(7) with a periodic driving as in
Eq. (13). This directly leads to the set of equations

$\displaystyle\frac{{\rm d}}{{\rm d}t}\langle a^{\dagger}_{k}a_{k}\rangle$
$\displaystyle=-{\rm i}\sum_{j,j\neq k}\bigl{(}g_{kj}\langle
a^{\dagger}_{k}a_{j}\rangle-g_{jk}\langle a_{k}a^{\dagger}_{j}\rangle\bigr{)}$
$\displaystyle\quad-2\kappa_{k}\langle
a^{\dagger}_{k}a_{k}\rangle+2\kappa_{k}n_{k},$ (33) $\displaystyle\frac{{\rm
d}}{{\rm d}t}\langle a^{\dagger}_{k}a_{l}\rangle$
$\displaystyle=\tilde{\Omega}_{kl}\langle a^{\dagger}_{k}a_{l}\rangle-{\rm
i}\sum_{j\neq k;j\neq l}\bigl{(}g_{lj}\langle a^{\dagger}_{k}a_{j}\rangle-
g_{jk}\langle a^{\dagger}_{j}a_{l}\rangle\bigr{)}$ $\displaystyle\quad-{\rm
i}g_{lk}\bigl{(}\langle a_{k}^{\dagger}a_{k}\rangle-\langle
a_{l}^{\dagger}a_{l}\rangle\bigr{)}$ (34)

with

$\begin{split}\tilde{\Omega}_{kl}&=+{\rm
i}(\omega_{k}-\omega_{l})-\kappa_{k}-\kappa_{l}\\\ &\quad+{\rm
i}\beta\bigl{[}\cos(\Omega t+\theta_{k})-\cos(\Omega
t+\theta_{l})\bigr{]}.\end{split}$ (35)

To solve the equations we make the Fourier series ansatz for the expectation
values of each observable $O$ such that

$\langle O\rangle=\sum_{n}{\rm e}^{-{\rm i}n\Omega t}\langle O\rangle_{n}.$
(36)

Then we note that

$\begin{split}\sum_{n}{\rm e}^{-{\rm i}n\Omega t}&\langle
O\rangle_{n}\bigl{(}\cos(\Omega t+\theta_{k})-\cos(\Omega
t+\theta_{l})\bigr{)}\\\ &=\sum_{n}{\rm e}^{-{\rm i}n\Omega
t}\biggl{[}\frac{\eta_{kl}}{2}\langle
O\rangle_{n+1}+\frac{\eta^{*}_{kl}}{2}\langle
O\rangle_{n-1}\biggr{]}\end{split}$ (37)

with

$\eta_{kl}=(m_{k}{\rm e}^{{\rm i}\theta_{k}}-m_{l}{\rm e}^{{\rm
i}\theta_{l}}).$ (38)

Inserting this ansatz in the set of Eqs. (33)-(34) gives the following set of
equations for the Fourier components

$\displaystyle\bigl{[}-{\rm i}n\Omega+2\kappa_{k}\bigr{]}\langle
a^{\dagger}_{k}a_{k}\rangle_{n}$ $\displaystyle=-{\rm i}\sum_{j,j\neq
k}\bigl{(}g_{kj}\langle a^{\dagger}_{k}a_{j}\rangle_{n}$
$\displaystyle\qquad\quad-g_{jk}\langle
a_{k}a^{\dagger}_{j}\rangle_{n}\bigr{)}$
$\displaystyle\quad+2\kappa_{k}{n}_{k}\delta_{n0}$ (39)
$\displaystyle\bigl{[}-{\rm i}n\Omega-\Omega_{kl}\bigr{]}\langle
a^{\dagger}_{k}a_{l}\rangle_{n}$ $\displaystyle=-{\rm i}\sum_{j\neq k;j\neq
l}\bigl{(}g_{lj}\langle a^{\dagger}_{k}a_{j}\rangle_{n}$
$\displaystyle\qquad\quad-g_{jk}\langle
a^{\dagger}_{j}a_{l}\rangle_{n}\bigr{)}$ $\displaystyle\quad-{\rm
i}g_{lk}\bigl{(}\langle a_{k}^{\dagger}a_{k}\rangle_{n}-\langle
a_{l}^{\dagger}a_{l}\rangle_{n}\bigr{)}$ $\displaystyle\quad-\frac{{\rm
i}\beta\eta_{kl}}{2}\langle a^{\dagger}_{k}a_{l}\rangle_{n+1}$
$\displaystyle\quad-\frac{{\rm i}\beta\eta^{*}_{kl}}{2}\langle
a^{\dagger}_{k}a_{l}\rangle_{n-1}.$ (40)

The set of equations for the Fourier components can again be written in matrix
form

$\uuline{\mathds{L}}\uline{\boldsymbol{\psi}}=\uline{\boldsymbol{\kappa}}$
(41)

when introducing the block vector

$\uline{\boldsymbol{\psi}}=(\ldots,\boldsymbol{\psi}_{1},\boldsymbol{\psi}_{0},\boldsymbol{\psi}_{-1},\ldots)^{t}.$
(42)

with

$\begin{split}\boldsymbol{\psi}_{n}&=(\langle
a^{\dagger}_{1}a_{1}\rangle_{n},\ldots,\langle
a^{\dagger}_{N}a_{N}\rangle_{n},\\\ &\quad\langle
a^{\dagger}_{1}a_{2}\rangle_{n},\langle
a^{\dagger}_{2}a_{1}\rangle_{n},\ldots,\langle
a^{\dagger}_{1}a_{N}\rangle_{n},\langle a^{\dagger}_{N}a_{1}\rangle_{n},\\\
&\quad\langle a^{\dagger}_{2}a_{3}\rangle_{n},\langle
a^{\dagger}_{3}a_{2}\rangle_{n},\ldots,\langle
a^{\dagger}_{2}a_{N}\rangle_{n},\langle a^{\dagger}_{N}a_{2}\rangle_{n},\\\
&\quad\ldots\langle a^{\dagger}_{N-1}a_{N}\rangle_{n},\langle
a^{\dagger}_{N}a_{N-1}\rangle)^{t}\end{split}$ (43)

as well as the block vector

$\uline{\boldsymbol{\kappa}}=(\ldots,0,0,+2\kappa_{1}{n}_{1},\ldots,+2\kappa_{N}{n}_{N},0,0,\ldots)^{t}.$
(44)

The block matrix $\uline{\mathds{L}}$ then takes the form of a tri-diagonal
block matrix

$\uline{\mathds{L}}=\begin{pmatrix}\ldots&\ldots&\ldots&\ldots&\ldots\\\
\ldots&\mathds{M}_{1}&\mathds{G}^{-}&\mathds{O}&\ldots\\\
\ldots&\mathds{G}^{+}&\mathds{M}_{0}&\mathds{G}^{-}&\ldots\\\
\ldots&\mathds{O}&\mathds{G}^{+}&\mathds{M}_{-1}&\ldots\\\
\ldots&\ldots&\ldots&\ldots&\ldots\end{pmatrix}$ (45)

with

$\mathds{M}_{n}=\begin{pmatrix}-{\rm i}n\Omega+2\kappa_{1}&0&\ldots&0&+{\rm
i}g_{12}&-{\rm i}g_{21}&\ldots&0&0\\\ 0&-{\rm
i}n\Omega+2\kappa_{2}&\ldots&0&-{\rm i}g_{12}&{\rm
i}g_{21}&\ldots&\ldots&\vdots\\\
\vdots&\ldots&\ldots&\ldots&\ldots&\ldots&\ldots&\ldots&\vdots\\\
0&0&\ldots&-{\rm i}n\Omega+2\kappa_{N}&0&0&\ldots&-{\rm i}g_{N-1,N}&{\rm
i}g_{N,N-1}\\\ -{\rm i}g_{21}&{\rm i}g_{12}&\ldots&0&-{\rm
i}n\Omega-\Omega_{12}&0&\ldots&0&0\\\ {\rm i}g_{21}&-{\rm
i}g_{12}&\ldots&0&0&-{\rm i}n\Omega-\Omega_{21}&\ldots&0&0\\\
\vdots&\ldots&\ldots&\ldots&\ldots&\ldots&\ldots&\ldots&\vdots\\\ 0&0&-{\rm
i}g_{3N}&{\rm i}g_{3N}&\ldots&\ldots&\ldots&0&-{\rm
i}n\Omega-\Omega_{N,N-1}\end{pmatrix}$ (46)

and

$\mathds{G}^{+}=\frac{{\rm i}\beta}{2}{\rm
diag}(0,\ldots,0,\eta_{12},-\eta_{12},\ldots,\eta_{N-1,N},-\eta_{N-1,N})$ (47)

and $\mathds{G}^{-}$ defined as the matrix obtained from $\mathds{G}$ when
complex conjugating $\eta_{kl}$. The different “perturbation orders” $n$ can
be obtained by using $2n+1$ subblocks in the matrix $\uline{\mathds{L}}$.
Note, that even though we use the same notation as in the qLE approach, the
used vectors and matrices are different and also have a different dimension.
Here the dimension of the block vectors and matrices is $N^{2}(2n+1)$ and
$N^{4}(2n+1)^{2}2$.

The mean heat flux (transferred power over one oscillation period) from
oscillator $k$ at temperature $T_{k}$ to an oscillator $l$ at temperature
$T_{l}=0\,{\rm K}$ is defined by SABGSA2023

$P_{k\rightarrow l}=\hbar\omega_{k}2\kappa_{l}\langle
a^{\dagger}_{l}a_{l}\rangle_{0}$ (48)

taking $n_{i}=0$ for all other resonators. Again the total emitted mean power
by particle $k$ is given by

$P_{k}^{\rm em}=\hbar\omega_{k}2\kappa_{k}(n_{k}-\langle
a^{\dagger}_{k}a_{k}\rangle_{0})$ (49)

and we have energy conservation, i.e. $P_{k}^{\rm em}=\sum_{l\neq
k}P_{k\rightarrow l}$. The advantage of the qME approach is that, differently
from the qLE (31)-(30), a frequency integration is not necessary. On the other
hand, the size of the matrices for a given perturbation order is much larger
than for the qLE approach.

##  V Fours resonators case: nonreciprocal heat flux with synthetic fields

We consider here the heat flux in a chain of four resonators as depicted in
Fig. 2. We assume that all resonators are identical and we further assume
reciprocal nearest-neighbour coupling with identical coupling strength $g$ so
that the non-zero coupling constants are
$g_{12}=g_{21}=g_{32}=g_{23}=g_{34}=g_{43}=g$. The resonance frequencies
$\omega_{1}$ and $\omega_{4}$ of the resonator $1$ and $4$ are fixed to
$\omega_{0}$, whereas the resonance frequencies of the resonators in the
middle are modulated as

$\displaystyle\omega_{2}$ $\displaystyle=\omega_{0}+\beta\cos(\Omega t),$ (50)
$\displaystyle\omega_{3}$ $\displaystyle=\omega_{0}+\beta\cos(\Omega
t+\theta).$ (51)

In this configuration, we first determine the power $P_{14}$ transferred from
resonator $1$ to resonator $4$ with $T_{1}=300\,{\rm K}$ and
$T_{2}=T_{3}=T_{4}=0\,{\rm K}$. Then we compare with the heat flow in backward
direction by calculating the power $P_{41}$ transferred from resonator $4$ to
resonator $1$ with $T_{4}=300\,{\rm K}$ and $T_{1}=T_{2}=T_{3}=0\,{\rm K}$.
Hence, only the first and the last resonator are in our configuration coupled
to a heat bath. Therefore the modulation frequency $\Omega$ and the modulation
strength $\beta$ are here in principle not limited by the constraint due to
the white noise assumption because the two resonators in the middle have zero
temperature. Nonetheless, we will restrict ourselves to values which fulfill
the above criteria for the white noise approximation. For our numerical
calculations we use $\omega_{0}=1.69\times 10^{14}\,{\rm rad/s}$ and
$\kappa=0.013\omega_{0}$ which are the the values taken from those for a
graphene flake with $E_{F}=0.4\,{\rm eV}$ from Ref. Graphene . The coupling
constant $g$ is determined by the near-field heat flux value which depends on
the relative distance between the graphene flakes. For a distance $d=100\,{\rm
nm}$ between two graphene flakes a fitting of the resonator model with the
results from fluctuating electrodynamics SABGSA2023  gives $g=0.011\kappa$.
Hence, we are in the weak coupling regime.

Figure 2:  Sketch of a chain of four resonators $1$, $2$, $3$, $4$ with equal
nearest-neighbor couplings $g$ and resonance frequencies $\omega_{0}$. The
oscillators in the middle are modulated with a modulation strength $\beta$ a
relative phase shift $\theta$ resulting in synthetic electric and magnetic
fields.

In Fig. 3(a) we show the results for the transferred power as function of the
modulation strength $\beta$ and for two different values of $\theta$. We also
show the numerical results obtained with the qME method with Eq. (48) and the
qLE approach with Eq. (30). First of all, we can see that both methods provide
the same values for the exchanged power. Furthermore, it can be seen that the
heat flux is clearly nonreciprocal in contrast to the case of two resonators
or two graphene flakes where the heat flux is reciprocal despite the
nonreciprocal spectra SABGSA2023 .

Figure 3: (a) $P_{14}$ (solid line) and $P_{41}$ (dashed line) from the qME
approach (48) at perturbation order $n=15$ normalized to the value
$P_{14}(\beta=0)=P_{41}(\beta=0)=5.88\times 10^{-22}\,{\rm W}$ for
$g=0.011\kappa$ and $\Omega=0.05\omega_{0}$ for $\theta=0.1\pi$ and
$\theta=0.5\pi$. The filled and open symbols are the results for $P_{14}$ and
$P_{41}$ resulting from integration of spectra as in Fig. 5 from qLE approach
according to Eq. (30) at perturbation order $n=10$. (b) Comparison of exact
numerical results (solid lines) for the difference $P_{14}-P_{41}$ normalized
to $P_{14}(\beta=0)=P_{41}(\beta=0)=5.88\times 10^{-22}\,{\rm W}$ with the
corresponding power difference from the approximate expression (dashed lines)
from Eq. (53).  Figure 4:  We show the relative power transmission $E$ defined
in Eq. (54) as function of the dephasing $\theta$ for $\Omega=0.05\omega_{0}$
and different values of modulation strength $\beta$ using the qME approach in
order $n=15$.

As detailed in Ref. PetersonEtAl2019 , for instance, the nonreciprocity in
transmission as sketched in Fig. 2 can be understood in second-order
perturbation theory as an interference of different transmission paths. The
energy at $\omega_{0}$ provided by resonator $1$ can go through the chain in
second order via the upper and lower sideband at $\omega_{0}\pm\Omega$ by two
”scattering events” $\omega_{0}\rightarrow\omega_{0}+\Omega$ and
$\omega_{0}+\Omega\rightarrow\omega_{0}$ or
$\omega_{0}\rightarrow\omega_{0}-\Omega$ and
$\omega_{0}-\Omega\rightarrow\omega_{0}$ as sketched in Fig. 2. Due to the
presence of the synthetic magnetic field a phase is picked up in this process
which is not the same in forward transmission from resonator $1\rightarrow 4$
and backward transmission from resonator $4\rightarrow 1$. This symmetry
breaking of the synthetic magnetic field can be directly understood from Eq.
(14) which shows that upward and downward transitions in the Floquet side-
bands is connected with picking up a positive or negative phase. Hence the
forward and backward transmission along the upper or lower sidebands results
in different phase factors. For a plane wave with frequency $\omega$ being
transmitted through the coupled resonators $2$ and $3$ the difference in the
transmission is explicitely given by PetersonEtAl2019

$\tau_{23}-\tau_{32}=-2{\rm
i}\frac{\beta^{2}}{4}\bigl{[}\tau(\omega+\Omega))-\tau(\omega-\Omega))\bigr{]}\sin(\theta)$
(52)

where $\tau(\omega)$ is the transmission coefficient without modulation. This
transmission coefficient shows that there is a nonreciprocal transmission for
any phase difference $\theta\neq m\pi$ with integer $m$. From this expression
it can be expected that at least in second-order perturbation theory, i.e.
when $\beta$ is sufficiently small, the largest difference can be expected for
$\theta=\pi/2$. For the four resonator configuration depicted in Fig. 2 a
similar expression can be derived using a second-order perturbation theory for
the qME approach as detailed in appendix A. In the weak coupling limit
$g\ll\kappa$ we find for the difference of heat flux in forward and backward
direction

$\begin{split}\frac{P_{14}-P_{41}}{\hbar\omega_{0}ng}&=\beta^{2}\frac{g^{5}}{\kappa^{5}}\biggl{[}\frac{7}{8}\frac{{\rm
Im}(A^{2})}{|A|^{4}}+\frac{\kappa{\rm Im}(A^{3})}{|A|^{6}}\\\
&\qquad-\frac{\kappa^{3}{\rm
Im}(A^{5})}{|A|^{10}}\biggr{]}\sin(\theta)\end{split}$ (53)

where $A=2\kappa-{\rm i}\Omega$ and $n\equiv n_{1}=n_{4}$ is the mean
occupation number of the resonators $1$ in forward or resonator $4$ in
backward direction. In Fig. 3(b) we compare its predictions with the numerical
exact results from Fig. 3(a) clearly showing its validity in the small $\beta$
limit. This expression has a similar structure as Eq. (52) indicating the same
dependence on $\theta$ in the limit of small driving amplitudes $\beta$. To
see this effect, we show in Fig. 4 relative power transmission

$E\equiv\frac{P_{14}-P_{41}}{P_{14}+P_{41}}.$ (54)

It can be seen that indeed for $\beta<0.05\omega_{0}$ the maximum difference
in forward and backward heat flow happens at $\theta=\pm\pi/2$. For larger
modulation strengths higher order effects play a role so that this maximum
shifts to slightly larger or smaller values of the dephasing.

Figure 5:  Spectra for mean power $P_{14,\omega}=2\kappa\hbar\omega_{0}\langle
a_{4}^{\dagger}a_{4}\rangle_{\omega}$ for the forward heat flow and
$P_{41,\omega}=2\kappa\hbar\omega_{0}\langle
a_{1}^{\dagger}a_{1}\rangle_{\omega}$ for the backward heat flow calculated
from the spectra for the mean occupations numbers $\langle
a_{4}^{\dagger}a_{4}\rangle_{\omega}$ and $\langle
a_{1}^{\dagger}a_{1}\rangle_{\omega}$ in Eq. (29). The modulation parameters
are $\Omega=0.05\omega_{0}$, $\beta=0.05\omega_{0}$, and $\theta=\pi/2$ and we
use perturbation order of $n=10$.

Finally, in Fig. 5 the spectra of the power $P_{14,\omega}$ and
$P_{41,\omega}$ obtained with the qLE approach in the forward and backward
direction are shown using $\Omega=0.05\omega_{0}$, $\beta=0.05\omega_{0}$, and
$\theta=\pi/2$. It can be seen that the spectra for the heat flow in forward
and backward dirction are not the same as also found for two graphene flakes
only SABGSA2023 . Furthermore, it can be seen that the side-band contribution
is very small so that the main nonreciprocity stems from frequencies around
the resonance $\omega_{0}$. Integrating these spectra according to Eq. (30)
gives the full transferred power for the forward and backward direction shown
in Fig. 3(a).

Hence, in our four resonator system we clearly find a nonreciprocal heat flow
due to electric and magnetic synthetic fields. Even though our example might
be difficult to realize in practice, it clearly shows that synthetic electric
and magnetic fields can generate a nonreciprocal heat flux. We emphasize that
this result is not limited to near-field heat transfer between graphene flakes
but it is generally valid for any configuration and any heat transfer channel
which can be described by four coupled resonators with synthetic fields.

##  VI Conclusion

To summarize, based on the local qME we have introduced a formalism for a qLE
and qME approach for $N$ coupled resonators with electric and magnetic
synthetic fields. The qLE approach is the natural choice when heat flux
spectra are studied whereas for the heat flow the qME approach is a better
choice, because it is faster. As a very important example, we have used both
approaches to show for a system of four linearly coupled resonators that the
heat flow is nonreciprocal when electric and magnetic synthetic fields are
present. We have also verified numerically that both approach give the same
values for the heat flux. Even though for the numerical evaluation we have
considered the near-field heat transfer in a system of four coupled graphene
flakes our findings are very general and applicable to any system and any heat
flux channel which can be described by coupled resonators. Hence, our
formalism provides the fundament for further studies on the heat flux and
other physical effects in coupled many resonator systems with synthetic
fields.

###### Acknowledgements.

S.-A. B. acknowledges support from Heisenberg Programme of the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under the project No.
461632548. S.-A. B. and P. R.-L. thank the University of Montpellier and the
group Theory of Light-Matter and Quantum Phenomena of the Laboratoire Charles
Coulomb for hospitality during his stay in Montpellier where parts of this
work has been done. S.-A. B., P. R.-L., and M. A. acknowledge QuantUM program
of the University of Montpellier. G.S. A. thanks the kind support of The Air
Force Office of Scientific Research [AFOSR award no. FA9550-20-1-0366] and The
Robert A. Welch Foundation [grant no. A-1943]. P. R.-L. acknowledges support
from AYUDA PUENTE 2022, URJC.

##  Appendix A Perturbation theory for qME approach

In this section we derive the second order expression in Eq. (53). To this
end, we start with Fourier equations for the qME in (41) taking terms with
$n=0,1,-1$. Then we have

$\displaystyle\mathds{M}_{0}\boldsymbol{\psi}_{0}$
$\displaystyle=\boldsymbol{\kappa}-\mathds{G}^{+}\boldsymbol{\psi}_{+1}-\mathds{G}^{-}\boldsymbol{\psi}_{-1},$
(55) $\displaystyle\mathds{M}_{+1}\boldsymbol{\psi}_{+1}$
$\displaystyle=-\mathds{G}^{+}\boldsymbol{\psi}_{2}-\mathds{G}^{-}\boldsymbol{\psi}_{0},$
(56) $\displaystyle\mathds{M}_{-1}\boldsymbol{\psi}_{-1}$
$\displaystyle=-\mathds{G}^{+}\boldsymbol{\psi}_{0}-\mathds{G}^{-}\boldsymbol{\psi}_{-2}.$
(57)

By inserting the expressions for $\boldsymbol{\psi}_{+1/-1}$ into the equation
for $\boldsymbol{\psi}_{0}$ and neclecting terms from $|n|\geq 2$ we arrive at

$\mathds{N}\boldsymbol{\psi}_{0}=\boldsymbol{\kappa}\quad\Rightarrow\quad\boldsymbol{\psi}_{0}=\mathds{N}^{-1}\boldsymbol{\kappa}$
(58)

with

$\mathds{N}=\bigl{[}\mathds{M}_{0}-\mathds{G}^{+}\mathds{M}_{+1}^{-1}\mathds{G}^{-}-\mathds{G}^{-}\mathds{M}_{-1}^{-1}\mathds{G}^{+}\bigr{]}.$
(59)

By defining

$\mathds{G}^{+}=\frac{{\rm
i}\beta}{2}\tilde{\mathds{G}}\quad\text{and}\quad\mathds{G}^{-}=\frac{{\rm
i}\beta}{2}\tilde{\mathds{G}}^{*}$ (60)

with $\tilde{G}={\rm
diag}(0,\ldots,0,\eta_{12},-\eta_{21},\ldots,\eta_{N-1,N},-\eta_{N-1,N})$ we
have

$\mathds{N}=\biggl{[}\mathds{M}_{0}+\frac{\beta^{2}}{4}\bigl{(}\tilde{\mathds{G}}\mathds{M}_{+1}^{-1}\tilde{\mathds{G}}^{*}+\tilde{\mathds{G}}^{*}\mathds{M}_{-1}^{-1}\tilde{\mathds{G}}\bigr{)}\biggr{]}.$
(61)

From this expressions it becomes more obvious that the first non-vanishing
contributions to the zeroth order are stemming from the second-order terms,
i.e. there is no contribution linear in $\beta$.

Figure 6:  Comparison of exact numerical results for $P_{14}$ (black lines)
with second order perturpation approach from Eq. (A8) (PA 1) and with those
from Eq. (A9) (PA 2) using the same parameters as in Fig. 3(a) and
$\theta=\pi/2$. The approximations for $P_{41}$ are similar (not shown).

For the tight-binding model of the four identical resonators the involved
vectors have 16 components and the matrices have a size of $16\times 16$. By
definition of $\boldsymbol{\psi}_{0}$ we are interested in the terms
$\mathds{N}^{-1}_{14}$ and $\mathds{N}^{-1}_{41}$ which are determining the
transferred power $P_{4\rightarrow 1}$ and $P_{1\rightarrow 4}$. Obviously,
there can only be a non-reciprocity if
$\mathds{N}^{-1}\neq{\mathds{N}^{-1}}^{t}$. From the equation for $\mathds{N}$
it can be seen that due to the phase terms $\tilde{\mathds{G}}$ and
$\tilde{\mathds{G}}^{*}$ in the second-order contribution, in general, we have
$\mathds{N}\neq\mathds{N}^{t}$ so that also $P_{4\rightarrow 1}\neq
P_{1\rightarrow 4}$ in general. Hence, the synthetic magnetic field results in
an asymmetry for $\mathds{N}$ and hence for $\mathds{N}^{-1}$.

For small $\beta$ we can further simplify the inverse of $\mathds{N}$ in Eq.
(A8)

$\begin{split}\mathds{N}^{-1}&=\biggl{[}\mathds{M}_{0}+\frac{\beta^{2}}{4}\bigl{(}\tilde{\mathds{G}}\mathds{M}_{+1}^{-1}\tilde{\mathds{G}}^{*}+\tilde{\mathds{G}}^{*}\mathds{M}_{-1}^{-1}\tilde{\mathds{G}}\bigr{)}\biggr{]}^{-1}\\\
&=\biggl{[}\mathds{1}+\frac{\beta^{2}}{4}\mathds{M}_{0}^{-1}\bigl{(}\tilde{\mathds{G}}\mathds{M}_{+1}^{-1}\tilde{\mathds{G}}^{*}+\tilde{\mathds{G}}^{*}\mathds{M}_{-1}^{-1}\tilde{\mathds{G}}\bigr{)}\biggr{]}^{-1}\mathds{M}_{0}^{-1}\\\
&\approx\biggl{[}\mathds{1}-\frac{\beta^{2}}{4}\mathds{M}_{0}^{-1}\bigl{(}\tilde{\mathds{G}}\mathds{M}_{+1}^{-1}\tilde{\mathds{G}}^{*}+\tilde{\mathds{G}}^{*}\mathds{M}_{-1}^{-1}\tilde{\mathds{G}}\bigr{)}\biggr{]}\mathds{M}_{0}^{-1}.\end{split}$
(62)

In Fig. 6 we show a comparison of the second-order results using Eq. (A8) and
Eq. (A9) with the numrical exact results. As expected the second-order
expansion is only reliable for small enough values of $\beta$ and the
perturbation expression in Eq. (A8) is valid for a larger range than the
perturbative expression in Eq. (A9).

Now, we want to derive an analytical expression for the heat flux difference.
Note, that the heat flux difference for the forward and backward case is in
our example given by

$P_{14}-P_{41}=4\hbar\omega_{0}n\kappa^{2}\Delta N_{14}$ (63)

where $\Delta N_{14}=\mathds{N}_{14}^{-1}-\mathds{N}_{41}^{-1}$ and $n\equiv
n_{1}=n_{4}$. That means we can focus on $\Delta N_{14}$ and add the
prefactors later. Starting with the approximate expression in Eq. (A9) and
making a Taylor expansion for $g\ll\kappa$ we obtain with Mathematica for
$\Delta N_{14}$ the relatively long expression

$\begin{split}\Delta
N_{14}\approx\frac{\beta^{2}g^{2}}{8|A_{1}|^{6}}\frac{g^{4}}{\kappa^{4}}&\biggl{[}\frac{|A_{1}|^{2}{\rm
Im}(A_{1}^{2})}{A_{0}^{3}}\biggl{(}4\bigl{[}{\rm
Im}(\eta_{13}\eta_{12}^{*})+{\rm Im}(\eta_{34}\eta_{24}^{*})]+3\bigl{[}{\rm
Im}(\eta_{23}\eta_{13}^{*})+{\rm Im}(\eta_{24}\eta_{23}^{*})\bigr{]}+{\rm
Im}(\eta_{14}\eta_{13}^{*})+{\rm Im}(\eta_{24}\eta_{14}^{*})\biggr{)}\\\
&+\frac{{\rm Im}(A_{1}^{3})}{A_{0}^{2}}\biggl{(}{\rm
Im}(\eta_{14}\eta_{12}^{*})+2{\rm Im}(\eta_{24}\eta_{13}^{*})+{\rm
Im}(\eta_{34}\eta_{14}^{*})-3{\rm Im}(\eta_{12}\eta_{23}^{*})-3{\rm
Im}(\eta_{23}\eta_{34}^{*})\biggr{)}\\\ &+2\frac{{\rm
Im}(A_{1}^{4})}{|A_{1}|^{2}A_{0}}\biggl{(}{\rm
Im}(\eta_{24}\eta_{12}^{*})+{\rm
Im}(\eta_{34}\eta_{13}^{*})\biggr{)}-\frac{2{\rm
Im}(A_{1}^{5})}{|A_{1}|^{4}}{\rm
Im}(\eta_{12}\eta_{34}^{*})\biggr{]}\end{split}$ (64)

where we have introduced $A_{n}=2\kappa-{\rm i}n\Omega$. From this expression
it can be seen that only for complex $\eta_{ij}$ there is a non-reciprocity.
It can be further observed that there seem to be plenty of combinations which
give a non-reciprocal heat flux. In our four oscillator example the resonator
$3$ is the only one with a nonzero phase $\theta\equiv\theta_{3}\neq 0$ and
resonator $1$ and $4$ are not modulated at all so that $\eta_{12}=-1$,
$\eta_{14}=0$, $\eta_{24}=1$ and $\eta_{34}={\rm e}^{{\rm
i}\theta}=-\eta_{13}$ and $\eta_{23}=1-{\rm e}^{{\rm i}\theta}$. With these
specific values we get

$\begin{split}\Delta
N_{14}\approx\frac{\beta^{2}g^{6}}{4\kappa^{4}}\sin(\theta)&\biggl{[}\frac{7{\rm
Im}(A_{1}^{2})}{|A_{1}|^{4}A_{0}^{3}}+\frac{4{\rm
Im}(A_{1}^{3})}{|A_{1}|^{6}A_{0}^{2}}-\frac{{\rm
Im}(A_{1}^{5})}{|A_{1}|^{10}}\biggr{]}.\end{split}$ (65)

By adding the correspondig factors as defined in Eq. (63) and realizing that
$A_{0}=2\kappa$ we obtain the approximative analytical expression for the heat
flux difference in Eq. (53).

## References

  * (1) L. Hu, A. Narayanaswamy, X. Chen, G. Chen, Appl. Phys. Lett. 92, 133106 (2008).
  * (2) R. S. Ottens, V. Quetschke, S. Wise, A. A. Alemi, R. Lundock, G. Mueller, D. H. Reitze, D. B. Tanner, B. F. Whiting, Phys. Rev. Lett. 107, 014301 (2011).
  * (3) T. Kralik, P. Hanzelka, M. Zobac, V. Musilova, T. Fort, M. Horak, Phys. Rev. Lett. 109, 224302 (2012).
  * (4) M. Lim, S. S. Lee, B. J. Lee, Phys. Rev. B 91, 195136 (2015).
  * (5) J. I. Watjen, B. Zhao, Z. M. Zhang, Appl. Phys. Lett. 109, 203112 (2016).
  * (6) M. P. Bernardi, D. Milovich. M. Francoeur, Nature Comm. 7, 12900 (2016).
  * (7) B. Song, D. Thompson, A. Fiorino, Y. Ganjeh, P. Reddy, E. Meyhofer, Nature Nanotechn. 11, 509 (2016).
  * (8) A. Fiorino, D. Thompson, L. Zhu, B. Song, P. Reddy, E. Meyhofer, Nano Letters 18, 3711 (2018).
  * (9) R. Mittapally, J. W. Lim, L. Zhang, O. D. Miller, P. Reddy, and E. Meyhofer, Nano Letters 23, 2187 (2023).
  * (10) K. Ito, K. Nishikawa, A. Miura, H. Toshiyoshi, and H. Iizuka, Nano Lett. 17, 4347 (2017).
  * (11) A. Fiorino, D. Thompson, L. Zhu, R. Mittapally, S.-A. Biehs, O. Bezencenet, N. El-Bondry, S. Bansropun, P. Ben-Abdallah, E. Meyhofer, P. Reddy, ACS Nano 12, 5774 (2018).
  * (12) K. Ito, K. Nishikawa, and H. Iizuka, Appl. Phys. Lett. 108, 108, 053507 (2016).
  * (13) J. Kou and A. J. Minnich, Opt. Exp. 26, A730 (2018).
  * (14) N. H. Thomas, M. C. Sherrott, J. Broulliet, H. A. Atwater, and A. J. Minnich, Nano Lett. 19, 3898 (2019).
  * (15) K. Shi, Z. Chen, Y. Xing, J. Yang, X. Xu, J. S. Evans, and S. He, Nano Lett. 22 22, 7753 (2022).
  * (16) I. Latella and P. Ben-Abdallah, Phys. Rev. Lett. 118, 173902, (2017).
  * (17) R. M. Abraham Ekeroth, P. Ben-Abdallah, J.C. Cuevas, and A. Garcia Martin, ACS Photonics 5, 705 (2017).
  * (18) M.-J. He, H. Qi, Y.-X Su, Y.-T. Ren, Y.-J. Zhao, and M. Antezza, Appl. Phys. Lett. 117, 113104 (2020).
  * (19) A. Ott, R. Messina, P. Ben-Abdallah, S.-A. Biehs, J. Photon. Energy 9, 032711 (2019).
  * (20) E. Moncada-Villa and J. C. Cuevas, Phys. Rev. B 103, 075432 (2021).
  * (21) E. Moncada-Villa and J. C. Cuevas, Phys. Rev. B 106, 235430 (2022).
  * (22) L. Ge, K. Gong, Y. Cang, Y. Lue, X. Shi, and Y. Wu, Phys. Rev. B 100, 035414 (2019).
  * (23) H. Wu, Y. Huang, L. Cui, and K. Zhu, Phys. Rev. Appl. 11, 054020 (2019).
  * (24) L. Zhu and S. Fan, Phys. Rev. Lett. 117, 134303 (2016).
  * (25) A. Ott, P. Ben-Abdallah, and S.-A. Biehs, Phys. Rev. B 97, 205414 (2018).
  * (26) M. G. Silveirinha, Phys. Rev. B 95, 115103 (2017).
  * (27) C. Khandekar, Z. Jacob, New J. Phys. 21, 103030 (2019).
  * (28) P. Ben-Abdallah, Phys. Rev. Lett. 116, 084301, (2016).
  * (29) A. Ott, S.-A. Biehs, P. Ben-Abdallah, Phys. Rev. B 101, 241411(R) (2020).
  * (30) P. Doyeux, S. A. Hassani Gangaraj, G. W. Hanson, and M. Antezza, Phys. Rev. Lett. 119, 173901 (2017).
  * (31) A. Ott, R. Messina, P. Ben-Abdallah and S.-A. Biehs, Appl. Phys. Lett. 114, 163105 (2019).
  * (32) Y. Hu, H. Liu, B. Yang, K. Shi, M. Antezza, X. Wu, and Y. Sun, Phys. Rev. Mat, 7, 035201 (2023).
  * (33) L. Fan , Y. Guo, G. T. Papadakis, B. Zhao, Z. Zhao, S. Buddhiraju, M. Orenstein, and S. Fan, Phys. Rev. B 101,085407(2020).
  * (34) J. Dong, W. Zhang, and L. Liu, Appl. Phys. Lett. 119, 021104 (2021).
  * (35) H. Lira, Z. Yu, S. Fan, and M. Lipson, Phys. Rev. Lett. 109, 033901 (2012).
  * (36) N. Chamanara, S. Taravati, Z.-L. Deck-Léger, and C. Caloz, Phys. Rev. B 96, 155409 (2017).
  * (37) N. A. Estep, D. L. Sounas, J. Soric and A. Alu, Nat. Phot. 10, 923 (2014).
  * (38) L. D. Tzuang, K. Fang, P. Nussenzveig, S. Fan, and M. Lipson, Nat. Photonics 8, 701 (2014).
  * (39) K. Fang, Z. Yu, and S. Fan, Phys. Rev. Lett. 108, 153901 (2012).
  * (40) A. Ghanekar, J. Wang, S. Fan, and M. L. Povinelli, ACS Photonics 9, 1157 (2022).
  * (41) S. Buddhiraju, W. Li ,and S. Fan, Phys. Rev. Lett. 124, 077402 (2020).
  * (42) H. Li, L. J. Fernández-Alcázar, F. Ellis, B. Shapiro, T. Kottos, Phys. Rev. Lett. 123, 165901 (2019).
  * (43) L. J. Fernández-Alcázar , R. Kononchuk , H. Li, T. Kottos, Phys. Rev. Lett. 126, 204101 (2021).
  * (44) R. Messina and P. Ben-Abdallah, Phys. Rev. B 101, 165435 (2020).
  * (45) N. Li, P. Hänggi, and B. Li, EPL 84, 40009(2008).
  * (46) N. Li, F. Zhan, P. Hänggi, and B. Li, Phys. Rev. E 80, 011125(2009)
  * (47) I. Latella, R. Messina, J. M. Rubi, and P. Ben-Abdallah, Phys. Rev. Lett. 121, 023903 (2018).
  * (48) S.-A. Biehs and P. Ben-Abdallah, Phys. Rev. B 106, 235412 (2022).
  * (49) R. Messina, A. Ott, C. Kathmann, S.-A. Biehs, and P. Ben-Abdallah, Phys. Rev. B 103, 115440 (2021).
  * (50) E. Lustig and M. Segev, Adv. Opt. Phot. 13, 426 (2021).
  * (51) C. W. Peterson, W. A. Benalcazar, M. Lin, T. L. Hughes, and G. Bahl, Phys. Rev. Lett. 123, 063901 (2019).
  * (52) S.-A. Biehs and G. S. Agarwal, Phys. Rev. Lett. 130, 110401 (2023).
  * (53) S.-A. Biehs and G. S. Agarwal, J. Opt. Soc. Am. B 30, 700 (2013).
  * (54) G. Barton, J. Stat. Phys. 165, 1153 (2016).
  * (55) K. Sasihithlu and G.S. Agarwal, Nanophotonics 7, 1581(2018).
  * (56) G. De Chiara, G. Landi, A. Hewgill, B. Reid, A. Ferraro, A. J. Roncaglia, M. Antezza, New J. Phys. 20, 113024 (2018).
  * (57) F. Herz and S.-A. Biehs, EPL 127, 44001 (2019).
  * (58) S.-A. Biehs, R. Messina, P. S. Venkataram, A. W. Rodriguez, J. C. Cuevas, and P. Ben-Abdallah, Rev. Mod. Phys. 93, 025009 (2021).
  * (59) G. T. Landi, D. Poletti, G. Schaller, Rev. Mod. Phys. 94, 045006 (2022).
  * (60) N. Kalantar, B. K. Agarwalla, D. Segal, Phys. Rev. E 103, 052130 (2021).
  * (61) A. Purkayastha, A. Dhar, and M. Kulkarn, Phys. Rev. A 94, 052134 (2016).
  * (62) S. Thongrattanasiri, F. H.L. Koppens, F. J. G. de Abajo, Phys. Rev. Lett. 108, 047401 (2012).
  * (63) J. H. Shirley, Phys. Rev. 138, B979 (1965).
  * (64) G. S. Agarwal, Quantum Optics, (Cambridge Univ. Press, Cambridge, 2012) .