# NNQS-Transformer: an Efficient and Scalable Neural Network
Quantum States Approach for Ab initio Quantum Chemistry

Yangjun Wu  Institute of Computing Technology, Chinese Academy of Sciences,
UCASBeijingChina<EMAIL_ADDRESS>,  Chu Guo  Hunan Normal
UniversityChangshaChina<EMAIL_ADDRESS>,  Yi Fan  University of Science
and Technology of ChinaHefeiChina<EMAIL_ADDRESS>,  Pengyu Zhou
Institute of Computing Technology, Chinese Academy of Sciences,
UCASBeijingChina<EMAIL_ADDRESS>and  Honghui Shang  University
of Science and Technology of ChinaHefeiChina<EMAIL_ADDRESS>
(2023)

###### Abstract.

Neural network quantum state (NNQS) has emerged as a promising candidate for
quantum many-body problems, but its practical applications are often hindered
by the high cost of sampling and local energy calculation. We develop a high-
performance NNQS method for ab initio electronic structure calculations. The
major innovations include: (1) A transformer based architecture as the quantum
wave function ansatz; (2) A data-centric parallelization scheme for the
variational Monte Carlo (VMC) algorithm which preserves data locality and well
adapts for different computing architectures; (3) A parallel batch sampling
strategy which reduces the sampling cost and achieves good load balance; (4) A
parallel local energy evaluation scheme which is both memory and
computationally efficient; (5) Study of real chemical systems demonstrates
both the superior accuracy of our method compared to state-of-the-art and the
strong and weak scalability for large molecular systems with up to $120$ spin
orbitals.

Quantum chemistry, many-body Schr$\ddot{\text{o}}$dinger equation, neural
network quantum state, transformer based architecture, autoregressive sampling

††journalyear: 2023††copyright: acmlicensed††conference: The International
Conference for High Performance Computing, Networking, Storage and Analysis;
November 12–17, 2023; Denver, CO, USA††booktitle: The International Conference
for High Performance Computing, Networking, Storage and Analysis (SC ’23),
November 12–17, 2023, Denver, CO, USA††price: 15.00††doi:
10.1145/3581784.3607061††isbn: 979-8-4007-0109-2/23/11

##  1\. Introduction

Electronic structure calculation based on first-principle quantum mechanics is
an elementary tool for predicting the chemical and physical properties of
matter. The grand challenge of ab initio electronic structure calculations is
to solve the Schr$\ddot{\text{o}}$dinger equation $H\psi=E\psi$ of an
interacting many-body system of $n$ electrons and $m$ nuclei, whose
Hamiltonian $H$ can be written as (under the Born-Oppenheimer approximation to
freeze the nuclear degrees of freedom),

(1) $\displaystyle
H=-\displaystyle\sum_{i=1}^{n}\frac{1}{2}\nabla_{i}^{2}-\displaystyle\sum_{i=1}^{n}\displaystyle\sum_{A=1}^{m}\frac{Z_{A}}{r_{iA}}+\displaystyle\sum_{i=1}^{n}\displaystyle\sum_{j>i}^{n}\frac{1}{r_{ij}},$

where $r_{ij}=|r_{i}-r_{j}|$ denotes the distance between two electrons at
positions $i$ and $j$, $r_{iA}$ denotes the distance between an electron at
position $i$ and a nucleus at position $A$, and $Z_{A}$ is the atomic number
of the nucleus at position $A$. The state of interest is the ground state
$|\psi\rangle$ that minimizes the energy: $E=\langle\psi|H|\psi\rangle$.

However, the exact solution for the ground state is generally impossible since
the dimension of the state space for wave functions grows exponentially with
the number of electrons. As a result, the full configuration interaction (FCI)
method, which takes into account all the possible states, is currently limited
within $24$ electrons and $24$ orbitals even with the best supercomputer
(VogiatzisJong2017, ). In order to deal with this problem, various
approximation techniques have been developed, such as the Hartree-Fock (HF)
theory and density functional theory (DFT), which are able to scale up to
thousands (or even more) of electrons by reducing the problem to a single-
electron problem instead, but they could introduce errors in the calculation
of physical properties of the system, and may be inaccurate for certain types
of systems. The truncated configuration interaction (CI) (DAVIDSHERRILL1999, )
considers only excitations above the HF reference state up to a fixed order,
while the coupled cluster (CC) methods consider certain nonlinear combinations
of the excitation operators up to a fixed order. In particular, the CC methods
could often give reasonable solutions with a limited computational cost, thus
referred to as the gold standard. However, they could still fail in presence
of strong static correlations (BulikScuseria2015, ). The density matrix
renormalization group (DMRG) method (White1992, ; White1993, ) could give
solutions with comparable precision to FCI and has been used to study
molecular systems with up to $30$ electrons ($108$ orbitals) (BrabecVeis2020,
; LarssonChan2022, ), however it is mainly limited by the expressivity of the
underlying matrix product state ansatz which can only efficiently parameterize
finitely correlated states (GarciaCirac2007, ). DMRG is also very unfriendly
for large-scale parallelization. The conventional quantum Monte Carlo (QMC)
methods (FoulkesRajagopal2001, ; NeedsRios2010, ; AustinLester2012, ),
including the variational Monte Carlo (VMC) and diffusion Monte Carlo (DMC),
can directly deal with real space calculations by using an anti-symmetric wave
function ansatz which accounts for the Pauli exclusion principle of electrons
(the first quantized formalism). However high precision calculations using
first quantized QMC methods in real space are still quite limited due to the
high computational cost. Prominent approaches in this direction include the
FermiNet (PfauFoulkes2020, ) and PauliNet (HermannNoe2020, ) which are
currently limited within $30$ electrons.

In 2017, Carleo and Troyer propose the neural network quantum state (NNQS)
algorithm for many-spin systems, which parameterizes the wave function as a
neural network and uses the VMC algorithm to optimize the parameters
(CarleoTroyer2017, ). In particular they demonstrate that using a very simple
neural network ansatz, namely the restricted Boltzmann machine (RBM) that is
essentially a dense layer plus a nonlinear activation, could already achieve
competitive accuracy to existing tensor network algorithms. Till now the
performance of NNQS has been demonstrated in a wide variety of many-spin
problems (CarleoTroyer2017, ; ChooCarleo2019, ; SharirShashua2020, ;
SchmittHeyl2020, ; YuanDeng2021, ; ZhaoLiang2022, ) as well as many-fermion
problems (MorenoStokes2022, ). NNQS is a promising alternative for quantum
many-body problems since 1) there is no fundamental limit on the expressive
power of (deep) neural network ansatz for many-body quantum state
(DengSarma2017, ; GlasserCirac2018, ; SharirCarleo2022, ; GaoDuan2017, ); 2)
the computational cost scales polynomially in general and 3) it often allows
transparent large-scale parallelization as a general feature of Monte Carlo
methods. In 2020, RBM was first introduced for electronic structure
calculations by converting the fermionic problem in the second quantized
formalism into a many-spin problem using Jordan-Wigner (JW) transformation
(ChooCarleo2020, ), where it is demonstrated that RBM could achieve higher
accuracy than CC with excitations up to second order (CCSD) and third order
(CCSD(T)) for a wide variety of molecular systems within $24$ spin orbitals
(qubits). Two later works pushed the simulation scale to $30$ qubits
(BarrettLvovsky2022, ) and $76$ qubits (Zhao_2023, ) respectively by using
autoregressive neural network ansatz to reduce the sampling cost.

Despite the polynomial scaling of the computational cost of NNQS, for
electronic structure calculations, it is still extremely difficult to scale up
the application size to $100$ qubits or more (while maintaining comparable or
higher accuracy than CC or DMRG) up to date. The challenge is mainly
threefold: 1) For conventional neural networks such as RBM or the
convolutional neural network (CNN), a Markov chain (MC) sampling algorithm is
used to generate samples for energy and gradient estimation, which could be
very inefficient in certain situations with low acception rate. Moreover, for
these neural networks one also has to use the stochastic reconfiguration (SR)
technique (SorellaCapriotti2000, ; SorellaRocca2007, ) for stable convergence
to global minimum, for which one needs to (approximately) compute the inverse
of the $M\times M$ SR matrix for a neural network with $M$ parameters, thus
greatly prohibiting the usage of very deep neural networks as well as the
scalability to a large number of processes. As an example a recent application
of NNQS on the new Sunway supercomputer used a CNN with only around $50,000$
parameters (ZhaoLiang2022, ); 2) The chemical bonding in molecular systems can
be very different in nature and of very different strengths, while for
practical applications it is important to use a model which is easily
extensible to deep architectures for better expressivity and thus can be
versatilely used for different molecular systems; 3) For molecular systems in
the second quantized formalism, the number of Pauli strings after the JW
transformation often scales as $O(N^{4})$ for $N$ qubits, therefore in large
applications, evaluation of the local energy in VMC could be very time-
consuming. For example 3) is the major limiting factor for Ref.
(BarrettLvovsky2022, ) to go beyond $30$ qubits.

In view of these challenges, we believe a highly efficient and scalable NNQS
method is in need, which is an important step towards large-scale applications
in electronic structure calculations and also acts as a powerful tool for
cross validation with existing methods.

The major innovations of this work can be summarized as follows.

  * •

A deep neural network wave function ansatz (QiankunNet (shang2023solving, ))
is used for VMC optimization, which consists of a transformer made of stacked
decoders for the amplitude and a multilevel perceptron (MLP) for the phase,
inspired by the enormous success of transformer based architectures in natural
language processing (NLP) and other classical machine learning tasks
(Transformer2017, ; radford2018gpt, ; radford2019gpt2, ).

  * •

A data-centric parallelization scheme based on MPI communication for VMC
algorithms using autoregressive wave function ansatz, which maximally
preserves the data locality and well adapts to different computing
architectures. The total amount of data communication in this scheme only
scales as $O(N_{u})$ for $N_{u}$ the total number of unique samples and $O(M)$
for $M$ parameters.

  * •

A parallel batch sampling strategy which makes use of the autoregressive
property of our wave function ansatz and achieves low sampling cost and good
load balance.

  * •

A local energy calculation scheme parallelized over batches of samples,
powered by a highly compressed data structure for the Hamiltonian together
with a fused design of nonzero Hamiltonian entry evaluation and local energy
calculation to reduce memory cost, and a sample-aware evaluation scheme
together with an efficient storage of the unique samples as a lookup table to
reduce the computational cost.

  * •

Study of real chemical systems shows that our method could reach superior
accuracy compared to state-of-the-art NNQS methods, and achieve both strong
and weak scaling for the benzene molecular system with $120$ qubits.

Our parallelized implementation of the transformer based NNQS can be easily
transplanted to different high-performance computing architectures as long as
the pytorch library (pytorch, ) is supported and the local energy evaluation
function is implemented on each process. Our work opens a verge of applying
state-of-the-art transformer based architectures for large-scale electronic
structure calculations.

##  2\. Background

In this section, we review the basic ideas of the VMC algorithm,
autoregressive neural networks, typical quantum chemistry Hamiltonians, as
well as state-of-the-art NNQS approaches for electronic structure
calculations.

###  2.1. Variational Monte Carlo algorithm

Figure 1. Schematic of the workflow for the serial NNQS calculation.

Assuming that one has chosen a variational ansatz for the ground state denoted
as $|\psi_{\vec{\theta}}\rangle$, with $\vec{\theta}$ the parameters to be
optimized, then the energy of the system can be written as a function of
$\vec{\theta}$:

(2) $\displaystyle
E(\vec{\theta})=\frac{\langle\psi_{\vec{\theta}}|H|\psi_{\vec{\theta}}\rangle}{\langle\psi_{\vec{\theta}}|\psi_{\vec{\theta}}\rangle},$

where the denominator indicates that the ansatz $|\psi_{\vec{\theta}}\rangle$
may not be normalized in general. The ground state can be approached by
minimizing $E(\vec{\theta})$ against $\vec{\theta}$. Expanding Eq.(2) in the
computational basis, we obtain

$\displaystyle E(\vec{\theta})=$
$\displaystyle\frac{\sum_{\mathbb{x},\mathbb{x}^{\prime}}\langle\psi_{\vec{\theta}}|\mathbb{x}\rangle\langle\mathbb{x}|H|\mathbb{x}^{\prime}\rangle\langle\mathbb{x}^{\prime}|\psi_{\vec{\theta}}\rangle}{\sum_{\mathbb{y}}\langle\psi_{\vec{\theta}}|\mathbb{y}\rangle\langle\mathbb{y}|\psi_{\vec{\theta}}\rangle}$
(3) $\displaystyle=$
$\displaystyle\frac{\sum_{\mathbb{x}}\left(\sum_{\mathbb{x}^{\prime}}H_{\mathbb{x}\mathbb{x}^{\prime}}\psi_{\vec{\theta}}(\mathbb{x}^{\prime})/\psi_{\vec{\theta}}(\mathbb{x})\right)|\psi_{\vec{\theta}}(\mathbb{x})|^{2}}{\sum_{\mathbb{y}}|\psi_{\vec{\theta}}(\mathbb{y})|^{2}},$

where $\mathbb{x}$, $\mathbb{x}^{\prime}$, $\mathbb{y}$ denote specific
bitstrings. We have written
$H_{\mathbb{x}\mathbb{x}^{\prime}}=\langle\mathbb{x}|H|\mathbb{x}^{\prime}\rangle$
and
$\psi_{\vec{\theta}}(\mathbb{x})=\langle\mathbb{x}|\psi_{\vec{\theta}}\rangle$
as the probability amplitude of the wave function ansatz
$|\psi_{\vec{\theta}}\rangle$ in basis $|\mathbb{x}\rangle$ in the second line
of Eq.(2.1). Defining the local energy $E_{loc}(\mathbb{x})$ as

(4) $\displaystyle
E_{loc}(\mathbb{x})=\sum_{\mathbb{x}^{\prime}}H_{\mathbb{x}\mathbb{x}^{\prime}}\psi_{\vec{\theta}}(\mathbb{x}^{\prime})/\psi_{\vec{\theta}}(\mathbb{x}),$

and the (unnormalized) probability
$p_{\vec{\theta}}(\mathbb{x})=|\psi_{\vec{\theta}}(\mathbb{x})|^{2}$, Eq.(2.1)
can be rewritten as

(5) $\displaystyle
E(\vec{\theta})=\frac{\sum_{\mathbb{x}}E_{loc}(\mathbb{x})p_{\vec{\theta}}(\mathbb{x})}{\sum_{\mathbb{y}}p_{\vec{\theta}}(\mathbb{y})}=\mathbb{E}_{p}[E_{loc}(\mathbb{x})].$

It is impossible to exactly evaluate Eq.(5) in general as there exists an
exponential number of different bitstrings. However, one could approximately
evaluate Eq.(5) by sampling from the probability distribution
$p_{\vec{\theta}}(\mathbb{x})$ to obtain a set of $N_{s}$ samples, denoted as
$\\{\mathbb{x}^{1},\mathbb{x}^{2},\dots,\mathbb{x}^{N_{s}}\\}$, and then
averaging over them

(6)
$\displaystyle\tilde{E}(\vec{\theta})=\frac{1}{N_{s}}\sum_{i=1}^{N_{s}}E_{loc}(\mathbb{x}^{i}),$

where we have used $\tilde{E}(\vec{\theta})$ instead of $E(\vec{\theta})$ to
stress that the former is only an approximation of the latter. To this end, we
can see that as long as one can efficiently sample from
$p_{\vec{\theta}}(\mathbb{x})$ (which is the case if we can evaluate
$\psi_{\vec{\theta}}(\mathbb{x})$ efficiently for each $\mathbb{x}$), and
efficiently find those nonzero entry $H_{\mathbb{x}\mathbb{x}^{\prime}}$
together with $\mathbb{x}^{\prime}$, then Eq.(6) can be efficiently evaluated.
To speed up the calculation, it is often advantageous to use a gradient-based
optimizer instead of gradient-free ones. With the obtained samples, one could
compute the approximate gradient of Eq.(5) using automatic differentiation
(BarrettLvovsky2022, ):

(7) $\displaystyle\nabla_{\vec{\theta}}\tilde{E}=2{\rm
Re}\left(\mathbb{E}_{p}\left[\left(E_{loc}(\mathbb{x})-\mathbb{E}_{p}\left[E_{loc}(\mathbb{x})\right]\right)\nabla_{\vec{\theta}}\ln\left(\Psi_{\vec{\theta}}^{\ast}(\mathbb{x})\right)\right]\right).$

Similarly, $\nabla_{\vec{\theta}}\tilde{E}$ is used to approximate the exact
gradient $\nabla_{\vec{\theta}}E$. After that, the parameters $\vec{\theta}$
can be updated based on $\nabla_{\vec{\theta}}\tilde{E}$ and the optimizer,
finishing one iteration of the VMC algorithm. The workflow of the serial NNQS
algorithm is demonstrated in Fig. 1.

###  2.2. Autoregressive neural networks

A high-dimension probability distribution $\pi(\mathbb{x})$ is said to be
autoregressive if it can be evaluated as follows:

(8)
$\displaystyle\pi(\mathbb{x})=\pi_{x_{1}}\pi_{x_{2}|x_{1}}\cdots\pi_{x_{N}|x_{N-1},\dots,x_{1}},$

where $N$ is the length of $\mathbb{x}$ and $x_{i}$ the $i$-th element of
$\mathbb{x}$, $\pi_{x_{i}|x_{i-1},\dots,x_{1}}$ denotes the conditional
probability distribution of $x_{i}$ on the observed sequence from $x_{1}$ to
$x_{i-1}$. The autoregressive property makes it very easy to sample from
$\pi(\mathbb{x})$: one can first sample from the single-variate probability
distribution $\pi(x)$ and get an outcome $x_{1}$, then one can sample from
another single-variate probability distribution $\pi(x|x_{1})$ and get an
outcome $x_{2}$, following this procedure one can obtain one sample
$\mathbb{x}=\\{x_{1},\dots,x_{N}\\}$ from $\pi(\mathbb{x})$ by only performing
$N$ local samplings from $N$ single-variate probability distributions. This
sampling algorithm will be referred to as the autoregressive sampling
(demonstrated in Fig. 3(a)).

The advantage of using autoregressive neural networks based NNQS is at least
two-fold: 1) It is efficient in the sense that only local samplings are
required (the local state space size is only $2$ for spin systems); 2) The
obtained samples are exact in the sense that they are free from sampling
inaccuracies such as the auto-correlation problem of the MC sampling, and
importantly, we observe in practice that we can often easily converge to the
ground state without using the SR technique, thus one can use very deep neural
networks with a large number of parameters and scale up the algorithm to a
large number of processes.

###  2.3. Quantum Chemistry Hamiltonians

Once a finite single-electron basis set has been chosen, one can expand the
continuous Hamiltonian in this basis and obtain a second quantized Hamiltonian
in the following form:

(9) $\displaystyle
H^{e}=\sum_{p,q}h^{p}_{q}\hat{a}^{\dagger}_{p}\hat{a}_{q}+\frac{1}{2}\sum_{p,q,r,s}g^{p,q}_{r,s}\hat{a}^{\dagger}_{p}\hat{a}^{\dagger}_{q}\hat{a}_{r}\hat{a}_{s},$

with $h^{p}_{q}$ and $g^{p,q}_{r,s}$ the one- and two-electron integrals,
$\hat{a}^{\dagger}_{p}$ and $\hat{a}_{q}$ the creation and annihilation
operators. For convenience of using arbitrary neural networks, the JW
transformation can be used to convert $H^{e}$ into a many-spin Hamiltonian,
which generally takes the form

(10) $\displaystyle H^{s}=\sum_{i=1}^{N_{h}}c_{i}P_{i},$

where each $P_{i}$ is the tensor product of Pauli spin operators
$\\{I,X,Y,Z\\}$ of length $N$, referred to as a Pauli string, and $c_{i}$ is a
real coefficient. $N_{h}$ denotes the total number of Pauli strings. We note
that the integrals in Eq.(9) can be calculated using the package PySCF (pyscf,
) while the JW transformation can be done using the package OpenFermion
(openfermion, ). For quantum chemistry Hamiltonians, $N_{h}$ often scales as
$O(N^{4})$, which means that for each input bitstring $\mathbb{x}$, there
could exist $O(N^{4})$ $\mathbb{x}^{\prime}$s with nonzero
$H^{s}_{\mathbb{x}\mathbb{x}^{\prime}}$. As a result for large $N$, evaluating
the local energy in Eq.(4) can be very expensive and storing all the
$\mathbb{x}^{\prime}$s that are nontrivially coupled to $\mathbb{x}$ could use
a large amount of memory. In this work, we only use the spin Hamiltonian
$H^{s}$, and we will denote it as $H$ for simplicity in the following.

###  2.4. Current state of the art

Electronic structure calculation using NNQS in the second quantized formalism
is still in its early stage. The RBM was first introduced for electronic
structure calculations in 2020, where the largest molecular system studied is
the H2O molecular system in the 6-31G basis ($26$ qubits) due to the high
sampling cost (ChooCarleo2020, ). The autoregressive neural network is first
introduced as the wave function ansatz in NNQS to study quantum many-spin
systems also in 2020 (SharirShashua2020, ), and a parallel implementation of
autoregressive neural network for many-spin systems is demonstrated in 2021
(ZhaoVeerapaneni2021, ). Autoregressive neural network is first introduced for
electronic structure calculations in 2022, which uses a multilevel perceptron
(MLP) with hard-coded pre- and postprocessing to ensure the autoregressive
property and is referred to as neural autoregressive quantum state (NAQS)
(BarrettLvovsky2022, ). It is demonstrated that NAQS can reach superior
accuracy for all the molecular systems studied with up to $30$ qubits compared
to CC methods and RBM. A parallel implementation of an autoregressive neural
network named MADE (masked autoencoder for distribution estimation (MADE, )),
with optimized local energy calculation on GPU, pushed the simulation scale to
$76$ qubits (the CNa2O3 molecular system) (Zhao_2023, ). However in the last
work the accuracy is about the same or only mildly improved against CCSD
results for all the test cases, and is worse than Ref. (BarrettLvovsky2022, )
for molecular systems within $30$ qubits.

To summarize, for electronic structure calculations using NNQS
(BarrettLvovsky2022, ) in the second quantized formalism, NAQS gives the most
accurate results up to date, but the simulation scale is currently limited
within $30$ qubits. MADE (Zhao_2023, ) with GPU optimized local energy
calculation reaches the largest simulation scale ($52$ electrons, $76$
qubits), but is not as accurate as NAQS.

##  3\. Innovations

###  3.1. QiankunNet: Transformer based deep neural network as wave function
ansatz

Figure 2. QiankunNet: The transformer based architecture for the wave function
ansatz, where the transformer is used for the amplitude and the multilevel
perception is used for the phase.

The difficulty of solving quantum many-body problems originates from the fact
that quantum states live in exponentially large Hilbert spaces ($2^{N}$ for
$N$ qubits). To tackle this difficulty we borrow ideas from classical machine
learning, where there also exist many problems for which the solutions live in
exponentially large state spaces. As an example, in NLP a sentence could
contain many words, each of which may be taken from a vocabulary of size
$30522$. If one sets a maximal sentence length to be $512$, then the possible
sentences constitute a state space of size $30522^{512}$ (devlin2018bert, ),
which is larger than the sizes of most quantum many-body problems that can be
numerically solved. Nevertheless, the patterns that different sentences can
appear are usually limited by a few (unknown) grammatical rules, and with
state-of-the-art language modeling techniques using the transformer-based
architectures (e.g. ChatGPT) (radford2018gpt, ; radford2019gpt2, ), these
rules can already be learned with high precision.

In the context of NLP, the language models are often designed such that they
are generative models which are able to accept variable-length inputs and/or
generate variable-length outputs, and they are also designed to be
autoregressive with a natural interpretation: given a sequence of observed
words, the model tries to guess the next word by sampling from the estimated
probability distribution of it conditioned on the observed ones (“given the
text so far, what should the next word be?”). From the measurement point of
view, a quantum state is similar to a language model: one could obtain a
sample (bitstring) from the quantum state by measuring each qubit of it
sequentially, and in general, the $i$-th measurement outcome will be dependent
on all the previous observed outcomes because of quantum correlations.

It is known that the transformer is able to learn distant correlations with a
constant number of operations (Transformer2017, ), while in comparison for
other neural networks such as CNN the number of required operations grows
linearly with the distance (JohnKremer2001, ). Inspired by the enormous
successes made by the transformer based architectures in NLP and other machine
learning tasks, we use a customized transformer based deep neural network as
the wave function ansatz. Concretely, we decompose the amplitude part and
phase part of the wave function ansatz

(11) $\displaystyle\Psi(\mathbb{x})=|\Psi(\mathbb{x})|e^{{\rm
i}\phi(\mathbb{x})},$

and then we use transformer to represent the probability
$|\Psi(\mathbb{x})|^{2}$, and use MLP to represent the phase
$\phi(\mathbb{x})$. In comparison, Refs. (BarrettLvovsky2022, ; Zhao_2023, )
have used MLP based neural networks for both the amplitude and the phase. Our
wave function ansatz is demonstrated in Fig. 2. In comparison with the general
transformer, we have only used decoders since our problem is unsupervised. Our
model can be easily made very deep by stacking a large number of decoders. In
addition, our ansatz is autoregressive since only the amplitude contributes to
the probability $|\Psi(\mathbb{x})|^{2}$ and the amplitude part is
autoregressive by design, therefore we can use very efficient sampling
algorithms designed for autoregressive models.

Figure 3. (a) The autoregressive sampling algorithm which generates one sample
per run. (b) The batch autoregressive sampling algorithm which generates
$N_{s}$ samples per run. We have assumed that two qubits are sampled in each
local sampling step, thus a quadtree. Each circle in (b) corresponds to a
particular local sampling outcome, and the number on the edge pointing to the
circle means the weight. $N_{s}$ can be chosen to be any number ($N_{s}=1000$
is used as an example).

The autoregressive sampling returns one sample per run, during which $N$ local
samplings are performed. For autoregressive wave function ansatz Ref.
(BarrettLvovsky2022, ) proposes a batch autoregressive sampling (BAS)
algorithm which is much more efficient than the autoregressive sampling. The
basic idea of BAS is that in each local sampling step of the autoregressive
sampling, one generates a batch of samples instead of a single one.
Concretely, one first generates $N_{s}$ samples ($N_{s}$ could easily be as
large as $10^{12}$) from the first local sampling step (namely sampling from
$\pi(x)$), and only stores the unique samples with their weights. Then in the
$i$-th step, one computes the probability distribution
$\pi_{x|x_{i-1},\dots,x_{1}}$ for each unique sample $x_{i-1}\cdots x_{1}$
generated in the previous step (assuming that $w_{x_{i-1},\dots,x_{1}}$ is the
corresponding weight), and performs local samplings from
$\pi_{x|x_{i-1},\dots,x_{1}}$ to obtain exactly $w_{x_{i-1},\dots,x_{1}}$
samples. The intermediate samples with zero weights are pruned in each layer.
In this way, one can generate $N_{s}$ samples at the $N$-th step, stored as
$N_{u}$ unique samples with their weights.

The autoregressive sampling and the BAS are demonstrated in Fig. 3. For
transformer based architectures, calculating the conditional probability
distribution $\pi(x|x_{i-1},\dots,x_{1})$ has a computational cost that
roughly scales as $O(i^{2})$. As a result the cost of the autoregressive
sampling for $N_{s}$ samples is $O(N_{s}N^{3}/3)$. For the BAS, assuming that
there are $N_{u,i}$ unique samples generated at the $i$-th step (which is the
width of the $i$-th layer in Fig. 3(b)), then the total computational cost is
$\sum_{n=1}^{N}i^{2}N_{u,i}$. If we assume that $N_{u,i}$ is roughly a
constant, denoted as $N_{u}$ ($N_{u}$ is usually less than $10^{6}$ even for
molecular systems with $N>100$), then the cost of BAS is only
$O(N_{u}N^{3}/3)$, which is independent of $N_{s}$. Existing researches on the
acceleration of training (lightseq2-sc22, ) and inference (lightseq-infer, ;
turbotransformers, ) for Transformer models on GPU platforms can be leveraged
to design the GPU-accelerated implementation for the sampling algorithm.

###  3.2. Data-centric parallelization strategy for NNQS

Figure 4. Schematic of the data-centric parallelization scheme for NNQS using
autoregressive wave function ansatz. Each process in the first level
parallelization takes charge of a number of unique samples with their weights
throughout the computation.

We propose a customized data-centric parallelization scheme for the VMC
algorithm which is valid for all autoregressive neural networks, as
demonstrated in Fig. 4 and in comparison with the serial scheme in Fig. 1. Our
parallelization scheme contains two levels. In the first level
parallelization, we divide all the unique samples ($N_{u}$) into $N_{p}$
batches and parallelize the sampling, local energy and gradient calculation
(backward propagation) over those batches, such that each functions processes
$N_{u}/N_{p}$ unique samples (instead of $N_{s}/N_{p}$ samples). In the second
level parallelization, the sampling and backward propagation functions for
each batch are further parallelized using pytorch multi-threading, while for
local energy we use our parallel local energy evaluation scheme (see Sec. 3.4)
over $N_{t}$ threads by further divide the $N_{u}/N_{p}$ samples into $N_{t}$
smaller batches.

More concretely, our parallelization scheme contains $6$ stages in each VMC
iteration: 1) We use a customized parallel BSA strategy (see Sec.3.3) over a
total of $N_{p}$ processes to obtain approximately $N_{u}/N_{p}$ unique
samples in each process; 2) We use MPI_Allgather to synchronize all the
samples on each process, with $N_{u}N_{p}(\lceil N/8\rceil+16)$ Byte data
communication; 3) Each process takes a chunk of $N_{u}/N_{p}$ unique samples
and performs the parallel local energy evaluation scheme; 4) Collecting all
the local energies together into the main process, computing the average
therein and then broadcasting it to all the processes using MPI_Allreduce,
with $16N_{p}$ Byte data communication; 5) Performing backward propagation on
each process independently; 6) Collecting all the local gradients from each
process into the main process and calculating the weighted average, then
updating the parameters using the gradient-based optimizer and send the new
parameters to all the processes using MPI_Allreduce, with $8MN_{p}$ Byte for
$M$ parameters data communication. Importantly, during the whole process of
our parallel VMC iteration, only steps 2), 4), 6) involve very minimal data
communication (Typically C2 with STO-3G basis set, $N=20$, $N_{u}=2.7\times
10^{4}$, $N_{p}=64$ and $M=2.7\times 10^{5}$, the total data communication
volume is about $173$ MB in our one iteration calculations), and the unique
samples are always kept in their processes until all calculations finish
(therefore data-centric).

###  3.3. Parallel batch autoregressive sampling algorithm

Figure 5. The parallel batch autoregressive sampling algorithm. The tasks are
divided at the $k$-th local sampling step ($k$-th layer of the tree) such that
the number of samples on each process are approximately equal. Pruning is done
in each step to remove leaves with zero occurrence and which break the number
conservation constrain.

The first challenge in NNQS is the sampling process. Autoregressive neural
networks allow one to generate exact samples without the need of a pre-
thermalization step and throwing away intermediate samples (to avoid auto-
correlations) necessary for MC sampling, thus are more efficient than MC
sampling in general. Moreover, the batch sampling algorithm can output a large
number of unique samples together with their occurrences in a single run (see
Fig. 3). However paralleling the BSA on distributed architectures is not
straightforward, mainly because 1) in BSA the calculation in one local
sampling step is dependent on all the previous steps and 2) the positions and
the total number of the unique samples in each step are not known before hand.

Figure 6. (a) Hamiltonian for the H2 molecular system with $4$ qubits and $2$
electrons stored as a list of Pauli strings. (b) The scheme used in Ref.
(Zhao_2023, ) to store the Hamiltonian in memory. (c) Our compressed data
structure to store the complete information of the Hamiltonian. The data
inside the dashed red boxes are the information stored in memory in both (b)
and (c).

To solve this issue, we propose a heuristic parallelization strategy for BSA
that proves to be effective in practice, which is shown in Fig. 5. We
represent the whole BSA process as a tree where each layer corresponds to a
local sampling step. We sample two qubits each time (thus a quadtree) since
they correspond to the same spatial orbital. We also sample in the reverse
order of the qubits after the JW transformation as in Ref. (Zhao_2023, ). We
first perform the serial BSA on each process using the same random seed (such
that we get exactly the same samples on each process) for the first $k$ steps,
and then we divide the unique samples (nodes in Fig. 5) at the $k$-th step
into $N_{p}$ processes such that each process gets approximately the same
number of samples (instead of unique samples) in this step. This does not
guarantee that each process generates the same number of unique samples in the
end (thus achieving good load balance since the later calculations are
directly determined by the number of unique samples), but from our numerical
experiments we observe that approximate load balance can be achieved with this
strategy if $k$ is chosen properly. In our implementation we use a simple
strategy to dynamically determine $k$: we set a threshold $N_{u}^{\ast}$ and
choose $k$ to be the first local sampling step such that the current number of
unique samples ( $N_{u,k}$) is larger than $N_{u}^{\ast}$. We note that in
this strategy the calculations from the $1$ to $k-1$-th steps are simply
repeated on each process, however, this will not significantly affect the
overall computational cost since the computational cost of the local sampling
grows as $i^{2}$ against step $i$.

During each step, we prune the leaves with $0$ occurrence. In addition, there
is a number conserving constrain which can be used to reduce the sample space
size: the total number of spin up and spin down electrons, denoted as
$n_{\uparrow}$ and $n_{\downarrow}$, are conserved separately. During the JW
transformation, the $i$-th spatial orbital is mapped to the two qubits at
positions $2i-1$ and $2i$, corresponding to the two spin orbitals with up and
down spins. Denoting each $x_{i}=y_{2i}y_{2i-1}$ with $y_{j}\in\\{0,1\\}$
(since we sample two qubits in each step), this number conservation constrain
can be imposed by regularizing the local conditional probability distributions
$\pi(x_{i}|x_{i-1},\dots,x_{1})$ as (Zhao_2023, ):

(12)
$\displaystyle\tilde{\pi}(x_{i}|x_{i-1},\dots,x_{1})=\begin{cases}0,\text{if
}\sum_{j=1}^{i}y_{2j-1}>n_{\uparrow}\\\ 0,\text{if
}\sum_{j=1}^{i}y_{2j}>n_{\downarrow}\\\
\pi(x_{i}|x_{i-1},\dots,x_{1}),\text{otherwise}\end{cases},$

and then using
$\tilde{\pi}(x_{i}|x_{i-1},\dots,x_{1})/\sum_{x_{i}}\tilde{\pi}(x_{i}|x_{i-1},\dots,x_{1})$
instead.

###  3.4. Accelerating local energy calculation on GPU

The second challenge in NNQS for electronic structure calculations is to
evaluate the local energy in Eq.(4). For molecular systems this step could be
very expensive due to the $O(N^{4})$ scaling of the total number of Pauli
strings. This could also lead to a memory issue since there are $O(N^{4})$
possible $\mathbb{x}^{\prime}$ coupled to an input $\mathbb{x}$. We have used
several techniques to parallelize the computation on a GPU while at the same
time limit the memory usage and reduce the computational cost, such that we
can scale up our calculations to more than $100$ qubits. These techniques are
shown as follows.

(1) A highly compressed data structure to store the Hamiltonian in memory. For
large molecular systems storing the Hamiltonian as a list of symbolic Pauli
strings (or integers, see Fig. 6(a)) would take a large amount of memory and
be inefficient for later calculations. In Ref. (Zhao_2023, ) a compressed data
structure is used to represent the Hamiltonian, where each Pauli string is
represented by three components: two boolean tuples of length $N$ where the
first one stores the occurrence of $X$ or $Y$ and is used to calculate the
states $\mathbb{x}^{\prime}$ coupled to each $\mathbb{x}$, while the second
one stores the occurrence of $Y$ or $Z$ and is used to calculate the
coefficient, plus an integer to store the occurrence of $Y$ used to calculate
the sign. This data structure is shown in Fig. 6(b). We propose a more memory-
efficient scheme as shown in Fig. 6(c), where we only store the unique ones of
the first list of boolean tuples and reorganize the second list accordingly.
The coefficient is precomputed using the occurrence of $Y$ and updated in-
place. The computational cost of the reorganization step can be neglected
since it only needs to be done once and can be pre-computed. The memory
reduction of our scheme in Fig. 6(c) compared to the scheme in Fig. 6(b) is
around $40\%$ in general (see Fig. 9 in Sec. 4.3). In addition, for each input
$\mathbb{x}$, all the Pauli strings with the same first boolean tuple will
couple $\mathbb{x}$ to the same bitstring $\mathbb{x}^{\prime}$, therefore our
data structure naturally allows one to calculate each unique
$\mathbb{x}^{\prime}$ only once (the coefficients for duplicate
$\mathbb{x}^{\prime}$ are simply summed over as shown in Fig. 7(b)). Algorithm
1 meticulously delineates the aforementioned optimizations, the initial step
entails the parsing of inputted Pauli strings, which serves to construct Pauli
matrices XY and YZ. Simultaneously, occurrences of Y are enumerated to compute
the new coefficients. Following this, a dictionary data structure (Dict) and
state encoding are employed to eliminate redundant states. Lastly, continuous
memory space is allocated. The dictionary-type data is compressed into this
continuous space for storage, and corresponding ranges are retained using an
index array (idxs).

(2) A fused nonzero Hamiltonian entry evaluation and local energy calculation
design. That is, when evaluating $E_{loc}(\mathbb{x})$ using Eq.(4), instead
of finding all the different $\mathbb{x}^{\prime}$ with nonzero
$H_{\mathbb{x}\mathbb{x}^{\prime}}$ (together with their coefficients) and
then using them to evaluate $E_{loc}(\mathbb{x})$, we directly calculate the
contribution of each unique $\mathbb{x}^{\prime}$ to $E_{loc}(\mathbb{x})$
(namely
$H_{\mathbb{x}\mathbb{x}^{\prime}}\Psi_{\vec{\theta}}(\mathbb{x}^{\prime})/\Psi_{\vec{\theta}}(\mathbb{x})$
which is a scalar) as soon as it is found, and then sum over them. In this
approach we can calculate $E_{loc}(\mathbb{x})$ by only storing two bitstrings
in memory.

Input: pauliStrings, n_qubits

Result: idxs, pmXY_buf, pmYZ_buf, coeffs_buf

1

2$N,K\leftarrow\text{n\\_qubits},\text{length(pauliStrings)}$

3 $\text{pmXY, pmYZ}\leftarrow\text{zeros}(N,K),\text{zeros}(N,K)$

4

5for  _( $i$, (coeff, pauliStr)) in enumerate(pauliStrings)_ do

$Y_{occ}\leftarrow 0$ // occurrence for Y

/* Construct Pauli mat XY and Pauli mat YZ */

6 for  _( $p$, pauli) in pauliStr_ do

7 if  _pauli $=="X"$_ then

8 pmXY[$p$, $i$] $\leftarrow 1$

9

10 else if  _pauli $=="Y"$_ then

11 pmXY[$p$, $i$] $\leftarrow$ pmYZ[$p$, $i$] $\leftarrow 1$

12 $Y_{occ}\leftarrow Y_{occ}+1$

13

14 else

15 pmYZ[$p$, $i$] $\leftarrow 1$

16

/* New coefficients: fuse const calculation */

17 coeff $\leftarrow\text{real}(coeff)\times\text{real}((-i)^{Y_{occ}})$

18

/* Compressed Pauli mat XY using Dict */

19 sid = _encode2id(pmXY[:,$i$])

20 if  _not haskey(coeffs_dict, sid)_ then

21 pmXY_dict[sid] $\leftarrow$ pmXY[:,$i$]

22

23 coeffs_dict[sid].append(coeff)

24 pmYZ_dict[sid].append([pmYZ[:,$i$]])

25

26

/* Compacting data into new continuous buffer */

27

28for  _( $i$, (sid, pm23)) in enumerate(pmYZ_dict)_ do

29 $l\leftarrow$ length(pm23)

30 pmYZ_buf[:,$j:j+l$] $\leftarrow$ pm23

31 pmXY_buf[:, $i$] $\leftarrow$ pmXY_dict[sid]

32 coeffs_buf[$j$:$j$+$l$] $\leftarrow$ coeffs_dict[sid]

33 $j\leftarrow j+l$; idxs[$i$] $\leftarrow j$

34

Algorithm 1 Preprocess: Hamiltonian data structure optimizations

Figure 7. (a) Parallelizing the calculation of local energies over the
samples. (b) Calculating the local energy for a given input sample.

(3) Parallelizing the calculation of local energies over unique samples on a
GPU. Each GPU is responsible for all the calculations related to a batch of
$l$ unique samples ($l=N_{u}/N_{p}$). GPU parallelization is implemented by
dividing the $l$ unique samples into $N_{t}$ smaller batches (each containing
$l/N_{t}$ unique samples) and calculating the local energies of each batch on
one thread, which is shown in Fig. 7.

(4) A sample-aware scheme for evaluating $E_{loc}(\mathbb{x})$. During the
local energy calculation, all the unique samples, denoted as
$\mathcal{S}=\\{\mathbb{x}^{1},\mathbb{x}^{2},\cdots,\mathbb{x}^{N_{u}}\\}$,
are stored on each GPU (see Sec. 3.2). When evaluating $E_{loc}(\mathbb{x})$,
we calculate the coefficient $H_{\mathbb{x}\mathbb{x}^{\prime}}$ only if
$\mathbb{x}^{\prime}\in\mathcal{S}$, as illustrated in Fig 7(b) (see Fig. 10
in Sec. 4.3 for actual speedup using this technique).

(5) Efficient storage of the unique samples $\mathcal{S}$ as a lookup table.
The $l$ samples generated on each GPU are stored as boolean lists for
computational efficiency as can be seen from Fig. 7(b). However, when
synchronizing the unique samples on all the GPUs, a memory efficient storage
scheme is preferable to reduce data communication, therefore we simply encode
the boolean tuple into the bits of a $64$-bit integer (in case $64\leq N<128$,
we use two integers). This storage format also allows us to use the very
efficient binary search during the sample-aware local energy evaluation, since
the integers can be naturally ordered according to their values.

Algorithm 2 includes optimization techniques (2)-(5), each of which assigns a
subset of samples to a single GPU thread. The computation process prioritizes
coupled states. Within the identification lookup table (id_lut), a binary
search is conducted. If a matching state exists, the computation proceeds to
calculate the coupling coefficients. By accessing the wavefunction lookup
table (wf_lut), the necessity for forward inference in the neural network is
bypassed, allowing for immediate contributions to be aggregated. In practical
implementation, complex numbers are separated into their real and imaginary
components for independent computation, and the resultant values are
subsequently written into their respective positions.

Input: n_qubits, $NK$, idxs, coeffs, pmXY, pmYZ, batch_size_cur_rank, $ist$,
states, id_lut, wf_lut

Result: res_elocs

/* Initialization of power of two table */

1 tbl_pow2$[0...63]\leftarrow\\{1,2,...,2^{63}\\}$

/* Partition for CUDA kernel */

2 index $\leftarrow$ blockIdx.x $\times$ blockDim.x $+$ threadIdx.x

3 stride $\leftarrow$ gridDim.x $\times$ blockDim.x

4

5for  _$ii\leftarrow$ index to batch_size_cur_rank step stride_ do

6 eloc_real, eloc_imag $\leftarrow 0$

7 for  _$k\leftarrow 0$ to $NK$_ do

/* Fuse coupled states and state2id */

8 id $\leftarrow 0$

9 for  _$j\leftarrow 0$ to $N$_ do

10 id $+=$ xor(states[$ii\times N+j$], pmXY[$k\times N+j$]) $\times$
tbl_pow2[$j]$

11

12 end for

/* Binary find id among the lookup table */

13 idx $\leftarrow$ binary_find(id_lut, id)

14 if  _idx does not exist_ then

15 continue

16

17 end if

/* Calculate coupled coef */

18 coef $\leftarrow 0.0$

19 for  _$i\leftarrow\rm{idxs}[k]$ to $\rm{idxs}[k+1]$_ do

20 $\rm{sum}\leftarrow 0$

21 for  _$j\leftarrow 0$ to $N$_ do

22 sum $+=$ states[$ii\times N+j$] & pmYZ[$i\times N+j]$

23

24 end for

25 coef $+=$ ( (sum & $1)?1:-1)\times$ coeffs[$i]$

26

27 end for

/* Accumulation */

28 eloc_real $+=$ coef $\times$ wf_lut[idx$\times 2]$

29 eloc_imag $+=$ coef $\times$ wf_lut[idx$\times 2+1]$

30

31 end for

/* Store the result number as return */

32 a, b $\leftarrow$ eloc_real, eloc_imag

33 c, d $\leftarrow$ wf_lut[$(ist+ii)\times 2$], wf_lut[$(ist+ii)\times 2+1$]

34 c2_d2 $\leftarrow$ c2 \+ d2

35 res_elocs[$ii\times 2]\leftarrow$ (a $\times$ c + b $\times$ d) $/$ c2_d2

36 res_elocs[$ii\times 2+1]\leftarrow-$(a $\times$ d $-$ b $\times$ c) $/$
c2_d2

37

38 end for

Algorithm 2 Optimizations of local energy calculation

##  4\. Evaluation

###  4.1. Experimental setup

The wave function ansatz used in all the simulations in this work is set as
follows. For the amplitude part, we have used two decoders with $d_{model}=16$
(embedding size) and $n_{heads}=4$ (number of heads). For the phase part, we
have used three dense layers in the MLP with sizes $N\times 512\times
512\times 1$. We achieved optimal outcomes exploring MLP layers (3-5), hidden
layers (256-1024), and decoder layers (2-4) with a $d_{model}$ range of 16-64.
We have used the gradient descent optimizer AdamW for training with the learn
rate schedule

(13) $\displaystyle\alpha_{i}=d_{model}^{-0.5}\times\min(i^{-0.5},i\times
S_{warmup}^{-1.5}),$

where $\alpha_{i}$ is the learn rate of the $i$-th training epoch and we set
the warm up steps as $S_{warmup}=4000$. The maximum number of VMC iterations
is set to $10^{5}$. All our simulations are performed on AMD EPYC 7742 CPU and
NVIDIA A100 PCIe 80GB for CPU and GPU computation environment respectively.

During the training process, we start from a wave function ansatz with
randomly initialized parameters. In the pre-training stage, we observe that a
large number of unique samples will be generated since the wavefunction is far
from exact, thus we set a lower threshold $N_{s}$ ($N_{s}=10^{5}$ in the first
$100$ VMC iterations) for the number of samples to increase efficiency. In the
later stage we set the threshold (gradually increase $N_{s}$ until it reaches
the maximum value of $10^{12}$) to be very large for accurate calculation (we
observe that in the later stage, only a smaller number of unique samples will
be generated). The efficacy of the model is assessed based on convergence
precision.

###  4.2. Simulation Validation

Molecule | $N$ | $N_{e}$ | $N_{h}$ | HF | CCSD | NAQS (BarrettLvovsky2022, ) | MADE (Zhao_2023, ) | QiankunNet | FCI
---|---|---|---|---|---|---|---|---|---
H2O | 14 | 10 | 1390 | $-74.964$ | -75.0154 | -75.0155 | -75.0155 | -75.0155 | -75.0155
N2 | 20 | 14 | 2239 | $-107.4990$ | -107.6560 | -107.6595 | -107.6567 | -107.6601 | -107.6602
O2 | 20 | 16 | 2879 | $-147.6319$ | -147.7027 | -147.7500 | -147.7499 | -147.7501 | -147.7502
H2S | 22 | 18 | 9558 | $-394.3114$ | -394.3545 | -394.3546 | -394.3545 | -394.3546 | -394.3546
PH3 | 24 | 18 | 24369 | $-338.6341$ | -338.6981 | -338.6984 | -338.6981 | -338.6982 | -338.6984
LiCl | 28 | 20 | 24255 | $-460.8273$ | -460.8475 | -460.8496 | -460.8481 | -460.8494 | -460.8496
Li2O | 30 | 14 | 20558 | $-87.7956$ | -87.8855 | -87.8909 | -87.8856 | -87.8907 | -87.8927
MAE (Hartree) |  |  |  |  |  | $3.8\times 10^{-4}$ | $1.8\times 10^{-3}$ | $3.7\times 10^{-4}$ |
Table 1. Ground state energies (in Hartree) calculated by our method. The
conventional HF and CCSD results, FCI results as well as existing NNQS results
including NAQS (BarrettLvovsky2022, ) and MADE (Zhao_2023, ) are shown for
comparison. $N$ is the number of qubits, $N_{e}$ the total number of electrons
(including spin up and spin down) and $N_{h}$ the total number of Pauli
strings of the Hamiltonian. The mean absolute errors (MAE) for various methods
compared with FCI are also listed.

To demonstrate the precision of our transformer-based NNQS method, we first
compute the ground state energies of several small-scale molecular systems and
compare them to existing results using NNQS, which are shown in Table. 1. The
mean absolute error (MAE) are listed for each method. We can see that our
method can reach about the same precision as NAQS (for N2 we are more precise)
and is generally more accurate than MADE.

Figure 8. (a) Potential energy surface of the BeH2 molecular system in the
STO-3G basis set calculated by our method (green circles), compared to the HF,
CCSD and FCI energies. (b) The absolute errors with respect to FCI results.

We also apply our method to calculate the potential energy surface of the BeH2
molecular system in the STO-3G basis set ($14$ qubits), which is shown in Fig.
8. We can see that chemical accuracy (within an absolute error of $1.6\times
10^{-3}$ Hartree compared to FCI) can be reached in the potential energy
surface simulation.

###  4.3. Memory saving and speedup of local energy calculation

Figure 9. Memory reduction using our compressed data structure to store the
Hamiltonian compared to the original scheme from Ref. (Zhao_2023, ) for
several exemplary molecular systems in the STO-3G basis set.

Now we show the memory reduction of using our compressed data structure for
the Hamiltonian (see Fig. 6(c) and Sec. 3.4), compared to the scheme proposed
in Ref. (Zhao_2023, ) (see Fig. 6(b)), which is plotted in Fig. 9 for a wide
variety of molecular systems. We can see that generally there is a more than
$40\%$ memory reduction.

Figure 10. The runtimes for local energy calculation by integrating different
optimization techniques, for the C2, LiCl and C2H4O molecular systems in the
STO-3G basis set with $N=20,28,38$. The total number of unique samples
produced in these simulations are $N_{u}=10553,9271,25622$ respectively. SA
means the sample-aware evaluation scheme (Sec. 3.4, method (4)), FUSE means
the fused design of nonzero Hamiltonian entry evaluation and local energy
calculation (Sec. 3.4, method (2)), LUT means to store the unique samples as
an efficient lookup table (Sec. 3.4, method (5)), GPU means the local energy
calculation implemented on GPU (Sec. 3.4, method (2)). The baseline is a bare
CPU version without using the optimization techniques mentioned here. The
other techniques in Sec. 3.4 which are not mentioned here are used in all
tests.

In Fig. 10, we show the speedups of the local energy calculation induced by
different optimization techniques, for the C2, LiCl, C2H4O molecular systems.
Taking C2 as an example, first, the sample aware evaluation scheme (reducing
the computational cost) plus the fused design of nonzero Hamiltonian entry
evaluation and local energy calculation (increasing data locality and reducing
memory usage), yields a performance gain of $24$x. Second, by storing the
unique samples as an ordered integer list, which reduces memory usage and
increases search efficiency, the performance gain further increases to $103$x.
Third, by using GPU parallelization, the overall performance gain becomes
$3768$x compared to the baseline CPU version where none of these optimizations
are used. Overall, implementing our parallel local energy calculation scheme
on GPU gives us around $4,000$x speed up compared to the baseline CPU version,
which is essential to scale up the simulation scale.

###  4.4. Scalability Results

Figure 11. Strong scaling of the computation for the benzene molecular system
in the 6-31G basis set ($120$ qubits). The green line shows the strong
parallelization efficiency. Three functions (sampling, local energy
calculation, backpropagation) are profiled in the computation time.

We use the benzene molecular system in the 6-31G basis set with $120$ qubits
for the scalability test, for which $N_{h}=2434919$. We measure the
performance after the pre-training stage such that the number of unique
samples in each iteration becomes more or less stable. We have used
$N_{u}^{\ast}=16384n$ for $n$ GPUs in the parallel BAS for both the strong and
weak scaling test.

Fig. 11 shows the strong scaling performance of our method, where we have used
$N_{s}=1.6\times 10^{6}$ and obtain $N_{u}\approx 6.5\times 10^{5}$. As the
number of used GPUs increases from $4$ to $32$, the strong parallelization
efficiency is still higher than $84\%$. While when increased to $64$ GPUs, the
efficiency decreases to $67.7\%$. This slightly lower parallel efficiencies
can be attributed to two factors: First the relatively small number of unique
samples ($\approx 10^{4}$) that is allocated on each GPU when distributing
$6.5\times 10^{5}$ unique samples over $64$ GPUs, the parallelization
efficiency of GPU is not well utilized at this scale; Second, during the
parallel BAS process, the imbalance of the number of unique samples generated
on each GPU could happen due to the pruning of the samples with physical
constrains. We note that the later the parallel BAS is performed, the more
evenly distributed the unique samples will be, yet the time for prior serial
sampling process will increase.

Figure 12. Weak scaling of the runtime for the benzene molecular system in the
6-31G basis set. To ensure that a approximately equal number of unique samples
($\approx 2.04\times 10^{4}$) is generated on each GPU, we have set
$N_{s}=5n\times 10^{4}$ for $n$ GPUs. The green line shows the weak
parallelization efficiency.

Fig. 12 shows the weak scaling performance of our method. We time the
computation time with an increasing number of GPUs from $4$ to $64$, with each
GPU calculating an approximately equal number of unique samples ($2.04\times
10^{4}$). Our results show good weak scaling performance. Taking $32$ GPUs as
a reference, we achieve a weak parallelization efficiency of approximately
$93\%$ compared to $4$ GPUs.

The batched sampling algorithm we adopted takes advantage of the
autoregressive properties to greatly reduce sampling overhead. But it could
sacrifice parallelization due to the tree structure of our sampling process.
The non-exact equivalence between the number of unique samples and the
computational cost in our algorithm will generally affect both the strong and
weak scaling at a larger scale when the fluctuations of the number of samples
on GPUs become larger. For even larger scale parallelization in the future
implementation, one could still take advantage of the conventional Monte Carlo
sampling by simply implementing several independent the batch sampling
algorithm, which will be effective as long as a larger number of unique
samples are going to be important for that problem (the latter approach will
have computational overhead if the number of unique samples is not large
enough).

##  5\. Implications

Figure 13. Potential energy surfaces of the H2 molecular system in the cc-pVTZ
(a) and aug-cc-pVTZ (c) basis sets calculated by our method, compared to the
HF, CCSD and FCI energies. The energies calculated in the complete basis set
limit (the blue dashed lines) are also plotted for comparison. The
corresponding absolute errors with respect to FCI results are shown in (b,d)
respectively.

Finally, we apply our method to calculate the potential energy surfaces of the
hydrogen molecular system in very large basis set: the cc-pVTZ ($56$ qubits)
and aug-cc-pVTZ ($92$ qubits) basis sets. The results are shown in Fig. 13.
This system is chosen because it can demonstrate the scalability of our method
while being exactly solvable by FCI (in our method we have not used the
knowledge that only a small number of determinants are necessary for this
system, thus the computational cost for each VMC iteration is almost as
difficult as for other large systems). For both the two large basis sets our
results are in excellent agreement with the FCI results, and chemical
accuracies are reached in all these simulations. We also show results obtained
with FCI in the complete basis set limit, which can be considered as the
ground truth for the potential energy curve of the hydrogen molecule. We can
see that in these larger basis sets the potential energy curve already well
agree the exact dissociation limit (the mean errors of our calculations with
respect to the ground truth are $2.27\times 10^{-3}$ Hartree for $56$ qubits
and $2.34\times 10^{-3}$ Hartree for $92$ qubits).

##  6\. Conclusion

In summary, we have demonstrated a high-performance NNQS method using the
transformer-based wave function ansatz for ab initio electronic structure
calculations. To our knowledge, this is the first work which introduces
transformer based architectures into NNQS to solve quantum chemistry problems.
To scale up the NNQS method on distributed computing architectures, we propose
a parallel batch autoregressive sampling strategy which works on multiple GPUs
and a highly efficient local energy evaluation scheme which is parallelized on
each GPU using multi-threading. The overall parallelization is designed to be
data centric, namely each GPU manages a batch of unique samples and is
responsible for all the calculations related to this batch throughout the
computation, therefore minimizing data communications among different
processes. The strong and weak scaling of our method up to $64$ GPUs
demonstrate the good scalability of our method. Our work paves the way of
using NNQS to study strongly correlated quantum chemistry systems with more
than $100$ qubits.

Our method is mostly based on the pytorch framework and is applicable for all
autoregressive wave function ansatz. It can also be easily portable to other
high-performance computing architectures as long as the local energy
calculation function is implemented on each process.

###### Acknowledgements.

This work was supported by the National Natural Science Foundation of China
(Grant No. T2222026, 22003073, 11805279). The corresponding authors of this
paper are Honghui Shang
<EMAIL_ADDRESS>and Chu Gu (guochu604b@gmail.com).

## References

  * [1] Konstantinos D. Vogiatzis, Dongxia Ma, Jeppe Olsen, Laura Gagliardi, and Wibe A. de Jong.  Pushing configuration-interaction to the limit: Towards massively parallel mcscf calculations.  The Journal of Chemical Physics, 147(18):184111, 2017.
  * [2] C. David Sherrill and Henry F. Schaefer.  The configuration interaction method: Advances in highly correlated approaches.  volume 34 of Advances in Quantum Chemistry, pages 143–269. Academic Press, 1999.
  * [3] Ireneusz W Bulik, Thomas M Henderson, and Gustavo E Scuseria.  Can single-reference coupled cluster theory describe static correlation?  Journal of chemical theory and computation, 11(7):3171–3179, 2015\.
  * [4] Steven R. White.  Density matrix formulation for quantum renormalization groups.  Physical Review Letters, 69:2863–2866, Nov 1992.
  * [5] Steven R. White.  Density-matrix algorithms for quantum renormalization groups.  Physical Review B, 48:10345–10356, Oct 1993.
  * [6] Jiri Brabec, Jan Brandejs, Karol Kowalski, Sotiris Xantheas, Örs Legeza, and Libor Veis.  Massively parallel quantum chemical density matrix renormalization group method.  Journal of Computational Chemistry, 42(8):534–544, 2021.
  * [7] Henrik R Larsson, Huanchen Zhai, Cyrus J Umrigar, and Garnet Kin-Lic Chan.  The chromium dimer: closing a chapter of quantum chemistry.  Journal of the American Chemical Society, 144(35):15932–15937, 2022\.
  * [8] D. Perez-Garcia, F. Verstraete, M. M. Wolf, and J. I. Cirac.  Matrix product state representations.  Quantum Information Computation, 7(5):401–430, jul 2007.
  * [9] W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal.  Quantum monte carlo simulations of solids.  Reviews of Modern Physics, 73:33–83, Jan 2001.
  * [10] R J Needs, M D Towler, N D Drummond, and P López Ríos.  Continuum variational and diffusion quantum monte carlo calculations.  Journal of Physics: Condensed Matter, 22(2):023201, dec 2009.
  * [11] Brian M. Austin, Dmitry Yu. Zubarev, and William A. Jr. Lester.  Quantum monte carlo and related approaches.  Chemical reviews, 112(1):263–288, 2012.
  * [12] David Pfau, James S. Spencer, Alexander G. D. G. Matthews, and W. M. C. Foulkes.  Ab initio solution of the many-electron schrödinger equation with deep neural networks.  Physical review research, 2:033429, Sep 2020.
  * [13] Jan Hermann, Zeno Schätzle, and Frank Noé.  Deep-neural-network solution of the electronic schrödinger equation.  Nature Chemistry, 12(10):891–897, 2020.
  * [14] Giuseppe Carleo and Matthias Troyer.  Solving the quantum many-body problem with artificial neural networks.  Science, 355(6325):602–606, 2017.
  * [15] Kenny Choo, Titus Neupert, and Giuseppe Carleo.  Two-dimensional frustrated ${J}_{1}\text{$-$}{J}_{2}$ model studied with neural network quantum states.  Physical Review B, 100:125124, Sep 2019.
  * [16] Or Sharir, Yoav Levine, Noam Wies, Giuseppe Carleo, and Amnon Shashua.  Deep autoregressive models for the efficient variational simulation of many-body quantum systems.  Physical Review Letters, 124:020503, Jan 2020.
  * [17] Markus Schmitt and Markus Heyl.  Quantum many-body dynamics in two dimensions with artificial neural networks.  Physical Review Letters, 125:100503, Sep 2020.
  * [18] Dong Yuan, He-Ran Wang, Zhong Wang, and Dong-Ling Deng.  Solving the liouvillian gap with artificial neural networks.  Physical Review Letters, 126:160401, Apr 2021.
  * [19] Xuncheng Zhao, Mingfan Li, Qian Xiao, Junshi Chen, Fei Wang, Li Shen, Meijia Zhao, Wenhao Wu, Hong An, Lixin He, and Xiao Liang.  Ai for quantum mechanics: High performance quantum many-body simulations via deep learning.  In SC22: International Conference for High Performance Computing, Networking, Storage and Analysis, pages 1–15, 2022.
  * [20] Javier Robledo Moreno, Giuseppe Carleo, Antoine Georges, and James Stokes.  Fermionic wave functions from neural-network constrained hidden states.  Proceedings of the National Academy of Sciences, 119(32):e2122059119, 2022.
  * [21] Dong-Ling Deng, Xiaopeng Li, and S. Das Sarma.  Quantum entanglement in neural network states.  Physical Review X, 7:021021, May 2017.
  * [22] Ivan Glasser, Nicola Pancotti, Moritz August, Ivan D. Rodriguez, and J. Ignacio Cirac.  Neural-network quantum states, string-bond states, and chiral topological states.  Physical Review X, 8:011006, Jan 2018.
  * [23] Or Sharir, Amnon Shashua, and Giuseppe Carleo.  Neural tensor contractions and the expressive power of deep neural quantum states.  Physical Review B, 106:205136, Nov 2022.
  * [24] Xun Gao and Lu-Ming Duan.  Efficient representation of quantum many-body states with deep neural networks.  Nature communications, 8(1):662, 2017.
  * [25] Kenny Choo, Antonio Mezzacapo, and Giuseppe Carleo.  Fermionic neural-network states for ab-initio electronic structure.  Nature communications, 11(1):2368, 2020.
  * [26] Thomas D Barrett, Aleksei Malyshev, and AI Lvovsky.  Autoregressive neural-network wavefunctions for ab initio quantum chemistry.  Nature Machine Intelligence, 4(4):351–358, 2022.
  * [27] Tianchen Zhao, James Stokes, and Shravan Veerapaneni.  Scalable neural quantum states architecture for quantum chemistry.  Machine Learning: Science and Technology, 4(2):025034, jun 2023\.
  * [28] Sandro Sorella and Luca Capriotti.  Green function monte carlo with stochastic reconfiguration: An effective remedy for the sign problem.  Physical Review B, 61:2599–2612, Jan 2000.
  * [29] Sandro Sorella, Michele Casula, and Dario Rocca.  Weak binding between two aromatic rings: Feeling the van der waals attraction by quantum monte carlo methods.  The Journal of Chemical Physics, 127(1):014105, 2007.
  * [30] Honghui Shang, Chu Guo, Yangjun Wu, Zhenyu Li, and Jinlong Yang.  Solving schrodinger equation with a language model.  arXiv: 2307.09343 [quant-ph].
  * [31] Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Ł ukasz Kaiser, and Illia Polosukhin.  Attention is all you need.  In I. Guyon, U. Von Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc., 2017.
  * [32] Alec Radford, Karthik Narasimhan, Tim Salimans, Ilya Sutskever, et al.  Improving language understanding by generative pre-training.  2018\.
  * [33] Alec Radford, Jeffrey Wu, Rewon Child, David Luan, Dario Amodei, Ilya Sutskever, et al.  Language models are unsupervised multitask learners.  OpenAI blog, 1(8):9, 2019.
  * [34] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Köpf, Edward Yang, Zach DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala.  PyTorch: An Imperative Style, High-Performance Deep Learning Library.  Curran Associates Inc., Red Hook, NY, USA, 2019.
  * [35] Qiming Sun, Timothy C. Berkelbach, Nick S. Blunt, George H. Booth, Sheng Guo, Zhendong Li, Junzi Liu, James D. McClain, Elvira R. Sayfutyarova, Sandeep Sharma, Sebastian Wouters, and Garnet Kin-Lic Chan.  Pyscf: the python-based simulations of chemistry framework.  WIREs Computational Molecular Science, 8(1):e1340, 2018.
  * [36] Jarrod R McClean, Nicholas C Rubin, Kevin J Sung, Ian D Kivlichan, Xavier Bonet-Monroig, Yudong Cao, Chengyu Dai, E Schuyler Fried, Craig Gidney, Brendan Gimby, Pranav Gokhale, Thomas Häner, Tarini Hardikar, Vojtěch Havlíček, Oscar Higgott, Cupjin Huang, Josh Izaac, Zhang Jiang, Xinle Liu, Sam McArdle, Matthew Neeley, Thomas O’Brien, Bryan O’Gorman, Isil Ozfidan, Maxwell D Radin, Jhonathan Romero, Nicolas P D Sawaya, Bruno Senjean, Kanav Setia, Sukin Sim, Damian S Steiger, Mark Steudtner, Qiming Sun, Wei Sun, Daochen Wang, Fang Zhang, and Ryan Babbush.  Openfermion: the electronic structure package for quantum computers.  Quantum Science and Technology, 5(3):034014, jun 2020.
  * [37] Tianchen Zhao, Saibal De, Brian Chen, James Stokes, and Shravan Veerapaneni.  Overcoming barriers to scalability in variational quantum monte carlo.  In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, SC ’21, New York, NY, USA, 2021. Association for Computing Machinery.
  * [38] Mathieu Germain, Karol Gregor, Iain Murray, and Hugo Larochelle.  Made: Masked autoencoder for distribution estimation.  In Francis Bach and David Blei, editors, Proceedings of the 32nd International Conference on Machine Learning, volume 37 of Proceedings of Machine Learning Research, pages 881–889, Lille, France, 07–09 Jul 2015\. PMLR.
  * [39] Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova.  BERT: Pre-training of deep bidirectional transformers for language understanding.  In Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers), pages 4171–4186, Minneapolis, Minnesota, June 2019. Association for Computational Linguistics.
  * [40] John F. Kolen and Stefan C. Kremer.  Gradient Flow in Recurrent Nets: The Difficulty of Learning LongTerm Dependencies, pages 237–243.  2001\.
  * [41] Xiaohui Wang, Yang Wei, Ying Xiong, Guyue Huang, Xian Qian, Yufei Ding, Mingxuan Wang, and Lei Li.  Lightseq2: Accelerated training for transformer-based models on gpus.  In Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis, SC ’22. IEEE Press, 2022\.
  * [42] Xiaohui Wang, Ying Xiong, Yang Wei, Mingxuan Wang, and Lei Li.  LightSeq: A high performance inference library for transformers.  In Proceedings of the 2021 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies: Industry Papers (NAACL-HLT), pages 113–120. Association for Computational Linguistics, June 2021.
  * [43] Jiarui Fang, Yang Yu, Chengduo Zhao, and Jie Zhou.  Turbotransformers: An efficient gpu serving system for transformer models.  In Proceedings of the 26th ACM SIGPLAN Symposium on Principles and Practice of Parallel Programming, PPoPP ’21, page 389–402, New York, NY, USA, 2021. Association for Computing Machinery.


# Directionality between driven-dissipative resonators

C. A. Downing<EMAIL_ADDRESS>Department of Physics and Astronomy,
University of Exeter, Exeter EX4 4QL, United Kingdom  T. J. Sturges  Institute
of Theoretical Solid State Physics, Karlsruhe Institute of Technology (KIT),
D-76131 Karlsruhe, Germany

###### Abstract

Abstract
The notion of nonreciprocity, in essence when going forwards is different from
going backwards, emerges in all branches of physics from cosmology to
electromagnetism. Intriguingly, the breakdown of reciprocity is typically
associated with extraordinary phenomena, which may be readily capitalized on
in the design of (for example) nontrivial electromagnetic devices when Lorentz
reciprocity is broken. However, in order to enable the exploitation of
nonreciprocal-like effects in the next generation of quantum technologies,
basic quantum optical theories are required. Here we present a versatile model
describing a pair of driven-dissipative quantum resonators, where the relative
phase difference between the coherent and incoherent couplings induces an
asymmetry. The interplay between the diverse dissipative landscape – which
encompasses both intrinsic losses and dissipative couplings – and the coherent
interactions leads to some remarkable consequences including highly
directional (or even one-way) energy transport. Our work proffers the
tantalizing prospect of observing dissipation-induced quantum directionality
in areas like photonics or cavity magnonics (spin waves), which may aid the
design of unconventional nanoscopic devices.

Introduction
In perhaps one of the most profound conversations in Confucianism, Tsze-kung
asks “Is there one word which may serve as a rule of practice for all one’s
life?” The Master replies “Is not shu [reciprocity] such a word? What you do
not want done to yourself, do not do to others” Legge1971 .

In physics, reciprocity is a famous hallmark across the entire discipline,
from the equal and opposite quality of Newton’s third law of motion
Feynman1989  to the Lorentz reciprocity in electromagnetism, which guarantees
the same response when the source and receiver are interchanged Caloz2018 .
Introducing various types of asymmetries can lead to a treasure trove of
curiosities, which in optics precipitated the exciting subfields of chiral
plasmonics Hentschel2017 ; Andrews2018  and chiral quantum optics Lodahl2017 ;
Downing2019 ; Sanchez2019 . Within device physics, it has already been shown
that exploiting nonreciprocity can give rise to practical applications like
high quality factor, large bandwidth devices Tsakmakidis2017  (which are
predicated upon the induced asymmetric transport properties), optical
isolators Jalas2013  (where the governing scattering matrix is inherently
asymmetric), and even magnetic diodes Prat2018  (where it was demonstrated
that the magnetic coupling between two coils above a conductor - moving with
constant velocity - may become asymmetric, leading to a diode for magnetic
fields).

Figure 1:  A sketch of a pair of driven-dissipative resonators. The first
resonator (red pillar) is driven by a laser with amplitude $\Omega$, while the
second resonator (green pillar) is undriven [cf. Eq. (Directionality between
driven-dissipative resonators)]. The coherent coupling (of magnitude $g$ and
phase $\theta$) between the resonators is represented by the gray rod, while
the dissipative coupling (of magnitude $\Gamma$ and phase $\phi$) is mediated
by the common bath (blue disk) [cf. Eq. (Directionality between driven-
dissipative resonators)]. The individual loss $\gamma$ of each resonator is
associated with the red and green disks.

Here we investigate a basic quantum optical model, namely a pair of driven-
dissipative resonators, with a view to inducing asymmetric behaviour
Metelmann2015 ; Metelmann2017 ; Kamal2017 ; Downing2020 . We consider two
oscillators, which are in general coupled both coherently and incoherently,
where the first resonator is additionally coherently driven by a laser as
sketched in Fig. 1. The admixture between the coherent and incoherent
couplings, which are in general complex quantities, has profound consequences
for directionality in the system. There are four principle coupling regimes of
our model: (i) coherent coupling, where the direct hopping between the
resonators dominates (gray rod in the figure), (ii) dissipative coupling,
where incoherent coupling between the resonators is of primary importance (as
mediated by their common heat bath, the blue disk in the figure), (iii)
unidirectional coupling, where the coupling between the resonators is
completely one-way, and (iv) asymmetric coupling, where the mixture of
coherent and incoherent coupling leads to asymmetries in the interactions
between the resonator pair [such that this case generalizes the more extreme
limiting case (iii)]. In what follows, we provide a simple analysis of the
population dynamics in these regimes.

Figure 2:  The populations of a pair of driven-dissipative resonators. First
column: sketches of the coherent, dissipative and unidirectional coupling
regimes. Second column: the scaled steady state populations
$\lim_{t\to\infty}\langle
b_{n}^{\dagger}b_{n}\rangle\times\left(\gamma/\Omega\right)^{2}$ of the pair
(red and green lines), as a function of the coherent coupling strength $g$ or
the dissipative coupling strength $\Gamma$. The population imbalance $\Delta$
is shown as thick orange lines [cf. Eq. (8)]. Third column: the dynamic
populations $\langle b_{n}^{\dagger}b_{n}\rangle$ of the pair, as a function
of time $t$ (in units of the inverse loss rate $\gamma^{-1}$). In the figure,
the first resonator is driven with an amplitude $\Omega=\gamma/10$, in panel
(c) $g=2\gamma$ and in panels (f) and (i) $\Gamma=(4/5)\gamma$.

Model
The driven coupled oscillators model represented in Fig. 1 may be described by
the Hamiltonian $\hat{H}$ (with $\hbar=1$)

$\displaystyle\hat{H}=$
$\displaystyle~{}\omega_{\Delta}\left(b_{1}^{\dagger}b_{1}+b_{2}^{\dagger}b_{2}\right)+\Omega\left(b_{1}^{\dagger}+b_{1}\right)$
$\displaystyle+g\mathrm{e}^{\mathrm{i}\theta}b_{1}^{\dagger}b_{2}+g\mathrm{e}^{-\mathrm{i}\theta}b_{2}^{\dagger}b_{1},$
(1)

where the $n$-th oscillator sustains bosonic excitations created (destroyed)
by the operator $b_{n}^{\dagger}$ ($b_{n}$). The coherent oscillator-
oscillator coupling is of strength $g\geq 0$ and phase $\theta$. The first
oscillator is driven by a laser of amplitude $\Omega$, and the detunings
$\omega_{\Delta}$ arise in the chosen rotating reference frame of the laser
SuppInfo .

We include dissipation in the model through a quantum master equation, which
describes the time evolution of the density matrix $\rho$ of the system via
SuppInfo ; Gardiner2014 ; Dung2002

$\displaystyle\partial_{t}\rho=$ $\displaystyle~{}~{}\mathrm{i}[\rho,\hat{H}]$
$\displaystyle+\sum_{n=1,2}\frac{\gamma}{2}\left(2b_{n}\rho
b_{n}^{\dagger}-b_{n}^{\dagger}b_{n}\rho-\rho b_{n}^{\dagger}b_{n}\right)$
$\displaystyle+\frac{\Gamma\mathrm{e}^{\mathrm{i}\phi}}{2}\left(2b_{2}\rho
b_{1}^{\dagger}-b_{1}^{\dagger}b_{2}\rho-\rho b_{1}^{\dagger}b_{2}\right)$
$\displaystyle+\frac{\Gamma\mathrm{e}^{-\mathrm{i}\phi}}{2}\left(2b_{1}\rho
b_{2}^{\dagger}-b_{2}^{\dagger}b_{1}\rho-\rho b_{2}^{\dagger}b_{1}\right).$
(2)

The first line of Eq. (Directionality between driven-dissipative resonators)
is the von Neumann equation, describing the unitary evolution of the closed
system as governed by the Hamiltonian $\hat{H}$ of Eq. (Directionality between
driven-dissipative resonators). The non-unitary evolution is captured by the
three lower lines of Eq. (Directionality between driven-dissipative
resonators), which are written as Lindblad terms and describe the open quantum
system sketched in Fig. 1. In particular, the second line on the right-hand-
side of Eq. (Directionality between driven-dissipative resonators) accounts
for the intrinsic loss $\gamma$ of each oscillator. The third and fourth lines
of Eq. (Directionality between driven-dissipative resonators) track the
dissipative (incoherent) coupling between the pair of oscillators due to their
shared heat bath, which can generally be regarded as a complex quantity of
magnitude $\Gamma$ and phase $\phi$ (subject to the condition
$0\leq\Gamma\leq\gamma$).

The first moments $\langle b_{n}\rangle$ of the system described by Eq.
(Directionality between driven-dissipative resonators) and Eq. (Directionality
between driven-dissipative resonators) may be found by the following pair of
coupled first-order equations SuppInfo

$\mathrm{i}\partial_{t}\begin{pmatrix}\langle b_{1}\rangle\\\ \langle
b_{2}\rangle\end{pmatrix}=\begin{pmatrix}\omega_{\Delta}-\mathrm{i}\frac{\gamma}{2}&G_{-}\\\
G_{+}^{\ast}&\omega_{\Delta}-\mathrm{i}\frac{\gamma}{2}\end{pmatrix}\begin{pmatrix}\langle
b_{1}\rangle\\\ \langle b_{2}\rangle\end{pmatrix}+\begin{pmatrix}\Omega\\\
0\end{pmatrix},$ (3)

where we have introduced two generalized coupling constants $G_{+}$ and
$G_{-}$, hereby defined as

$G_{\pm}=g\mathrm{e}^{\mathrm{i}\theta}\pm\tfrac{1}{2}\Gamma\mathrm{e}^{\mathrm{i}\phi},$
(4)

which accounts for the admixture between the competing coherent and
dissipative couplings, including their magnitudes and phases. Most notably,
this analysis reveals that the oscillator-oscillator coupling may be
completely unidirectional, as was first noticed in the celebrated works on
cascaded quantum systems Gardiner1993 ; Carmichael1993 . When $G_{-}=0$ in Eq.
(3) one-way coupling in the rightwards direction ($\rightarrow$) arises, and
likewise when $G_{+}^{\ast}=0$ in Eq. (3) the coupling is completely leftwards
($\leftarrow$). Amongst the entire space of possible couplings, these two
special circumstances occur when the following conditions on the coupling
magnitudes and relative phases hold [cf. Eq. (4)]

$\Gamma=2g,\quad\quad\quad\quad\theta-\phi=\begin{cases}\tfrac{\pi}{2}&\left(\rightarrow\right)\\\
\tfrac{3\pi}{2}&\left(\leftarrow\right)\end{cases}.$ (5)

Away from these twin conditions for unidirectionality, the overall coupling is
generally asymmetric. Wonderfully, the basic theoretical model encapsulated by
Eq. (Directionality between driven-dissipative resonators) and Eq.
(Directionality between driven-dissipative resonators) may be realized in an
eclectic range of systems. For example, in spin-photon systems such as cavity
magnons with coherent and dissipative couplings Wang2020 ; Harder2021 , in
circuit-QED setups utilizing superconducting qubits Wendin2017 ; Blais2021 ,
in plasmonic epsilon-near-zero waveguides Issah2021 , in metallic nanoparticle
architectures exploiting plasmonic responses Sturges2020 , and with coupled
cavity-based photonic devices Metelmann2015 ; Azcona2021 . In what follows we
keep our discussion of the model general, keeping in mind specializations are
readily obtainable experimentally.

Coherent coupling
Let us first consider the simplest case of purely coherent coupling between
the resonators (so that $\Gamma=0$), as is represented in the sketch of Fig. 2
(a). We are interested in the populations $\langle
b_{n}^{\dagger}b_{n}\rangle$ of the resonators, which may be calculated from
the second moment analogue of Eq. (3), as discussed in Ref. SuppInfo .

At long time scales, the competition between the driving and dissipation leads
to a well-defined steady state, described by the analytic expressions SuppInfo

$\lim_{t\to\infty}\langle
b_{1}^{\dagger}b_{1}\rangle=\left(\frac{2\gamma\Omega}{\gamma^{2}+4g^{2}}\right)^{2},$
(6) $\lim_{t\to\infty}\langle
b_{2}^{\dagger}b_{2}\rangle=\left(\frac{4g\Omega}{\gamma^{2}+4g^{2}}\right)^{2}.$
(7)

These populations, scaled by $\left(\gamma/\Omega\right)^{2}$, are plotted as
a function of the coherent coupling strength $g$ for the first resonator (red
line) and second resonator (green line) in Fig. 2 (b). Clearly, since only the
first resonator is driven, its steady state population is bounded by its
maximum of $\left(2\Omega/\gamma\right)^{2}$ when $g\ll\gamma$, and decreases
to $\\{\gamma\Omega/\left(2g^{2}\right)\\}^{2}$ when $g\gg\gamma$. In these
limits, the second (and undriven) resonator population is zero when
$g\ll\gamma$ and $\left(\Omega/g\right)^{2}$ when $g\gg\gamma$, its maximum
population of $\left(\Omega/\gamma\right)^{2}$ is instead met when
$g=\gamma/2$. Hence the population imbalance between the resonators is a
useful quantity to describe the system, in the steady state it reads

$\Delta=\lim_{t\to\infty}\frac{\langle b_{1}^{\dagger}b_{1}\rangle-\langle
b_{2}^{\dagger}b_{2}\rangle}{\langle b_{1}^{\dagger}b_{1}\rangle+\langle
b_{2}^{\dagger}b_{2}\rangle}.$ (8)

In this coherent coupling regime, $\Delta$ is a sign-changing quantity, which
is explicitly given by

$\Delta=1-\frac{8g^{2}}{\gamma^{2}+4g^{2}},$ (9)

which observes the bounds of $-1\leq\Delta\leq 1$, as displayed with the thick
orange line in Fig. 2 (b). The critical point of a completely balanced
populations across the resonators $\Delta=0$ is reached when $g=\gamma/2$, and
above this coupling strength the second resonator has a larger steady state
population despite being undriven.

The full dynamic populations $\langle b_{n}^{\dagger}b_{n}\rangle$ of the
coupled resonators are given by the exact equations SuppInfo

$\displaystyle\langle b_{1}^{\dagger}b_{1}\rangle=$
$\displaystyle\left(\frac{2\gamma\Omega}{\gamma^{2}+4g^{2}}\right)^{2}+2\left(\frac{2\gamma\Omega}{\gamma^{2}+4g^{2}}\right)^{2}\biggl{\\{}2g\sin\left(gt\right)-\gamma\cos\left(gt\right)\biggl{\\}}\mathrm{e}^{-\frac{\gamma
t}{2}}$ (10)
$\displaystyle+\biggl{\\{}\Bigr{[}\gamma^{2}+4g^{2}\Bigr{]}\Bigr{[}\gamma^{2}+4g^{2}+4\Omega^{2}\Bigr{]}+\Bigr{[}\left(\gamma^{2}+4g^{2}\right)^{2}+4\Omega^{2}\left(\gamma^{2}-4g^{2}\right)\Bigr{]}\cos\left(2gt\right)-16g\gamma\Omega^{2}\sin\left(2gt\right)\biggl{\\}}\frac{\mathrm{e}^{-\gamma
t}}{2\left(\gamma^{2}+4g^{2}\right)^{2}},$ $\displaystyle\langle
b_{2}^{\dagger}b_{2}\rangle=$
$\displaystyle\left(\frac{4g\Omega}{\gamma^{2}+4g^{2}}\right)^{2}-\left(\frac{4g\Omega}{\gamma^{2}+4g^{2}}\right)^{2}\biggl{\\{}2g\cos\left(gt\right)+\gamma\sin\left(gt\right)\biggl{\\}}\mathrm{e}^{-\frac{\gamma
t}{2}}$ (11)
$\displaystyle+\biggl{\\{}\Bigr{[}\gamma^{2}+4g^{2}\Bigr{]}\Bigr{[}\gamma^{2}+4g^{2}+4\Omega^{2}\Bigr{]}-\Bigr{[}\left(\gamma^{2}+4g^{2}\right)^{2}+4\Omega^{2}\left(\gamma^{2}-4g^{2}\right)\Bigr{]}\cos\left(2gt\right)+16g\gamma\Omega^{2}\sin\left(2gt\right)\biggl{\\}}\frac{\mathrm{e}^{-\gamma
t}}{2\left(\gamma^{2}+4g^{2}\right)^{2}},$

which are plotted in Fig. 2 (c) as a function of time (for the example case of
$g=2\gamma$). The characteristic damped Rabi oscillations are shown, and the
two different time constants appearing in the above expressions ($1/\gamma$
and $2/\gamma$) is characteristic of coherently driven systems. This sets out
the stall for the most well known coupling regime.

Dissipative coupling
When the coherent coupling is negligible (that is, $g=0$), the system is in
the dissipative coupling regime. This setup is sketched in Fig. 2 (d), where
the blue disk represents the common heat bath enabling the incoherent
coupling. In the steady state, the resonator populations are described by the
simple forms SuppInfo

$\lim_{t\to\infty}\langle
b_{1}^{\dagger}b_{1}\rangle=\left(\frac{2\gamma\Omega}{\gamma^{2}-\Gamma^{2}}\right)^{2},$
(12) $\lim_{t\to\infty}\langle
b_{2}^{\dagger}b_{2}\rangle=\left(\frac{2\Gamma\Omega}{\gamma^{2}-\Gamma^{2}}\right)^{2}.$
(13)

In the weak dissipative coupling limit $\Gamma\ll\gamma$ of course only the
first, driven resonator is populated, with $\left(2\Omega/\gamma\right)^{2}$.
In the opposing strong dissipative coupling $\Gamma\to\gamma$ limit the
bosonic nature of the resonators becomes readily apparent, since both
populations tend towards infinity, as is shown in in Fig. 2 (e). The
population imbalance [cf. Eq. (8)] is always singled-signed and reads

$\Delta=1-\frac{2\Gamma^{2}}{\Gamma^{2}+\gamma^{2}},$ (14)

which exposes the necessarily non-negative bounds of $0\leq\Delta\leq 1$, as
show by the thick orange line in Fig. 2 (e). This in stark contrast to the
coherent coupling regime, which allows for positive and negative imbalances
[cf. panel (b)].

The time-dependent populations are described by the following expressions
SuppInfo

$\displaystyle\langle b_{1}^{\dagger}b_{1}\rangle=$
$\displaystyle\left(\frac{2\gamma\Omega}{\gamma^{2}-\Gamma^{2}}\right)^{2}+\left(\cosh\left(\Gamma
t\right)+1+\frac{4\Omega^{2}}{\gamma^{2}-\Gamma^{2}}\right)\frac{\mathrm{e}^{-\gamma
t}}{2}+\Omega^{2}\biggl{\\{}\frac{\mathrm{e}^{-\left(\gamma+\Gamma\right)t}}{\left(\gamma+\Gamma\right)^{2}}+\frac{\mathrm{e}^{-\left(\gamma-\Gamma\right)t}}{\left(\gamma-\Gamma\right)^{2}}\biggl{\\}}$
$\displaystyle-\frac{4\gamma\Omega^{2}}{\gamma^{2}-\Gamma^{2}}\biggl{\\{}\frac{\mathrm{e}^{-\frac{\left(\gamma+\Gamma\right)t}{2}}}{\gamma+\Gamma}+\frac{\mathrm{e}^{-\frac{\left(\gamma-\Gamma\right)t}{2}}}{\gamma-\Gamma}\biggl{\\}},$
(15) $\displaystyle\langle b_{2}^{\dagger}b_{2}\rangle=$
$\displaystyle\left(\frac{2\Gamma\Omega}{\gamma^{2}-\Gamma^{2}}\right)^{2}+\left(\cosh\left(\Gamma
t\right)-1-\frac{4\Omega^{2}}{\gamma^{2}-\Gamma^{2}}\right)\frac{\mathrm{e}^{-\gamma
t}}{2}+\Omega^{2}\biggl{\\{}\frac{\mathrm{e}^{-\left(\gamma+\Gamma\right)t}}{\left(\gamma+\Gamma\right)^{2}}+\frac{\mathrm{e}^{-\left(\gamma-\Gamma\right)t}}{\left(\gamma-\Gamma\right)^{2}}\biggl{\\}}$
$\displaystyle+\frac{4\Gamma\Omega^{2}}{\gamma^{2}-\Gamma^{2}}\biggl{\\{}\frac{\mathrm{e}^{-\frac{\left(\gamma+\Gamma\right)t}{2}}}{\gamma+\Gamma}-\frac{\mathrm{e}^{-\frac{\left(\gamma-\Gamma\right)t}{2}}}{\gamma-\Gamma}\biggl{\\}},$
(16)

as displayed in Fig. 2 (f) for the example case of reasonably strong
dissipative coupling $\Gamma=(4/5)\gamma$. A hallmark of this coupling regime
is the lack of any Rabi oscillations due to the absence of any coherent
coupling, and the supremacy in population of the first resonator for any value
of the dissipative coupling strength $\Gamma$, as suggested by Eq. (14).

Unidirectional coupling
The final special case of coupling that we shall consider is that of
unidirectional coupling, and in particular when the conditions of Eq. (5) are
met for the rightwards ($\rightarrow$) direction only, as represented by the
picture in Fig. 2 (g). The lack of backaction ensures that the first resonator
population is coupling independent (in this regime, the dissipative coupling
strength $\Gamma=2g$ is fixed) while the second resonator population is
enhanced due to the one-way nature of the interaction. The steady state
results are simply SuppInfo

$\lim_{t\to\infty}\langle
b_{1}^{\dagger}b_{1}\rangle=\left(\frac{2\Omega}{\gamma^{2}}\right)^{2},$ (17)
$\lim_{t\to\infty}\langle
b_{2}^{\dagger}b_{2}\rangle=\left(\frac{4\Gamma\Omega}{\gamma^{2}}\right)^{2},$
(18)

which are plotted in Fig. 2 (h) as a function of $\Gamma$. As must be the
case, the first resonator population (red line) is exactly that of a single
driven-dissipative oscillator SuppInfo , while the second resonator resonator
presents a quadratic scaling with the dissipative coupling $\Gamma$.
Therefore, the population imbalance [cf. Eq. (8)] is given by

$\Delta=1-\frac{8\Gamma^{2}}{\gamma^{2}+4\Gamma^{2}},$ (19)

as plotted as the thick orange line in Fig. 2 (h). Unlike in the dissipative
coupling regime, this population imbalance is a sign-changing quantity, being
bounded by $-3/5\leq\Delta\leq 1$. The critical point of $\Delta=0$ is reached
when $\Gamma=\gamma/2$, such that above this dissipative coupling strength the
undriven second resonator is more highly populated that the driven and first
resonator.

Figure 3:  Phase-dependent population imbalance of coupled driven-dissipative
resonators. The population imbalance $\Delta$ of the pair in the steady state
[cf. Eq. (24)], as a function of the relative phase $\theta-\phi$ between the
coherent and dissipative couplings, and the magnitude $\Gamma$ of the
dissipative coupling (in units of the loss rate $\gamma$). In the figure, the
magnitude of coherent coupling $g=\gamma/2$.

The dynamic populations are given by the compact analytical expressions
SuppInfo

$\displaystyle\langle b_{1}^{\dagger}b_{1}\rangle=$
$\displaystyle\left(\frac{2\Omega}{\gamma}\right)^{2}-2\left(\frac{2\Omega}{\gamma}\right)^{2}\mathrm{e}^{-\frac{\gamma
t}{2}}+\biggl{\\{}1+\left(\frac{2\Omega}{\gamma}\right)^{2}\biggl{\\}}\mathrm{e}^{-\gamma
t},$ (20) $\displaystyle\langle b_{2}^{\dagger}b_{2}\rangle=$
$\displaystyle\left(\frac{4\Gamma\Omega}{\gamma}\right)^{2}-\left(\frac{4\Gamma\Omega}{\gamma}\right)^{2}\left(2+\gamma
t\right)\mathrm{e}^{-\frac{\gamma
t}{2}}+\left(\frac{\Gamma}{\gamma}\right)^{2}\biggl{\\{}\left(\gamma
t\right)^{2}+\left(2+\gamma
t\right)^{2}\left(\frac{2\Omega}{\gamma}\right)^{2}\biggl{\\}}\mathrm{e}^{-\gamma
t},$ (21)

where Eq. (20) is exactly the form for a solitary driven resonator, as if the
second resonator was not there, due to the completely supressed backaction. We
plot the expressions of Eq. (20) and Eq. (21) in Fig. 2 (i) for the case of
the rather strong coupling $\Gamma=(4/5)\gamma$. This coupling arrangement
allows for the second resonator population (green line) to become dominant
after only a short timescale $t\sim 1/\gamma$, which is maintained through to
the steady state and thus evermore.

Asymmetric coupling
In general, the coupling encompassed by the model of Eq. (Directionality
between driven-dissipative resonators) is asymmetric - with the preceding
unidirectional case being the most extreme example. Most generally then (when
$g\neq 0$ and $\Gamma\neq 0$), the resonator steady states become dependent on
the relative phase $\theta-\phi$ as follows SuppInfo

$\lim_{t\to\infty}\langle
b_{1}^{\dagger}b_{1}\rangle=\frac{\left(2\gamma\Omega\right)^{2}}{16g^{4}+8g^{2}\gamma^{2}+\left(\gamma^{2}-\Gamma^{2}\right)^{2}+8g^{2}\Gamma^{2}\cos\left(2\left[\theta-\phi\right]\right)},$
(22) $\lim_{t\to\infty}\langle
b_{2}^{\dagger}b_{2}\rangle=\frac{4\Omega^{2}\left(4g^{2}+\Gamma^{2}+4g\Gamma\sin\left[\theta-\phi\right]\right)^{2}}{16g^{4}+8g^{2}\gamma^{2}+\left(\gamma^{2}-\Gamma^{2}\right)^{2}+8g^{2}\Gamma^{2}\cos\left(2\left[\theta-\phi\right]\right)},$
(23)

from which the analogous results in other, more specialized coupling regimes
may be derived. The population imbalance measure [cf. Eq. (8)] is given by the
rich expression

$\Delta=\frac{2\gamma^{2}}{16g^{4}+8g^{2}\gamma^{2}+\left(\gamma^{2}-\Gamma^{2}\right)^{2}+8g^{2}\Gamma^{2}\cos\left(2\left[\theta-\phi\right]\right)}-1,$
(24)

and is plotted in Fig. 3, as a function of the relative phase $\theta-\phi$
and the dimensionless dissipative coupling strength $\Gamma/\gamma$, for the
example case where the magnitude of the coherent coupling $g=\gamma/2$.
Clearly, the map of Fig. 3 exposes the importance of the phase as the
determiner of the asymmetry of the coupling, since modulating $\theta-\phi$
allows for either positive (yellow to red) or negative (cyan to blue) steady
state population imbalances while holding the magnitudes of all other
parameters constant. In particular, the red area around $\theta-\phi=\pi/2$
and blue region around $\theta-\phi=3\pi/2$ are plausible from the knowledge
of the unidirectional phase conditions of Eq. (5).

Figure 4:  Phase-dependent population dynamics of coupled driven-dissipative
resonators. The dynamic populations $\langle b_{n}^{\dagger}b_{n}\rangle$ of
the pair, as a function of time $t$ (in units of the inverse loss rate
$\gamma^{-1}$). The relative phase $\theta-\phi$ between the coherent and
dissipative couplings is increased from $0$ to $\pi/2$ to $3\pi/2$ upon
descending the column of panels. In the figure, the first resonator is driven
with an amplitude $\Omega=\gamma/10$ and the magnitudes of the coherent and
dissipative couplings are $g=\gamma/2$ and $\Gamma=\gamma$, so that the
magnitude unidirectional coupling condition $\Gamma=2g$ is met [cf. Eq. (5)].

The dynamical populations are shown in Fig. 4 for the case of $g=\gamma/2$ and
maximal dissipative coupling $\Gamma=\gamma$, so that the magnitude
unidirectional coupling condition is met [cf. Eq. (5)]. In panel (a) the
relative phase is zero ($\theta-\phi=0$) and the large incoherent-to-coherent
coupling ratio of $\Gamma/g=2$ ensures that Rabi cycles are not discernible
within the damped population cycles. In panel (b) the rightwards
($\rightarrow$) unidirectional phase condition $\theta-\phi=\pi/2$ is
fulfilled, such that the lack of backaction sees a purely exponential decay of
the first resonator population (red line) and strong enhancement of the second
resonator population (green line). Finally, in panel (c) the leftwards
($\leftarrow$) unidirectional phase condition $\theta-\phi=3\pi/2$ is
satisfied, such that the undriven second resonator is never populated, in the
most dramatic realization of the directionality of the coupled system.

Conclusions
We have studied a simple yet explanatorily powerful model of a pair of driven-
dissipative resonators with both coherent and incoherent couplings. Our theory
acts as a prototypical example of how dissipation-induced directionality may
arise in quantum optical systems, with dramatic implications. In particular,
we have shown how tailoring the relative magntiude and phase of the coherent
and dissipative coupling can lead to highly directional and even one-way
quantum transport. Our results provide perspectives for the quantum
engineering of coupled resonators Chang2018 , with applications for
directional devices such as isolators, circulators and quantum batteries
Scheucher2016 ; Barzanjeh2017 ; Roushan2017 .

Acknowledgments
Funding: CAD is supported by the Royal Society via a University Research
Fellowship (URF/R1/201158) and a Royal Society Research Grant (RGS/R1
/211220). TJS acknowledges funding from the Alexander von Humboldt Foundation.
Discussions: we thank A. I. Fernández-Domínguez and E. del Valle for
discussions. Data and materials availability: There is no data in this wholly
theoretical work. All necessary information is available in the manuscript and
the Supplementary Information SuppInfo .

ORCID
C. A. Downing: 0000-0002-0058-9746.
T. J. Sturges: 0000-0003-1320-2843.

## References

  * (1) J. Legge, Confucius: Confucian Analects, the Great Learning, and the Doctrine of the Mean (Dover, New York, 1971).
  * (2) R. P. Feynman, R. Leighton and M. Sands, The Feynman Lectures on Physics (Dover, New York, 1989).
  * (3) C. Caloz, A. Alu, S. Tretyakov, D. Sounas, K. Achouri, and Z.-L. Deck-Leger, Electromagnetic nonreciprocity, Phys. Rev. Applied 10, 047001 (2018).
  * (4) M. Hentschel, M. Schaferling, X. Duan, H. Giessen, and N. Liu, Chiral plasmonics, Sci. Adv. 3, e1602735 (2017).
  * (5) D. L. Andrews, Quantum formulation for nanoscale optical and material chirality: symmetry issues, space and time parity, and observables, J. Opt. 20, 033003 (2018).
  * (6) P. Lodahl, S. Mahmoodian, S. Stobbe, A. Rauschenbeutel, P. Schneeweiss, J. Volz, H. Pichler, and P. Zoller, Chiral quantum optics, Nature 541, 473 (2017).
  * (7) C. A. Downing, J. C. Lopez Carreno, F. P. Laussy, E. del Valle, A. I. Fernandez-Dominguez, Quasichiral interactions between quantum emitters at the nanoscale, Phys. Rev. Lett. 122, 057401 (2019).
  * (8) E. Sanchez-Burillo, C. Wan, D. Zueco, and A. Gonzalez-Tudela, Chiral quantum optics in photonic sawtooth lattices, Phys. Rev. Research 2, 023003 (2020).
  * (9) K. L. Tsakmakidis, L. Shen, S. A. Schulz, X. Zheng, J. Upham, X. Deng, H. Altug, A. F. Vakakis, R. W. Boyd, Breaking Lorentz reciprocity to overcome the time-bandwidth limit in physics and engineering, Science 356, 1260 (2017).
  * (10) D. Jalas, A. Petrov, M. Eich, W. Freude, S. Fan, Z. Yu, R. Baets, M. Popovic, A. Melloni, J. D. Joannopoulos, M. Vanwolleghem, C. R. Doerr, and H. Renner, What is - and what is not - an optical isolator, Nature Photon 7, 579 (2013).
  * (11) J. Prat-Camps, P. Maurer, G. Kirchmair, and O. Romero-Isart, Circumventing magnetostatic reciprocity: a diode for magnetic fields, Phys. Rev. Lett. 121, 213903 (2018).
  * (12) A. Metelmann and A. A. Clerk, Nonreciprocal photon transmission and amplification via reservoir engineering, Phys. Rev. X 5, 021025 (2015).
  * (13) A. Metelmann and A. A. Clerk, Nonreciprocal quantum interactions and devices via autonomous feedforward, Phys. Rev. A 95, 013837 (2017).
  * (14) A. Kamal and A. Metelmann, Minimal models for nonreciprocal amplification using biharmonic drives, Phys. Rev. Applied 7, 034031 (2017).
  * (15) C. A. Downing, J. C. López Carreño, A. I. Fernández-Domínguez, and E del Valle, Asymmetric coupling between two quantum emitters, Phys. Rev. A 102, 013723 (2020).
  * (16) Please see the Supplementary Information for the background theory supporting the results reported in the main text. Included is a systematic review of a single driven-dissipative oscillator, before a detailed treatment of a pair of coupled driven-dissipative oscillators. It includes the Refs Gardiner2014 ; Metelmann2015 ; Downing2019 ; Downing2020 .
  * (17) C. Gardiner and P. Zoller, The Quantum World of Ultra-Cold Atoms and Light, Book I: Foundations of Quantum Optics (Imperial College Press, London, 2014).
  * (18) H. T. Dung, L. Knoll, and D. -G. Welsch, Resonant dipole-dipole interaction in the presence of dispersing and absorbing surroundings, Phys. Rev. A 66, 063810 (2002).
  * (19) C. W. Gardiner, Driving a quantum system with the output field from another driven quantum system, Phys. Rev. Lett. 70, 2269 (1993).
  * (20) H. J. Carmichael, Quantum trajectory theory for cascaded open systems, Phys. Rev. Lett. 70, 2273 (1993).
  * (21) Y.-P. Wang, and C.-M. Hu, Dissipative couplings in cavity magnonics, J. Appl. Phys. 127, 130901 (2020).
  * (22) M. Harder, B. Yao, Y. Gui, and C.-M. Hu, Coherent and dissipative cavity magnonics, J. Appl. Phys. 129, 201101 (2021).
  * (23) G. Wendin, Quantum information processing with superconducting circuits: a review, Rep. Prog. Phys. 80, 106001 (2017).
  * (24) A. Blais, A. L. Grimsmo, S. M. Girvin, and A. Wallraff, Circuit quantum electrodynamics, Rev. Mod. Phys. 93, 025005 (2021).
  * (25) I. Issah and H. Caglayan, Qubit-qubit entanglement mediated by epsilon-near-zero waveguide reservoirs, arxiv.
  * (26) T. J. Sturges, T. Repan, C. A. Downing, C. Rockstuhl, M. Stobinska, Extreme renormalisations of dimer eigenmodes by strong light-matter coupling, New J. Phys. 22, 103001 (2020).
  * (27) P. Martinez Azcona, and C. A. Downing, Doublons, topology and interactions in a one-dimensional lattice, Sci. Rep. 11, 12540 (2021).
  * (28) D. E. Chang, J. S. Douglas, A. Gonzalez-Tudela, C.-L. Hung, and H. J. Kimble, Colloquium: Quantum matter built from nanoscopic lattices of atoms and photons, Rev. Mod. Phys. 90, 031002 (2018).
  * (29) M. Scheucher, A. Hilico, E. Will, J. Volz, and A. Rauschenbeutel, Quantum optical circulator controlled by a single chirally coupled atom, Science 354, 1577 (2016).
  * (30) S. Barzanjeh, M. Wulf, M. Peruzzo, M. Kalaee, P. B. Dieterle, O. Painter, and J. M. Fink, Mechanical on-chip microwave circulator, Nat. Commun 8, 953 (2017).
  * (31) P. Roushan, 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. J. J. O’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).