# Dissipative variational quantum algorithms for Gibbs state preparation

Yigal Ilin1<EMAIL_ADDRESS>Itai Arad2 1Andrew and Erna Viterbi
Department of Electrical and Computer Engineering, Technion – Israel Institute
of Technology, 3200003, Haifa, Israel  2Centre for Quantum Technologies,
National University of Singapore, 117543 Singapore, Singapore

###### Abstract

In recent years, variational quantum algorithms (VQAs) have gained significant
attention due to their adaptability and efficiency on near-term quantum
hardware. They have shown potential in a variety of tasks, including linear
algebra, search problems, Gibbs and ground state preparation. Nevertheless,
the presence of noise in current day quantum hardware, severely limits their
performance. In this work, we introduce _dissipative variational quantum
algorithms_ (D-VQAs) by incorporating dissipative operations, such as qubit
RESET and stochastic gates, as an intrinsic part of a variational quantum
circuit. We argue that such dissipative variational algorithms posses some
natural resilience to dissipative noise. We demonstrate how such algorithms
can prepare Gibbs states over a wide range of quantum many-body Hamiltonians
and temperatures, while significantly reducing errors due to both coherent and
non-coherent noise. An additional advantage of our approach is that no ancilla
qubits are need. Our results highlight the potential of D-VQAs to enhance the
robustness and accuracy of variational quantum computations on NISQ devices.

##  I Introduction

In recent years, variational quantum algorithms (VQAs) have gathered
significant research interest due to their adaptability and efficiency on
near-term quantum hardware [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]. A typical VQA
involves preparing a trial state on a quantum computer and calculating a cost
function based on measurements taken from that state [11, 12, 13, 10, 1]. As
of today, VQAs are one of our best candidates for a useful early quantum
advantage. Firstly, the framework is very rich and adaptive, and can be
applied to many practically important problems. VQAs have been designed for
various linear algebra problems [14, 15, 16], search problems [17],
determining the ground state of a given Hamiltonian [18], singular value
decomposition [19], fidelity estimation [20] and many other tasks [13, 8, 21,
11]. Secondly, the VQA framework leverages a significant portion of the
computation to a classical side, allowing for much shallower quantum circuits
that are crucial for our current day noisy devices [22, 23]. Finally, the
variational nature of these algorithms provides some resilience against
coherent errors, which are prevalent in today’s quantum hardware [24, 22, 25,
26, 27].

Despite the above strengths, VQAs have not yet demonstrated a convincing
quantum advantage. Arguably, the main obstacle are the noise levels present in
current day quantum hardware (although other challenges exist even in
noiseless devices, see Ref. [28] and Refs. [29, 30, 31]). The accumulation of
noise in a VQA circuit degrades its performance. Specifically, non-coherent
noise, such as decoherence or amplitude damping, poses the greatest challenge,
as coherent noise can largely be mitigated by the variational ansatz. In such
cases, one must resort to error mitigation techniques (see Refs. [32, 33, 34,
35, 36, 37]) to rectify the effects of noise, which usually come at an
exponential cost, thereby again limiting the performance of the VQA, or wait
for error correction to become viable.

In this study, we expand the capabilities of variational quantum algorithms by
incorporating _dissipative_ operations in the variational circuits. We call
the resulting framework _dissipative VQA_ (D-VQA). Our main claim is that the
introduction of these new elements can alleviate the effect of non-coherent
noise, much like variational unitaries can fight coherent errors. A simple
example of a dissipative operation within a quantum circuit is a mid-circuit
measurement, potentially followed by a qubit reset. The inclusion of such
operations in quantum circuits has been investigated in various contexts,
including local quantum channel learning [38], measurement-induced
entanglement phase transitions [39], error mitigation through postselection
[40], quantum steering for state preparation [41], non-equilibrium phase
transitions [42], and others [43, 44, 45].

Here, we investigate the efficacy of dissipative VQA, where we introduce a
variational dissipative operation, denoted as $\mathcal{R}_{i}$, alongside
unitary gates. Our dissipative operation includes both a stochastic element (a
probabilistic gate), together with a RESET operation. Specifically, it is a
parametrized gate, denoted by $\mathcal{R}_{i}(p,{\bm{\phi}})$, which resets
qubit $i$ to the pure state
${|{{\bm{\phi}}}\rangle\\!\,\langle{{\bm{\phi}}}|}$ with probability $p$. We
demonstrate that incorporating such gates offers two significant advantages:
1) The inclusion of dissipative operations allows for the creation of mixed
states, eliminating the need for ancilla qubits when the trial state of the
VQA is mixed. 2) D-VQA circuits can mitigate some of the non-coherent errors,
similar to how unitary VQAs mitigate coherent errors, resulting in better
performance in the presence of noise. This resilience is particularly good
when targeting states with low purity.

We demonstrate the advantages of our D-VQA ansatz through a simple toy model
and a set of classical simulations of variational Gibbs state preparation with
high fidelity in both noisy and noiseless scenarios across a wide range of
quantum many-body Hamiltonians.

The structure of this paper is as follows: In Sec. II, we provide a
comprehensive background on Gibbs state preparation on quantum devices. In
Sec. III, we introduce a single-qubit toy model of a quantum circuit used in a
D-VQA scheme, demonstrating its noise resilience against non-coherent errors.
We also present our global ansatz, which employs the $\mathcal{R}$ gate and
general 2-local unitaries. Sec. IV outlines the numerical optimization scheme
and defines the noise model used. Sec. V presents the numerical results.
Finally, in Sec. VI we present our conclusions and discuss possible future
research directions.

##  II Background

Throughout this work we consider the problem of preparing a Gibbs state on a
quantum computer. Gibbs state preparation is a central problem in quantum
computation and information of both theoretical and practical importance.
Specifically, Gibbs states are used in quantum machine learning [46, 47, 48,
49, 50], quantum simulation [51, 52], and quantum optimization [53, 54], and
for training quantum Boltzmann machines by sampling from a well-prepared Gibbs
state [55, 56, 57].

Given a Hamiltonian $H$, its Gibbs state $\rho_{G}$ at inverse temperature
$\beta=1/k_{B}T$ is given by

$\displaystyle\rho_{G}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\frac{1}{Z}e^{-\beta
H},$ (1)

where
$Z\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\operatorname{Tr}(\exp{(-\beta
H)})$ is the partition function.

There are various approaches for preparing a Gibbs state on a quantum
computer, which can be broadly categorized into two types. The first type
includes non-variational algorithms that steer the system toward the target
Gibbs state. These algorithms may simulate the physical process of
thermalization, as seen in Refs. [58, 59, 60, 61, 62, 63, 64, 65], utilize
imaginary time evolution like in Refs. [66, 67], or use Monte Carlo style
Gibbs samplers as in Refs. [68, 69]. Other methods include thermal shadow
tomography [70], sampling from the Gibbs state using cluster expansions [71],
dissipative sampling [72], or various other approaches [73, 74, 75, 76, 77].

The second type consists of variational algorithms, where a classical
optimization routine iteratively updates the parameters of the variational
quantum circuit based on a cost function that measures the proximity of the
system state to the desired Gibbs state [4, 78, 79, 2]. These algorithms may
leverage approaches similar to non-variational algorithms, such as
thermalization [80, 26], imaginary time evolution [26, 81], or Monte Carlo
style methods [82], among others [3, 83]. Additionally, variational algorithms
have been successfully combined with classical neural network architectures to
enhance the classical component of Gibbs sampling algorithms, as shown in
Refs. [84, 85].

Some variational algorithms are designed specifically for NISQ-type computers.
These algorithms aim to prepare a Gibbs state using minimal quantum resources,
such as the number of gates and measurements, often relying on heuristic
approaches [84, 85, 86, 87, 88, 89, 90, 91, 92, 93]. Most of these NISQ-type
algorithms require ancillary qubits coupled to the system where the thermal
equilibrium state is prepared.

The natural cost function for the variational algorithm is the Helmholtz free-
energy function
$\mathcal{F}(\rho)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\operatorname{Tr}(H\rho)-\beta^{-1}S(\rho)$,
as it is minimized by the Gibbs state $\rho_{G}$ [94]. While the energy term
$\operatorname{Tr}(H\rho)$ is a local observable, which can be estimated
efficiently on a quantum computer, the von-Neumman entropy term
$S(\rho)\stackrel{{\scriptstyle\mathrm{def}}}{{=}}-\operatorname{Tr}(\rho\log\rho)$
is highly non-local, and cannot be estimated efficiently. Therefore, in
variational algorithms for the Gibbs states one has to make additional
assumptions and/or approximations when using the free energy cost function
[95, 96, 86, 87, 97, 88, 89, 90, 91, 92, 93, 98], or resort to a different
cost function altogether [99].

The main purpose of this work is to investigate which Gibbs states can be
prepared with dissipative variational circuits and how the accuracy changes in
the presence of noise. Specifically, we aim to study the _expressibility_ of
the proposed ansatz. We shall therefore ignore the question of efficiency in
the cost function and use the _infidelity_ of the resultant state with the
ideal Gibbs state as a cost function. Nevertheless, our ansatz can readily
accommodate other loss function choices, including those with gradients
calculated on quantum hardware, such as in Refs. [99, 89].

##  III Dissipative circuits and noise resilience

In this section we describe the central building block of our circuit, which
is the dissipative single-qubit $\mathcal{R}$ gate. We study its noise
resilience properties through a simple one qubit example, and then describe
our global ansatz that uses the $\mathcal{R}$ gate and general 2-local
unitaries.

###  III.1 Dissipative gates and noise resilience

Our variational circuit is based on a specific dissipative gate, which is a
single-qubit, probabilistic gate that combines a RESET gate and a general
1-qubit rotation. For brevity, we shall refer to it as the $\mathcal{R}$ gate.
When acting on qubit $i$, its action is described by the non-unital quantum
channel

$\displaystyle\mathcal{R}_{i}(p,{\bm{\phi}})[\rho]$
$\displaystyle\stackrel{{\scriptstyle\mathrm{def}}}{{=}}(1-p)\rho+p{|{{\bm{\phi}}}\rangle\\!\,\langle{{\bm{\phi}}}|}\cdot\operatorname{Tr}_{i}(\rho).$
(2)

$\mathcal{R}_{i}(p,{\bm{\phi}})$ resets qubit $i$ to the state
${|{{\bm{\phi}}}\rangle\\!\,\langle{{\bm{\phi}}}|}$ with probability
$p\in[0,1]$, and with probability $1-p$ leaves the state unchanged. Since
${|{{\bm{\phi}}}\rangle\\!\,\langle{{\bm{\phi}}}|}$ is a pure state on the
Bloch sphere, it is described by two angles. Together with $p$, it uses three
parameters.

One advantage of variational dissipative quantum circuits over unitary
circuits is their partial resilience to dissipative noise. To illustrate this,
we consider a toy model system with one qubit, assuming a simple depolarizing
noise model.

Let ${\bm{r}}$ be a vector in the Bloch ball, and
${\bm{\sigma}}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}(X,Y,Z)$ be the vector
of Pauli matrices. Consider a situation in which the ideal circuit applies a
single qubit unitary gate $U$ to a mixed input state
$\rho_{0}=\frac{1}{2}(\mathbb{I}+|{\bm{r}}|\cdot Z)$, and takes it to the
ideal output state $\rho_{1}$:

$\displaystyle\rho_{1}=\mathcal{U}[\rho_{0}]=U\rho_{0}U^{\dagger}=\frac{1}{2}(\mathbb{I}+{\bm{r}}\cdot{\bm{\sigma}}),$

where $\mathcal{U}$ is the channel corresponding to the unitary $U$, i.e.,
$\mathcal{U}[\rho]=U\rho U^{\dagger}$. For simplicity, we assume that a
depolarizing noise channel, given by

$\displaystyle\mathcal{N}[\rho]\stackrel{{\scriptstyle\mathrm{def}}}{{=}}(1-\lambda)\rho+\lambda\operatorname{Tr}(\rho)\cdot\frac{\mathbb{I}}{2},$
(3)

acts on the qubit after each gate. In such case, we can look for a unitary $V$
with a corresponding channel $\mathcal{V}$ such that
$\mathcal{N}\circ\mathcal{V}(\rho_{0})$ is as close as possible to the ideal
state $\rho_{1}$. For any one-qubit state
$\tau=\frac{1}{2}(\mathbb{I}+{\bm{n}}\cdot{\bm{\sigma}})$, the action of the
depolarizing channel $\mathcal{N}$ is to shrink its Bloch radius, i.e.,
$\tau\to\mathcal{N}[\tau]=\frac{1}{2}\big{(}\mathbb{I}+(1-\lambda){\bm{n}}\cdot{\bm{\sigma}}\big{)}$.
Therefore, it is easy to see that the optimal unitary $V$ will be $V=U$, and
the optimal noisy output $\rho^{\prime}_{1}$ is

$\displaystyle\rho^{\prime}_{1}=\mathcal{N}\circ\mathcal{U}[\rho_{0}]=\mathcal{N}[\rho_{1}]=\frac{1}{2}\big{(}\mathbb{I}+(1-\lambda){\bm{r}}\cdot{\bm{\sigma}}\big{)}.$

In such case, the optimized noisy state suffers an error of

$\displaystyle{\|\rho_{1}-\rho^{\prime}_{1}\|}_{1}=\big{\|}\frac{1}{2}\lambda{\bm{r}}\cdot{\bm{\sigma}}\big{\|}_{1}=\lambda\cdot|{\bm{r}}|.$

Now assume that we can use the dissipative $\mathcal{R}$ gate, and apply a
circuit of the form $\mathcal{R}\circ\mathcal{U}$. Together with the noise,
our circuit will apply the channel
$\mathcal{N}\circ\mathcal{R}\circ\mathcal{N}\circ\mathcal{U}$ so that

$\displaystyle\tilde{\rho}_{1}=\mathcal{N}\circ\mathcal{R}\circ\mathcal{N}\circ\mathcal{U}(\rho_{0})=\mathcal{N}\circ\mathcal{R}\big{[}\frac{1}{2}\big{(}\mathbb{I}+(1-\lambda){\bm{r}}\cdot{\bm{\sigma}}\big{)}\big{]}.$

To optimize the circuit, we first fix the ${\bm{\phi}}$ parameter of the
$\mathcal{R}$ gate to be such that
${|{{\bm{\phi}}}\rangle\\!\,\langle{{\bm{\phi}}}|}=\frac{1}{2}(\mathbb{I}+\hat{{\bm{r}}}\cdot{\bm{\sigma}})$,
where
$\hat{{\bm{r}}}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}{\bm{r}}/|{\bm{r}}|$,
and leave $p$ unfixed. A simple calculation shows that

$\displaystyle\mathcal{R}\big{[}\frac{1}{2}\big{(}\mathbb{I}+(1-\lambda){\bm{r}}\cdot{\bm{\sigma}}\big{)}\big{]}$
$\displaystyle=\frac{1}{2}\big{[}\mathbb{I}+\big{(}(1-p)(1-\lambda)+p/|{\bm{r}}|\big{)}{\bm{r}}\cdot{\bm{\sigma}}\big{]},$

and therefore,

$\displaystyle\tilde{\rho}_{1}$
$\displaystyle=\frac{1}{2}\big{[}\mathbb{I}+(1-\lambda)\big{(}(1-p)(1-\lambda)+p/|{\bm{r}}|\big{)}{\bm{r}}\cdot{\bm{\sigma}}\big{]}.$

It follows that as long as we can choose $p\in[0,1]$ such that

$\displaystyle(1-\lambda)\big{(}(1-p)(1-\lambda)+p/|{\bm{r}}|\big{)}=1,$

we can fully cancel the effect of the dissipative noise. A straightforward
calculation shows that this is possible iff $\lambda\leq 1-|{\bm{r}}|$.

Figure 1: Toy model illustrating a noisy single qubit circuit with the input
state $\rho_{0}=\frac{1}{2}(\mathbb{I}+|{\bm{r}}|\cdot Z)$, and the output
state
$\tilde{\rho}_{1}=\mathcal{N}\circ\mathcal{R}\circ\mathcal{N}\circ\mathcal{U}(\rho_{0})$.
The red arrows in the Bloch balls (a-e) represent the mixed state after each
operation. From left to right: (a) The input state $\rho_{0}$. (b) Unitary
channel $\mathcal{U}$ rotates the input state to the ideal output state
$\rho_{1}$. (c) The depolarizing noise channel $\mathcal{N}$ “shrinks” the
vector ${\bm{r}}$. (d) The dissipative $\mathcal{R}$ gate “expands” the vector
${\bm{r}}$ beyond what is required to mitigate only the noise channel in (c).
(e) The action of the noise channel $\mathcal{N}$ applied after the
$\mathcal{R}$ gate reverts the state back to the ideal output of (b), such
that in the best case $\tilde{\rho}_{1}=\rho_{1}$.

The above conclusion is further illustrated in Fig. 1, and can be understood
intuitively as follows. The application of the noise increases the entropy of
the state, or, equivalently, shrinks its Bloch radius $|{\bm{r}}|$. Therefore,
to counter it, we first inflate the Bloch radius using $\mathcal{R}$ so that
the subsequent decrease by $\mathcal{N}$ will take us back to the ideal
output. However, if the initial state is close to being pure (i.e.,
$|{\bm{r}}|$ is close to $1$), we cannot increase its radius by much, and we
cannot fully cancel the shrinkage due to $\mathcal{N}$. Therefore, for
$\mathcal{R}$ to cancel the noise, it must be that the distance of
$|{\bm{r}}|$ from $1$ is comparable to $\lambda$.

On a high level, the above analysis shows that the dissipative variational
circuit behaves like the probabilistic error mitigation (PEC) [100, 32, 101],
as it effectively inverts the noise channel. Unlike, PEC, however, it does
this automatically without a prior characterization of the noise model. On the
other hand, it can only apply inverse noise channels that are physical, and as
such, it cannot fully correct high noise levels.

We conclude by noting that another potential advantage of using the
dissipative $\mathcal{R}$ gate in variational quantum circuits is that it can
be seen as a stochastic, parameterized version of a mid-circuit measurement.
As shown in Ref. [39], adding dissipation through random measurements can
significantly increase the variance of gradients with respect to the circuit’s
parameters. This may help the classical optimization routine in VQA algorithms
to avoid local minima and prevent rapidly decaying gradients, which may lead
to barren plateaus [102, 79].

###  III.2 The dissipative circuit ansatz

Having defined the dissipative $\mathcal{R}$ gate, our dissipative circuit
anstaz follows a simple brick-wall architecture that is common in VQA
circuits.

Figure 2: Dissipative circuit model for system size $n=4$ qubits arranged on a
ring. The input state is initialized to $|0\rangle\langle 0|^{\otimes n}$. The
dissipative parametrized layer is repeated $D$ times, and followed by a single
unitary (coherent) layer. The unitary gates $U(\theta)$ are general two-qubit
gates. The parametrized probabilistic $\operatorname{RESET}$ gates are
represented as non-unital channels $\mathcal{R}$ given in Eq. (2). The
parameters $\theta,p,\phi$ are different from gate to gate.

Alongside the dissipative $\mathcal{R}$ gates, we use general two-qubit gates
$U$, parameterized by a vector of rotation angles ${\bm{\theta}}$, applied on
neighboring qubits. For a general unitary $U\in SU(4)$, there exists a
standard KAK decomposition into three $\operatorname{CX}$ gates and $15$
elementary single qubit rotations between them [103]. We used this
decomposition and let each two-qubit unitary to have its own set of $15$
rotation angles.

Using the dissipative $\mathcal{R}$ gates and general $SU(4)$ gates, our
ansatz uses a brick-wall structure, described in Fig. 2. While we focus on the
periodic 1D case, our results can be readily generalized to other types of
lattices. From a high-level perspective, the circuit consists of two parts.
The first part is the dissipative part, which consists of $D$ even-odd brick-
wall layers, where each layer is made of entangling $SU(4)$ gates and
dissipative gates (see Fig. 2). The second part is purely coherent, consisting
of a single even-odd brick-wall layer of $SU(4)$ gates. For a system size $n$,
with $D$ dissipative layers, our numerical experiments were performed on
variational circuits with $(18D+15)n$ parameters.

The division of the circuit into two parts, loosely follows the approach taken
in Refs. [86, 87, 91, 92]. There, the first part prepares a mixed state of
orthogonal states (usually in the computational basis) with Boltzmann weights
that approximate the weights of the target Gibbs state. This mixed state is
then rotated to the actual Gibbs state using several layers of unitary gates.
Similarly, the dissipative part of our circuit can be viewed as preparing a
mixed state with the target Boltzmann weights, which is then rotated to the
target Gibbs state by the coherent part. We note, however, that in our
approach the dissipative part also contains many unitary rotations, which help
align the eigenbasis of the initial mixed state with the target Gibbs state.
This allows the coherent part to be relatively shallow at the expense of
increasing the depth of the dissipative part, which, in view of Sec. III.1, is
important for achieving dissipative noise resilience.

We conclude this section by noting that our anstaz can be used on current day
quantum computers. Indeed, our circuit is equivalent to an ensemble of circuit
instances, each instance representing a different branch of the various
$\mathcal{R}_{i}(p,{\bm{\phi}})$ probabilistic gates — either the identity
operator or the reset operator
${|{{\bm{\phi}}}\rangle\\!\,\langle{{\bm{\phi}}}|}\cdot\operatorname{Tr}_{i}(\rho)$.
The latter can be implemented by a standard RESET to
${|{0}\rangle\\!\,\langle{0}|}$ gate, followed by a local rotation
${|{0}\rangle}\to{|{{\bm{\phi}}}\rangle}$. Therefore, to run the circuit, we
first sample a sufficient number of instances, and then execute each of these
instances on a quantum computer. The final output state generated by the
ansatz is obtained by averaging the outputs of the sampled circuit instances.

##  IV Numerical Simulations

To evaluate the effectiveness of our dissipative variational circuit in
preparing Gibbs states, we conducted several numerical simulations of it,
together with simulations of other coherent anstazes. These simulations were
performed both under ideal, noiseless conditions and using a simple noise
model typically assumed for gate-based quantum computers like those offered by
IBMQ [104]. In this section we describe the systems we simulated and their
noise model, together with technical details on the optimization algorithm,
and finally the numerical results.

###  IV.1 Simulated systems

We simulated system sizes between $n=2$ to $n=6$ qubits arranged on a ring,
with the initial state $\rho_{0}=|0\rangle\langle 0|^{\otimes n}$. In all our
simulations, we considered the Gibbs states of $2$-local translation invariant
Hamiltonians defined on a ring of $n$ qubits as:

$\displaystyle H=\sum_{i=1}^{n}h_{i,i+1},$ (4)

where $h_{i,i+1}$ is the same operator, with ${\|h_{i,i+1}\|}\leq 1$, acting
on qubits $i,i+1$, with the $n+1$ qubit identified as qubit $1$. We note that
while our target Gibbs Hamiltonians were translationally invaraint, we did not
impose translation invariance on the underlying parameters of the circuit, and
let the optimizer pick the optimal parameters independently.

###  IV.2 Noise model

For the noisy simulations of our ansatz, we assumed a local Markovian noise
model without crosstalks. Since we have two distinct types of gates, there are
two types of noises to consider: the noise in the two-qubit unitary gates $U$
(using the KAK formula, which is a combination of $\operatorname{CX}$ and
rotation gates [103]) and the noise in the non-unital single-qubit channel
$\mathcal{R}$.

For the two-qubit unitary gates, we used a simplified noise model similar to
those in quantum noise simulators like IBM’s Qiskit [104, 105]. In this model,
each qubit $j$ has its own noise channel $\mathcal{N}_{j}$, so the global
noise channel on $n$ qubits factors into a tensor product of single-qubit
noise channels $\bigotimes^{n}_{j=1}\mathcal{N}^{(j)}$. Each
$\mathcal{N}^{(j)}$ combines two noise sources: dephasing (phase damping) and
amplitude damping [106]. The dephasing noise model with parameter
$\lambda_{j}$ acting on a qubit $j$ is defined as a Pauli noise channel:

$\displaystyle\mathcal{N}_{dep}^{(j)}[\rho]\stackrel{{\scriptstyle\mathrm{def}}}{{=}}(1-\lambda_{j})\rho+\lambda_{j}Z\rho
Z.$ (5)

The amplitude damping noise model with parameter $\omega_{j}$ acting on qubit
$j$ is defined by:

$\displaystyle\mathcal{N}^{(j)}_{amp}[\rho]\stackrel{{\scriptstyle\mathrm{def}}}{{=}}K_{0}(\omega_{j})\rho
K^{\dagger}_{0}(\omega_{j})+K_{1}(\omega_{j})\rho_{j}K^{\dagger}_{1}(\omega_{j}),$
(6)

where

$\displaystyle K_{0}(\omega_{j})=\begin{pmatrix}1&0\\\
0&\sqrt{1-\omega_{j}}\end{pmatrix}\,\,\,K_{1}(\omega_{j})=\begin{pmatrix}0&\sqrt{\omega_{j}}\\\
0&0\end{pmatrix}.$ (7)

The dephasing parameters $\\{\lambda_{j}\\}_{j=1}^{n}$ and amplitude damping
parameters $\\{\omega_{j}\\}_{j=1}^{n}$ were chosen to correspond to a low
noise regime with values $\sim 10^{-3}$ and were randomly generated before
each optimization run: each $\lambda_{j}$ and $\omega_{j}$ were picked
uniformly at random inside the range $[1,2]\times 10^{-3}$.

Given that the $\operatorname{CX}$ gates, which have the longest duration and
highest error rates, dominate the noise in two-qubit unitaries [107, 108], we
assumed that most noise originates from these gates. Consequently, each
application of a $\operatorname{CX}$ gate in the general two-qubit unitaries
$U$ on qubits $j$ and $j+1$ was followed by the noise channel
$\mathcal{N}_{j}\otimes\mathcal{N}_{j+1}$, and each 1-qubit rotation was
assumed ideal.

As the non-unital channel $\mathcal{R}$ from Eq. (2) is accomplished by a
regular RESET gate

$\displaystyle\mathcal{R}_{i}(p,{\bm{\phi}})[\rho]=(1-p)\rho+pU({\bm{\phi}}){|{0}\rangle\\!\,\langle{0}|}U^{\dagger}({\bm{\phi}})\cdot\operatorname{Tr}_{i}(\rho),$

where $U({\bm{\phi}})$ is a single-qubit rotation taking
${|{0}\rangle}\to{|{{\bm{\phi}}}\rangle}$, we modeled its noise by a noise
model of the standard RESET gate. Following Ref. [109], we assumed a
phenomenological noise model where the reinitialization of the qubit to the
${|{0}\rangle\\!\,\langle{0}|}$ state might fail. Specifically, we
upperbounded the activation probability $p$ by a global value $p^{*}=0.99$,
i.e., $p\leq p^{*}<1$, so that a prefect RESET is never possible.

###  IV.3 Simulation method

Our simulations were done by evolving the full density matrix according to the
dissipative variational circuit from Fig. 2, starting from the initial state
$\rho_{0}={|{0}\rangle\\!\,\langle{0}|}^{\otimes n}$. This produced an output
variational state $\rho({\bm{\theta}},{\bm{p}})$, where ${\bm{\theta}}$
represents the unitary rotation angles for both the $\mathcal{R}$ gates and
the two-qubit unitary gates, and ${\bm{p}}$ are the activation probabilities
of the $\mathcal{R}$ gates. We then optimized the variational parameters to
minimize the infidelity cost function

$\displaystyle\Phi({\bm{\theta}},{\bm{p}})\stackrel{{\scriptstyle\mathrm{def}}}{{=}}1-F\big{(}\rho({\bm{\theta}},{\bm{p}}),\rho_{G}\big{)},$
(8)

where

$\displaystyle
F\big{(}\rho({\bm{\theta}},{\bm{p}}),\rho_{G}\big{)}\stackrel{{\scriptstyle\mathrm{def}}}{{=}}\left[\operatorname{Tr}\sqrt{\sqrt{\rho({\bm{\theta}},{\bm{p}})}\rho_{G}\sqrt{\rho({\bm{\theta}},{\bm{p}})}}\right]^{2}$
(9)

is the fidelity between the variational state $\rho({\bm{\theta}},{\bm{p}})$
and the target Gibbs state $\rho_{G}$ [110].

As said by the end of Sec. II, the main focus of this work is to check the
_expressibility_ of our dissipative circuit and its resilience to noise.
Consequently, we allowed ourselves to choose the infidelity as the loss
function, even though it is hard to estimate it efficiently on a quantum
computer. Similarly, other measures of distinguishability, such as trace
distance or relative entropy between the target and variational Gibbs states,
can be equally employed to define the loss function [106].

After each simulation, gradients with respect to the variational parameters
were computed using PyTorch’s automatic differentiation engine [111]. The
parameters were then updated using the Adam [112] stochastic optimization
algorithm based on the infidelity cost function
$\Phi({\bm{\theta}},{\bm{p}})$.

The initial values of the variational parameters ${\bm{\theta}}$ were randomly
drawn at uniform in range $[-\pi,\pi]$, and the initial activation
probabilities ${\bm{p}}$ were drawn at uniform in range $[0,1]$.

Heuristically, we fixed the maximum number of optimization steps to be 2000.
We added another heuristic termination criterion for the optimization process,
wherein the optimization ceased if at step $i$, the value of the loss function
$\Phi<10^{-3}$.

##  V Results

In this section we present the results of our numerical simulations for Gibbs
state preparation using our dissipative ansatz for the ideal and noisy cases.

###  V.1 Representation capabilities of the ansatz

We first investigated the expressibility of our ansatz, by checking how well
it prepares the Gibbs state of a generic $2$-local translation-invariant
Hamiltonian, as defined by Eq. (4), across a broad range of inverse
temperatures $\beta$. To this end, we fixed the system size to $n=4$ qubits
and set the number of repetitions of the dissipative parametrized layer to
$D=n$ (see Fig. 2).

To evaluate the expressibility of our ansatz, we sampled $N=1000$ $2$-local
translation-invariant Hamiltonians defined by Eq. (4) across five inverse
temperatures: $\beta=0.5,1,2,4,8$. This range allowed us to examine both high
and low temperature regimes. For each $\beta$ value, 200 $2$-local
translation-invariant Hamiltonians were sampled.

For each sample, we ran the optimization procedure, described in Sec. IV.3, to
minimize the cost function $\Phi({\bm{\theta}},{\bm{p}})$ five times for the
noiseless case and five times for the noisy case, using the noise model
described in Sec. IV.2. Each run was initialized with different variational
parameters $\\{{\bm{\theta}},{\bm{p}}\\}$. We then selected the best result
among these five runs for both the noisy and noiseless cases, respectively.

Figure 3: The median fidelity between the estimated and target Gibbs states
for 2-local translation-invariant Hamiltonians is shown for both noiseless
(blue) and noisy (green) cases. Each data point represents the median fidelity
of 200 random Hamiltonians. For each generated Hamiltonian, five optimization
runs were performed with the best result taken for each case. The shaded area
indicates the uncertainty, represented as one standard deviation.

In Fig. 3 we present the median fidelity between the state estimated by our
ansatz and the target Gibbs state as a function of the inverse temperature
$\beta$. For each $\beta$ value the median was calculated over 200 sampled
Hamiltonians.

For the noiseless case, the most notable drop in the accuracy occurs at near
room temperature corresponding to $\beta\sim 1$. As we discuss in the Sec.
V.2, when considering specific Hamiltonians, such dip appears also in other
ansatzes [93, 113, 114] and may represent a fundamental limit on the
representation power of the variational ansatzes. Nevertheless, for the ideal
case, the median fidelity in Fig. 3 remains above $99\%$ with uncertainty of
at most $\sim 0.5\%$.

For the noisy case, the median fidelity gradually decreases towards $96\%$ as
$\beta$ increases. Such behavior is expected, as for higher $\beta$ values
representing lower temperatures, the Gibbs state becomes more and more pure,
and in light of Sec. III.1, our dissipative circuit becomes less resilient to
noise.

###  V.2 Ising and XY models

In this part we focused on two specific, well-studied models: the transverse
field Ising (TFI) and the XY models. As in previous sections, we considered
Gibbs states of translation invariant Hamiltonians arranged on a ring of $n$
qubits.

The Hamiltonian of the transverse field Ising is given by:

$\displaystyle H_{TFI}=-\sum_{j=1}^{n}X_{j}X_{j+1}-h\sum_{j=1}^{n}Z_{j},$ (10)

where $h$ represents the strength of the external field. In all our
simulations we consider three values of $h=0.5,1,1.5$. The Hamiltonian of the
XY model is:

$\displaystyle
H_{XY}=-\sum_{j=1}^{n}\left[\frac{1+\gamma}{2}X_{j}X_{j+1}+\frac{1-\gamma}{2}Y_{j}Y_{j+1}\right]-h\sum_{j=1}^{n}Z_{j},$
(11)

where $\gamma$ represents the degree of anisotropy with respect to the
$xy$-plane [115]. In our simulations, in this section, we fixed $h=0.5$ and
considered $\gamma=0.1,0.5,0.9$.

In all numerical trials that follow, for each data point, the optimization
procedure was ran ten times with different initial parameter values, and the
best result was selected.

####  V.2.1 Ising and XY models: system size dependence

In Fig. 4, we fixed the number of repetitions of the dissipative parametrized
layer to $D=n$ (see Fig. 2), and show how the precision of our ansatz varies
for system sizes $n=2,4,6$, for the target Gibbs states of TFI and XY
Hamiltonians for several inverse temperature values and for ideal and noisy
cases.

Figure 4: Fidelity between the exact Gibbs state and the state produced by our
ansatz, versus the inverse temperature $\beta$, for several system sizes
$n=2,4,6$ and fixed dissipative layer depth $D=n$, for each system size. (a,
c) Transverse field Ising model for ideal and noisy cases, respectively. (b,
d) XY model for ideal and noise cases, respectively. For each point, ten
optimization runs were performed with the best result taken.

In the noiseless case, Fig. 4 shows that our ansatz is capable of preparing
the Gibbs state of the target Hamiltonian with fidelity $>99\%$ for all values
of $\beta$ considered, for up to six qubits. Our noiseless ansatz results
closely resemble the numerical results in Refs. [113, 114] for the TFI and XY
Hamiltonians, respectively. The most notable drop in the fidelity for
different system sizes occurs at intermediate temperatures where $\beta\sim
1$. Such dip in accuracy was also observed in Refs. [113, 114], where the
authors hinted that this can represent the limit of the ansatz’s
representation power. This limitation arises from the higher number of
eigenstates needed to accurately represent the desired Gibbs state at
intermediate temperatures.

For the noisy simulations, the noise parameters were set as described in Sec.
IV.2. Fig. 4 illustrates that the fidelity between the estimated and target
Gibbs states gradually decreases with increasing values of $\beta$,
corresponding to lower temperatures, similar to the trend observed in Fig. 3.
As system size increases, the accuracy of the noisy ansatz significantly
drops, particularly at lower temperatures. This is expected in light of Sec.
III.1, which demonstrates that the noise resilience of our dissipative ansatz
drops as the target state becomes more and more pure. As the system size
increases, these errors naturally accumulate.

####  V.2.2 Ising and XY models: circuit depth

In Fig. 5, we fixed the system size at $n=6$ and examined how the depth $D$ of
the dissipative parametrized layer affects the accuracy of our ansatz in the
noiseless and noisy cases, for both the TFI Hamiltonian with $h=1$ and the XY
Hamiltonian with $\gamma=0.5$.

For the noiseless case, shown in Fig. 5(a,b), the accuracy of our ansatz
significantly increases as the depth $D$ approaches the system size $n$. For
the TFI Hamiltonian, a fidelity of over $99\%$ is achieved at $D=4$, whereas
for the XY Hamiltonian, this fidelity is reached only at $D=6$. This
difference is due to the higher entanglement in the Gibbs state of the XY
model, requiring deeper circuits to generate the necessary entanglement. For
$D=8$, the fidelity only slightly increases for both models.

For the noisy case, shown in Fig. 5(c,d), the best accuracy is reached at
$D=n$. However, unlike the noiseless case, where increasing the circuit depth
from $D=6$ to $D=8$ resulted in a slight increase in accuracy, in the noisy
case, increasing the depth from $D=6$ to $D=8$ may lead to a noticeable
decrease in accuracy, as particularly evident from Fig. 5(d). This behavior is
expected in the noisy case, as adding more layers may lead to an accumulation
of errors that cannot be mitigated by the dissipative circuit, which might
lead to an overall degradation in the fidelity.

Figure 5: Fidelity between the exact Gibbs state and the state produced by our
ansatz, versus the inverse temperature $\beta$, for different depth values of
the dissipative layer $D=2,4,6,8$ and for fixed system size $n=6$. (a,c)
Transverse field Ising model for ideal and noisy cases with $h=1$,
respectively. (b,d) XY model for ideal and noisy cases with $\gamma=0.5$,
respectively. For each point, ten optimization runs were performed with the
best result taken.

####  V.2.3 Ising and XY models: comparison with unitary ansatzes

In Fig. 6, we compare the accuracy of our method with the full unitary ansatz
proposed in Ref. [89]. The comparison uses the XY model Hamiltonian defined in
Eq. (11) on $n=6$ qubits with anisotropic parameter $\gamma=0.5$ and $h=0.5$.

For our circuits, we set $D=6$, and for the ansatz from Ref. [89], we also
chose a depth parameter $d=6$, as shown in Fig. 7 of Ref. [89]. The ansatz
from Ref. [89] uses a single ancilla qubit, requiring a total of 7 qubits for
the Gibbs state preparation on $n=6$ qubits.

We tested both ansatzes under noiseless and noisy conditions, assuming a low
noise level of $\sim 10^{-3}$ for the noisy case. Both ansatzes were tested at
inverse temperatures $\beta=2,3,4,5$ because the ansatz from Ref. [89]
underperformed at lower $\beta$ values. This underperformance can be
attributed to the fact that the specific ansatz from Fig. 7 in Ref. [89] was
tested on a smaller system. Fig. 6 clearly demonstrates the advantage of using
dissipative variational circuits over unitary circuits in the noisy case, due
to their higher noise resilience.

Another comparison can be made between our results for the noisy transverse
field Ising (TFI) model presented in Fig. 4(c) and the noisy simulation
results presented in Ref. [91], for the case of the same TFI model. Although
the noise model used for the TFI model simulations in Ref. [91] is not
identical to the one we used, it is a simplified version derived from the IBMQ
device backend. As demonstrated in Ref. [38], this noise model primarily
consists of dephasing and amplitude damping processes. Therefore, although a
direct comparison between the noisy TFI results in Ref. [91] and our noisy TFI
model results in Fig. 4(c) is not possible, as it was in Fig. 6, we can still
conclude that, generally, unitary VQAs do not perform well under noisy
conditions.

Figure 6: Fidelity comparison between our ansatz and the ansatz presented in
Fig. 7 in Ref. [89], versus the inverse temperature $\beta$, for the fixed
system size $n=6$. (a) In the noiseless case, the ansatzes perform similarly
at low temperatures. (b) Results for the noisy case, assuming low noise
regime, the accuracy of our method gradually decreases for lower temperatures,
while the ansatz presented in Fig. 7 in Ref. [89] reaches $\leq 80\%$
fidelity. For each point, ten optimization runs were performed with the best
result taken.

##  VI Discussion and outlook

In this article, we introduced the dissipative variational quantum algorithm
(D-VQA) framework, which extends the capabilities of unitary variational
quantum algorithms (VQAs) by incorporating a dissipative operation as an
intrinsic component of the variational circuit. Specifically, we presented the
dissipative $\mathcal{R}$ gate, which re-initializes the single qubit state to
a variational pure state with variational probability $p$. Through analytical
analysis using a single-qubit toy model and various numerical examples, we
demonstrated that incorporating such dissipative operations enhances the
resilience of the variational circuit against both coherent and non-coherent
errors. The resilience to coherent errors arises naturally from the unitary
part of the variational circuit [24, 22, 25, 26, 27], while the non-coherent
noise resilience is directly attributed to the dissipative component.

Our numerical examples showed that a D-VQA-based ansatz can accurately prepare
Gibbs states for a broad range of quantum many-body Hamiltonians across a wide
spectrum of temperatures. We specifically focused on well-known quantum many-
body models, such as the transverse field Ising model and the XY model. In the
ideal noiseless case, our proposed D-VQA-based ansatz achieved at least the
same degree of accuracy as the latest existing VQA-based ansatzes [89, 91, 92,
93]. In the noisy case, within the low-noise regime, our ansatz significantly
outperformed the state-of-the-art unitary VQA-based ansatz from Ref. [89].

Although the low-noise regime is not yet available on publicly accessible
quantum platforms, significant research efforts over the past years have
achieved quantum gate fidelities exceeding $99.9\%$, suggesting that this
noise regime will soon be attainable on NISQ devices [116, 107, 117, 118].
Moreover, our framework can be combined with either error mitigation or error
correction techniques — and will require a lesser reduction in the noise
levels to achieve a certain fidelity, in comparison with unitary VQAs.

Additionally, using D-VQA-based ansatzes, where the dissipative operation
involves parametrized mid-circuit measurements, may help classical
optimization routines to avoid barren plateaus [102, 79] by significantly
increasing the variance of the gradients with respect to variational circuit
parameters, as suggested by Ref. [39]. In this context, it would be
interesting to investigate the extent to which D-VQA ansatzes exhibit
resilience to the different types of barren plateaus, specifically the noise-
induced barren plateaus identified in Ref. [119].

Our work leaves several open questions. It would be valuable to gain further
theoretical understanding of D-VQAs by exploring more complex scenarios beyond
the single-qubit toy model presented here. Specifically, understanding how
entanglement between qubits influences the D-VQA circuit’s resilience to non-
coherent errors would be of great interest.

It would also be interesting to implement our D-VQA ansatzes on real quantum
hardware, potentially demonstrating high-fidelity Gibbs state preparation on a
quantum device for systems larger than just a few qubits. Such demonstrations,
coupled with further analytical investigations of D-VQAs, could provide more
evidence for the importance of dissipative operations for noise resilience,
and potentially influence future error mitigation and correction approaches.

Another direction of both theoretical and practical significance could be
comparing our work with the probabilistic error correction (PEC) framework
introduced in Refs. [100, 32, 101]. It would be interesting to explore how our
D-VQA circuits can be integrated with the PEC framework to further enhance the
noise resilience of the resulting variational circuits. Although no mitigation
method, including PEC, is perfectly scalable, integrating the PEC with the
D-VQA framework might reduce the exponential cost associated with error
mitigation.

## Acknowledgements

We thank Amit Kam, Rotem Elimelech, Netanel Lindner and Yoav Sagi for
enlightening discussions.

This research project was supported by the National Research Foundation,
Singapore and A*STAR under its CQT Bridging Grant. It was also partially
supported by the Helen Diller Quantum Center at the Technion – Israel
Institute of Technology.

## References

  * Peruzzo  _et al._ [2014] A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’Brien, A variational eigenvalue solver on a photonic quantum processor, Nature Communications 5, 10.1038/ncomms5213 (2014).
  * Farhi  _et al._ [2014] E. Farhi, J. Goldstone, and S. Gutmann, A quantum approximate optimization algorithm (2014).
  * Wecker  _et al._ [2015] D. Wecker, M. B. Hastings, and M. Troyer, Progress towards practical quantum variational algorithms, Phys. Rev. A 92, 042303 (2015).
  * Endo  _et al._ [2020] S. Endo, J. Sun, Y. Li, S. C. Benjamin, and X. Yuan, Variational quantum simulation of general processes, Phys. Rev. Lett. 125, 010501 (2020).
  * Bittel and Kliesch [2021] L. Bittel and M. Kliesch, Training variational quantum algorithms is np-hard, Phys. Rev. Lett. 127, 120502 (2021).
  * Benedetti  _et al._ [2021] M. Benedetti, M. Fiorentini, and M. Lubasch, Hardware-efficient variational quantum algorithms for time evolution, Phys. Rev. Res. 3, 033083 (2021).
  * Bonet-Monroig  _et al._ [2023] X. Bonet-Monroig, H. Wang, D. Vermetten, B. Senjean, C. Moussa, T. Bäck, V. Dunjko, and T. E. O’Brien, Performance comparison of optimization methods on variational quantum algorithms, Phys. Rev. A 107, 032407 (2023).
  * Yuan  _et al._ [2019] X. Yuan, S. Endo, Q. Zhao, Y. Li, and S. C. Benjamin, Theory of variational quantum simulation, Quantum 3, 191 (2019).
  * Hastings  _et al._ [2014] M. B. Hastings, D. Wecker, B. Bauer, and M. Troyer, Improving quantum algorithms for quantum chemistry (2014).
  * Kandala  _et al._ [2017] A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets, Nature 549, 242–246 (2017).
  * Li and Benjamin [2017] Y. Li and S. C. Benjamin, Efficient variational quantum simulator incorporating active error minimization, Phys. Rev. X 7, 021050 (2017).
  * Hempel  _et al._ [2018] C. Hempel, C. Maier, J. Romero, J. McClean, T. Monz, H. Shen, P. Jurcevic, B. P. Lanyon, P. Love, R. Babbush, A. Aspuru-Guzik, R. Blatt, and C. F. Roos, Quantum chemistry calculations on a trapped-ion quantum simulator, Phys. Rev. X 8, 031022 (2018).
  * Romero  _et al._ [2018] J. Romero, R. Babbush, J. R. McClean, C. Hempel, P. J. Love, and A. Aspuru-Guzik, Strategies for quantum computing molecular energies using the unitary coupled cluster ansatz, Quantum Science and Technology 4, 014008 (2018).
  * Xu  _et al._ [2021] X. Xu, J. Sun, S. Endo, Y. Li, S. C. Benjamin, and X. Yuan, Variational algorithms for linear algebra, Science Bulletin 66, 2181–2188 (2021).
  * Bravo-Prieto  _et al._ [2023] C. Bravo-Prieto, R. LaRose, M. Cerezo, Y. Subasi, L. Cincio, and P. J. Coles, Variational quantum linear solver, Quantum 7, 1188 (2023).
  * Huang  _et al._ [2019] H.-Y. Huang, K. Bharti, and P. Rebentrost, Near-term quantum algorithms for linear systems of equations (2019).
  * Nakanishi  _et al._ [2019] K. M. Nakanishi, K. Mitarai, and K. Fujii, Subspace-search variational quantum eigensolver for excited states, Phys. Rev. Res. 1, 033062 (2019).
  * Cerezo  _et al._ [2022] M. Cerezo, K. Sharma, A. Arrasmith, and P. J. Coles, Variational quantum state eigensolver, npj Quantum Information 8, 10.1038/s41534-022-00611-6 (2022).
  * Wang  _et al._ [2021a] X. Wang, Z. Song, and Y. Wang, Variational quantum singular value decomposition, Quantum 5, 483 (2021a).
  * Chen  _et al._ [2021] R. Chen, Z. Song, X. Zhao, and X. Wang, Variational quantum algorithms for trace distance and fidelity estimation, Quantum Science and Technology 7, 015019 (2021).
  * McArdle  _et al._ [2019] S. McArdle, T. Jones, S. Endo, Y. Li, S. C. Benjamin, and X. Yuan, Variational ansatz-based quantum simulation of imaginary time evolution, npj Quantum Information 5, 10.1038/s41534-019-0187-2 (2019).
  * McClean  _et al._ [2016] J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik, The theory of variational hybrid quantum-classical algorithms, New Journal of Physics 18, 023023 (2016).
  * Grimsley  _et al._ [2019] H. R. Grimsley, S. E. Economou, E. Barnes, and N. J. Mayhall, An adaptive variational algorithm for exact molecular simulations on a quantum computer, Nature Communications 10, 10.1038/s41467-019-10988-2 (2019).
  * Khodjasteh and Viola [2009] K. Khodjasteh and L. Viola, Dynamically error-corrected gates for universal quantum computation, Phys. Rev. Lett. 102, 080501 (2009).
  * Sharma  _et al._ [2020] K. Sharma, S. Khatri, M. Cerezo, and P. J. Coles, Noise resilience of variational quantum compiling, New Journal of Physics 22, 043006 (2020).
  * Shtanko and Movassagh [2021] O. Shtanko and R. Movassagh, Preparing thermal states on noiseless and noisy programmable quantum processors (2021).
  * Berberich  _et al._ [2024] J. Berberich, D. Fink, and C. Holm, Robustness of quantum algorithms against coherent control errors, Phys. Rev. A 109, 012417 (2024).
  * Lee  _et al._ [2023] S. Lee, J. Lee, H. Zhai, Y. Tong, A. M. Dalzell, A. Kumar, P. Helms, J. Gray, Z.-H. Cui, W. Liu, M. Kastoryano, R. Babbush, J. Preskill, D. R. Reichman, E. T. Campbell, E. F. Valeev, L. Lin, and G. K.-L. Chan, Evaluating the evidence for exponential quantum advantage in ground-state quantum chemistry, Nature Communications 14, 1952 (2023).
  * Larocca  _et al._ [2024] M. Larocca, S. Thanasilp, S. Wang, K. Sharma, J. Biamonte, P. J. Coles, L. Cincio, J. R. McClean, Z. Holmes, and M. Cerezo, A review of barren plateaus in variational quantum computing (2024).
  * Anschuetz and Kiani [2022] E. R. Anschuetz and B. T. Kiani, Quantum variational algorithms are swamped with traps, Nature Communications 13, 10.1038/s41467-022-35364-5 (2022).
  * Uvarov and Biamonte [2021] A. V. Uvarov and J. D. Biamonte, On barren plateaus and cost function locality in variational quantum algorithms, Journal of Physics A: Mathematical and Theoretical 54, 245301 (2021).
  * van den Berg  _et al._ [2023] E. van den Berg, Z. K. Minev, A. Kandala, and K. Temme, Probabilistic error cancellation with sparse pauli–lindblad models on noisy quantum processors, Nature Physics 19, 1116–1121 (2023).
  * Seif  _et al._ [2023] A. Seif, Z.-P. Cian, S. Zhou, S. Chen, and L. Jiang, Shadow distillation: Quantum error mitigation with classical shadows for near-term quantum processors, PRX Quantum 4, 010303 (2023).
  * Strikis  _et al._ [2021] A. Strikis, D. Qin, Y. Chen, S. C. Benjamin, and Y. Li, Learning-based quantum error mitigation, PRX Quantum 2, 040330 (2021).
  * Cai [2021] Z. Cai, Quantum error mitigation using symmetry expansion, Quantum 5, 548 (2021).
  * Huggins  _et al._ [2021] W. J. Huggins, S. McArdle, T. E. O’Brien, J. Lee, N. C. Rubin, S. Boixo, K. B. Whaley, R. Babbush, and J. R. McClean, Virtual distillation for quantum error mitigation, Phys. Rev. X 11, 041036 (2021).
  * Ravi  _et al._ [2022] G. S. Ravi, K. N. Smith, P. Gokhale, A. Mari, N. Earnest, A. Javadi-Abhari, and F. T. Chong, Vaqem: A variational approach to quantum error mitigation, in  _2022 IEEE International Symposium on High-Performance Computer Architecture (HPCA)_ (IEEE, 2022).
  * Ilin and Arad [2024] Y. Ilin and I. Arad, Learning a quantum channel from its steady-state, New Journal of Physics 26, 073003 (2024).
  * Wiersema  _et al._ [2023] R. Wiersema, C. Zhou, J. F. Carrasquilla, and Y. B. Kim, Measurement-induced entanglement phase transitions in variational quantum circuits, SciPost Phys. 14, 147 (2023).
  * Botelho  _et al._ [2022] L. Botelho, A. Glos, A. Kundu, J. A. Miszczak, O. Salehi, and Z. Zimborás, Error mitigation for variational quantum algorithms through mid-circuit measurements, Phys. Rev. A 105, 022441 (2022).
  * Volya and Mishra [2024] D. Volya and P. Mishra, State preparation on quantum computers via quantum steering, IEEE Transactions on Quantum Engineering 5, 1–14 (2024).
  * Chertkov  _et al._ [2023] E. Chertkov, Z. Cheng, A. C. Potter, S. Gopalakrishnan, T. M. Gatterman, J. A. Gerber, K. Gilmore, D. Gresh, A. Hall, A. Hankin, M. Matheny, T. Mengle, D. Hayes, B. Neyenhuis, R. Stutz, and M. Foss-Feig, Characterizing a non-equilibrium phase transition on a quantum computer, Nature Physics 19, 1799–1804 (2023).
  * DeCross  _et al._ [2023] M. DeCross, E. Chertkov, M. Kohagen, and M. Foss-Feig, Qubit-reuse compilation with mid-circuit measurement and reset, Phys. Rev. X 13, 041057 (2023).
  * Hua  _et al._ [2022] F. Hua, Y. Jin, Y. Chen, S. Vittal, K. Krsulich, L. S. Bishop, J. Lapeyre, A. Javadi-Abhari, and E. Z. Zhang, Exploiting qubit reuse through mid-circuit measurement and reset (2022).
  * Brandhofer  _et al._ [2023] S. Brandhofer, I. Polian, and K. Krsulich, Optimal qubit reuse for near-term quantum computers, in  _2023 IEEE International Conference on Quantum Computing and Engineering (QCE)_ (IEEE, 2023).
  * Biamonte  _et al._ [2017] J. Biamonte, P. Wittek, N. Pancotti, P. Rebentrost, N. Wiebe, and S. Lloyd, Quantum machine learning, Nature 549, 195–202 (2017).
  * Lifshitz  _et al._ [2021] Y. Y. Lifshitz, E. Bairey, E. Arbel, G. Aleksandrowicz, H. Landa, and I. Arad, Practical quantum state tomography for gibbs states (2021).
  * Zoufal  _et al._ [2021] C. Zoufal, A. Lucchi, and S. Woerner, Variational quantum boltzmann machines, Quantum Machine Intelligence 3, 10.1007/s42484-020-00033-7 (2021).
  * Zoufal [2021] C. Zoufal, Generative quantum machine learning (2021).
  * Torlai and Melko [2020] G. Torlai and R. G. Melko, Machine-learning quantum states in the nisq era, Annual Review of Condensed Matter Physics 11, 325–344 (2020).
  * Childs  _et al._ [2018] A. M. Childs, D. Maslov, Y. Nam, N. J. Ross, and Y. Su, Toward the first quantum simulation with quantum speedup, Proceedings of the National Academy of Sciences 115, 9456–9461 (2018).
  * Vernier  _et al._ [2023] E. Vernier, B. Bertini, G. Giudici, and L. Piroli, Integrable digital quantum simulation: Generalized gibbs ensembles and trotter transitions, Phys. Rev. Lett. 130, 260401 (2023).
  * Somma  _et al._ [2008] R. D. Somma, S. Boixo, H. Barnum, and E. Knill, Quantum simulations of classical annealing processes, Phys. Rev. Lett. 101, 130504 (2008).
  * Mohseni  _et al._ [2022] N. Mohseni, P. L. McMahon, and T. Byrnes, Ising machines as hardware solvers of combinatorial optimization problems, Nature Reviews Physics 4, 363–379 (2022).
  * Wiebe  _et al._ [2014] N. Wiebe, A. Kapoor, and K. M. Svore, Quantum deep learning (2014).
  * Kieferová and Wiebe [2017] M. Kieferová and N. Wiebe, Tomography and generative training with quantum boltzmann machines, Phys. Rev. A 96, 062327 (2017).
  * Amin  _et al._ [2018] M. H. Amin, E. Andriyash, J. Rolfe, B. Kulchytskyy, and R. Melko, Quantum boltzmann machine, Phys. Rev. X 8, 021050 (2018).
  * Poulin and Wocjan [2009] D. Poulin and P. Wocjan, Sampling from the thermal quantum gibbs state and evaluating partition functions with a quantum computer, Phys. Rev. Lett. 103, 220502 (2009).
  * Bilgin and Boixo [2010] E. Bilgin and S. Boixo, Preparing thermal states of quantum systems by dimension reduction, Phys. Rev. Lett. 105, 170405 (2010).
  * Riera  _et al._ [2012] A. Riera, C. Gogolin, and J. Eisert, Thermalization in nature and on a quantum computer, Phys. Rev. Lett. 108, 080402 (2012).
  * Su and Li [2020] H.-Y. Su and Y. Li, Quantum algorithm for the simulation of open-system dynamics and thermalization, Phys. Rev. A 101, 012328 (2020).
  * Holmes  _et al._ [2022] Z. Holmes, G. Muraleedharan, R. D. Somma, Y. Subasi, and B. Şahinoğlu, Quantum algorithms from fluctuation theorems: Thermal-state preparation, Quantum 6, 825 (2022).
  * Temme  _et al._ [2011] K. Temme, T. J. Osborne, K. G. Vollbrecht, D. Poulin, and F. Verstraete, Quantum metropolis sampling, Nature 471, 87–90 (2011).
  * Yung and Aspuru-Guzik [2012] M.-H. Yung and A. Aspuru-Guzik, A quantum–quantum metropolis algorithm, Proceedings of the National Academy of Sciences 109, 754–759 (2012).
  * Metcalf  _et al._ [2020] M. Metcalf, J. E. Moussa, W. A. de Jong, and M. Sarovar, Engineered thermalization and cooling of quantum many-body systems, Phys. Rev. Res. 2, 023214 (2020).
  * Verstraete  _et al._ [2004] F. Verstraete, J. J. García-Ripoll, and J. I. Cirac, Matrix product density operators: Simulation of finite-temperature and dissipative systems, Phys. Rev. Lett. 93, 207204 (2004).
  * Motta  _et al._ [2019] M. Motta, C. Sun, A. T. K. Tan, M. J. O’Rourke, E. Ye, A. J. Minnich, F. G. S. L. Brandão, and G. K.-L. Chan, Determining eigenstates and thermal states on a quantum computer using quantum imaginary time evolution, Nature Physics 16, 205–210 (2019).
  * Chen  _et al._ [2023] C.-F. Chen, M. J. Kastoryano, F. G. S. L. Brandão, and A. Gilyén, Quantum thermal state preparation (2023).
  * Layden  _et al._ [2023] D. Layden, G. Mazzola, R. V. Mishmash, M. Motta, P. Wocjan, J.-S. Kim, and S. Sheldon, Quantum-enhanced markov chain monte carlo, Nature 619, 282–287 (2023).
  * Coopmans  _et al._ [2023] L. Coopmans, Y. Kikuchi, and M. Benedetti, Predicting gibbs-state expectation values with pure thermal shadows, PRX Quantum 4, 010305 (2023).
  * Eassa  _et al._ [2023] N. M. Eassa, M. M. Moustafa, A. Banerjee, and J. Cohn, Gibbs state sampling via cluster expansions (2023).
  * Zhang  _et al._ [2023] D. Zhang, J. L. Bosse, and T. Cubitt, Dissipative quantum gibbs sampling (2023).
  * Brandão and Kastoryano [2018] F. G. S. L. Brandão and M. J. Kastoryano, Finite correlation length implies efficient preparation of quantum thermal states, Communications in Mathematical Physics 365, 1–16 (2018).
  * Chowdhury and Somma [2016] A. N. Chowdhury and R. D. Somma, Quantum algorithms for gibbs sampling and hitting-time estimation (2016), arXiv:1603.02940 [quant-ph] .
  * Gilyén  _et al._ [2022] A. Gilyén, S. Lloyd, I. Marvian, Y. Quek, and M. M. Wilde, Quantum algorithm for petz recovery channels and pretty good measurements, Phys. Rev. Lett. 128, 220502 (2022).
  * Lu  _et al._ [2021] S. Lu, M. C. Bañuls, and J. I. Cirac, Algorithms for quantum simulation at finite energies, PRX Quantum 2, 020321 (2021).
  * Alhambra [2023] A. M. Alhambra, Quantum many-body systems in thermal equilibrium, PRX Quantum 4, 040201 (2023).
  * Tilly  _et al._ [2022] J. Tilly, H. Chen, S. Cao, D. Picozzi, K. Setia, Y. Li, E. Grant, L. Wossnig, I. Rungger, G. H. Booth, and J. Tennyson, The variational quantum eigensolver: A review of methods and best practices, Physics Reports 986, 1–128 (2022).
  * Cerezo  _et al._ [2021] M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, and P. J. Coles, Variational quantum algorithms, Nature Reviews Physics 3, 625–644 (2021).
  * Sagastizabal  _et al._ [2021] R. Sagastizabal, S. P. Premaratne, B. A. Klaver, M. A. Rol, V. Negîrneac, M. S. Moreira, X. Zou, S. Johri, N. Muthusubramanian, M. Beekman, C. Zachariadis, V. P. Ostroukh, N. Haider, A. Bruno, A. Y. Matsuura, and L. DiCarlo, Variational preparation of finite-temperature states on a quantum computer, npj Quantum Information 7, 10.1038/s41534-021-00468-1 (2021).
  * Wang  _et al._ [2023] X. Wang, X. Feng, T. Hartung, K. Jansen, and P. Stornati, Critical behavior of the ising model by preparing the thermal state on a quantum computer, Physical Review A 108, 10.1103/physreva.108.022612 (2023).
  * Patti  _et al._ [2022] T. L. Patti, O. Shehab, K. Najafi, and S. F. Yelin, Markov chain monte carlo enhanced variational quantum algorithms, Quantum Science and Technology 8, 015019 (2022).
  * Cohn  _et al._ [2020] J. Cohn, F. Yang, K. Najafi, B. Jones, and J. K. Freericks, Minimal effective gibbs ansatz: A simple protocol for extracting an accurate thermal representation for quantum simulation, Phys. Rev. A 102, 022622 (2020).
  * Verdon  _et al._ [2017] G. Verdon, M. Broughton, and J. Biamonte, A quantum algorithm to train neural networks using low-depth circuits (2017).
  * Liu  _et al._ [2021] J.-G. Liu, L. Mao, P. Zhang, and L. Wang, Solving quantum statistical mechanics with variational autoregressive networks and quantum circuits, Machine Learning: Science and Technology 2, 025011 (2021).
  * Martyn and Swingle [2019] J. Martyn and B. Swingle, Product spectrum ansatz and the simplicity of thermal states, Phys. Rev. A 100, 032107 (2019).
  * Verdon  _et al._ [2019] G. Verdon, J. Marks, S. Nanda, S. Leichenauer, and J. Hidary, Quantum hamiltonian-based models and the variational quantum thermalizer algorithm (2019).
  * Zhu  _et al._ [2020] D. Zhu, S. Johri, N. M. Linke, K. A. Landsman, C. Huerta Alderete, N. H. Nguyen, A. Y. Matsuura, T. H. Hsieh, and C. Monroe, Generation of thermofield double states and critical ground states with a quantum computer, Proceedings of the National Academy of Sciences 117, 25402–25406 (2020).
  * Wang  _et al._ [2021b] Y. Wang, G. Li, and X. Wang, Variational quantum gibbs state preparation with a truncated taylor series, Phys. Rev. Appl. 16, 054035 (2021b).
  * Foldager  _et al._ [2022] J. Foldager, A. Pesah, and L. K. Hansen, Noise-assisted variational quantum thermalization, Scientific Reports 12, 10.1038/s41598-022-07296-z (2022).
  * Consiglio  _et al._ [2023a] M. Consiglio, J. Settino, A. Giordano, C. Mastroianni, F. Plastina, S. Lorenzo, S. Maniscalco, J. Goold, and T. J. G. Apollaro, Variational gibbs state preparation on nisq devices (2023a).
  * Consiglio [2023a] M. Consiglio, Variational quantum algorithms for gibbs state preparation (2023a).
  * Selisko  _et al._ [2023] J. Selisko, M. Amsler, T. Hammerschmidt, R. Drautz, and T. Eckl, Extending the variational quantum eigensolver to finite temperatures, Quantum Science and Technology 9, 015026 (2023).
  * Reif [1965] F. Reif,  _Fundamentals of Statistical and Thermal Physics_ (McGraw-Hill, New York, 1965).
  * Wu and Hsieh [2019] J. Wu and T. H. Hsieh, Variational thermal quantum simulation via thermofield double states, Phys. Rev. Lett. 123, 220502 (2019).
  * Sewell  _et al._ [2022] T. J. Sewell, C. D. White, and B. Swingle, Thermal multi-scale entanglement renormalization ansatz for variational gibbs state preparation (2022).
  * Chowdhury  _et al._ [2020] A. N. Chowdhury, G. H. Low, and N. Wiebe, A variational quantum algorithm for preparing quantum gibbs states (2020).
  * Araz  _et al._ [2024] J. Y. Araz, R. G. Jha, F. Ringer, and B. Sambasivam, Thermal state preparation of the syk model using a variational quantum algorithm (2024).
  * Warren  _et al._ [2022] A. Warren, L. Zhu, E. Barnes, and S. Economou, Adaptive variational algorithms for quantum gibbs state preparation, Bulletin of the American Physical Society 67 (2022).
  * Temme  _et al._ [2017] K. Temme, S. Bravyi, and J. M. Gambetta, Error mitigation for short-depth quantum circuits, Phys. Rev. Lett. 119, 180509 (2017).
  * Gupta  _et al._ [2023] R. S. Gupta, E. v. d. Berg, M. Takita, D. Riste, K. Temme, and A. Kandala, Probabilistic error cancellation for dynamic quantum circuits (2023).
  * McClean  _et al._ [2018] J. R. McClean, S. Boixo, V. N. Smelyanskiy, R. Babbush, and H. Neven, Barren plateaus in quantum neural network training landscapes, Nature Communications 9, 10.1038/s41467-018-07090-4 (2018).
  * Vatan and Williams [2004] F. Vatan and C. Williams, Optimal quantum circuits for general two-qubit gates, Phys. Rev. A 69, 032315 (2004).
  * ref [2022a] Ibm quantum. https://quantum-computing.ibm.com/, 2022 (2022a).
  * ref [2022b] Qiskit: An open-source framework for quantum computing (2022b).
  * Nielsen and Chuang [2010] M. A. Nielsen and I. L. Chuang,  _Quantum computation and quantum information_ (Cambridge university press, 2010).
  * Kjaergaard  _et al._ [2020] M. Kjaergaard, M. E. Schwartz, J. Braumuller, P. Krantz, J. I.-J. Wang, S. Gustavsson, and W. D. Oliver, Superconducting qubits: Current state of play, Annual Review of Condensed Matter Physics 11, 369 (2020), https://doi.org/10.1146/annurev-conmatphys-031119-050605 .
  * Jurcevic  _et al._ [2021] P. Jurcevic, A. Javadi-Abhari, L. S. Bishop, I. Lauer, D. F. Bogorin, M. Brink, L. Capelluto, O. Günlük, T. Itoko, and N. Kanazawa, Demonstration of quantum volume 64 on a superconducting quantum computing system, Quantum Science and Technology 6, 025020 (2021).
  * Rost  _et al._ [2021] B. Rost, L. Del Re, N. Earnest, A. F. Kemper, B. Jones, and J. K. Freericks, Demonstrating robust simulation of driven-dissipative problems on near-term quantum computers, arXiv , arXiv:2108.01183 (2021).
  * Uhlmann [2010] A. Uhlmann, Transition probability (fidelity) and its relatives, Foundations of Physics 41, 288–298 (2010).
  * Paszke  _et al._ [2019] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala, Pytorch:an imperative style, high-performance deep learning library, in  _Advances in Neural Information Processing Systems 32_ (Curran Associates, Inc., 2019) pp. 8024–8035.
  * Kingma and Ba [2015] D. P. Kingma and J. L. Ba, Adam: A method for stochastic optimization, 3rd International Conference on Learning Representations, ICLR 2015 - Conference Track Proceedings , 1 (2015), arXiv:1412.6980 .
  * Consiglio  _et al._ [2023b] M. Consiglio, J. Settino, A. Giordano, C. Mastroianni, F. Plastina, S. Lorenzo, S. Maniscalco, J. Goold, and T. J. G. Apollaro, Variational Gibbs State Preparation on NISQ devices, arXiv (2023b), arXiv:2303.11276 .
  * Consiglio [2023b] M. Consiglio, Variational Quantum Algorithms for Gibbs State Preparation, arXiv (2023b), arXiv:2305.17713 .
  * Lieb  _et al._ [1961] E. Lieb, T. Schultz, and D. Mattis, Two soluble models of an antiferromagnetic chain, Annals of Physics 16, 407 (1961).
  * Zahedinejad  _et al._ [2015] E. Zahedinejad, J. Ghosh, and B. C. Sanders, High-fidelity single-shot toffoli gate via quantum control, Phys. Rev. Lett. 114, 200502 (2015).
  * Xie  _et al._ [2023] T. Xie, Z. Zhao, S. Xu, X. Kong, Z. Yang, M. Wang, Y. Wang, F. Shi, and J. Du, 99.92%-fidelity cnot gates in solids by noise filtering, Phys. Rev. Lett. 130, 030601 (2023).
  * Zhang  _et al._ [2020] E. J. Zhang, S. Srinivasan, N. Sundaresan, D. F. Bogorin, Y. Martin, J. B. Hertzberg, J. Timmerwilke, E. J. Pritchett, J.-B. Yau, C. Wang, W. Landers, E. P. Lewandowski, A. Narasgond, S. Rosenblatt, G. A. Keefe, I. Lauer, M. B. Rothwell, D. T. McClure, O. E. Dial, J. S. Orcutt, M. Brink, and J. M. Chow, High-fidelity superconducting quantum processors via laser-annealing of transmon qubits (2020).
  * Wang  _et al._ [2021c] S. Wang, E. Fontana, M. Cerezo, K. Sharma, A. Sone, L. Cincio, and P. J. Coles, Noise-induced barren plateaus in variational quantum algorithms, Nature Communications 12, 10.1038/s41467-021-27045-6 (2021c).


# N-channel parity-time symmetry

Ege Özgün<EMAIL_ADDRESS>Department of Physics Engineering,
Hacettepe University, 06800, Ankara, Türkiye

###### Abstract

We calculated the eigenvalues for a general N-channel coupled system with
parity-time symmetry due to equal loss/gain. We found that the eigenspectrum
displays a mixing of parity-time symmetric and broken phases, with $N-2$ of
the eigenvalues being parity-time broken whereas the remaining two being
either parity-time symmetric or broken depending on the loss/gain and coupling
parameters. Our results also show that mixing of parity-time symmetric and
parity-time broken phases can only be obtained for at least four-channels if
other degrees of freedom like polarization is not taken into account.

Parity-time symmetry, equal loss/gain, non-Hermitian systems

##  I Introduction

Ignited with Bender and Boettcher’s seminal paper [1], parity-time $({\cal
PT})$ symmetry became a significant theoretical and experimental area of
interest. In their work, Bender and Boettcher showed that non-Hermitian
Hamiltonians in quantum mechanics can have a partial or a full real
eigenspectrum provided that the Hamiltonian is ${\cal PT}$-symmetric. Then in
a series of papers, Mostafazadeh generalized the idea and theoretically
demonstrated that the Hamiltonians belonging to the class of pseudo-Hermitian
Hamiltonians display real spectrum and ${\cal PT}$-symmetric Hamiltonians also
belong to that class [2, 3, 4]. Following the theoretical achievements, ${\cal
PT}$ symmetry has found numerous applications, including implementations in
optical systems [5, 6, 7], waveguides [8] and single-mode lasers [9], study of
polarization depending scattering in photonic systems [10] and scattering in
spin-$1/2$ systems in quantum mechanics [11], implementations in coupled RLC
circuits [12] and optoelectronic oscillators (OEO)s [13, 14]. Recently, a
generalized four-channel ${\cal PT}$ symmetry scheme was theoretically
suggested for the equal loss/gain setting, that utilizes two different
coupling constants [15]. Here we will generalize this scheme two N-channels by
assuming the same coupling between all channels and theoretically show that
$N-2$ of the eigenvalues of the coupled system are ${\cal PT}$-broken and the
remaining two eigenvalues can be either ${\cal PT}$-symmetric or ${\cal
PT}$-broken depending on the coupling parameter and the loss/gain value. Thus,
for a wide range of parameters the N-channel case displays a mixed ${\cal PT}$
spectrum with coexisting ${\cal PT}$-symmetric and ${\cal PT}$-broken phases.

The remaining part of the manuscript is organized as the following. In Section
II we will go over the simplest scenario i.e., two-channel case as a warm up,
then in Section III we will study the four-channel case. In Section IV we will
show our results for the most general case of N-channels, which is the main
result of this manuscript. We finally conclude with Section V.

##  II Two-channel case

A system consisting of two modes that couple to each other can be expressed in
terms of two coupled differential equations, which are known as the coupled-
mode equations:

$\displaystyle i\dot{a}_{1}$ $\displaystyle=$
$\displaystyle\omega_{1}a_{1}-iga_{1}-\kappa_{12}a_{2}$ $\displaystyle
i\dot{a}_{2}$ $\displaystyle=$
$\displaystyle\omega_{2}a_{2}+ila_{2}-\kappa_{21}a_{1}$ (1)

In the above equation, $a_{i}$’s and $\omega_{i}$’s are the amplitudes and
(angular) frequencies of the respective channels $(i=1,2)$, dot denotes time
derivative, $\kappa_{12}$ and $\kappa_{21}$ are coupling constants between the
channels and $g$ and $l$ are gain and loss parameters, respectively. Figure 1
gives the schematic structure of such two-mode coupling. This formalism is
quite general, thus it can be used to describe different systems, such as
${\cal PT}$-symmetric OEOs, coupled waveguides or coupled RLC oscillators.

Figure 1: A general two-channel scheme. For the equal loss/gain configuration
$(g=l)$, ${\cal PT}$ transition can be obtained by tuning either the loss/gain
parameter or the coupling constant.

From now on, we will assume equal frequencies in all channels
$(\omega\coloneqq\omega_{1}=\omega_{2})$ and equal loss/gain parameters
$(\gamma\coloneqq l=g)$. Moreover we will assume real coupling constants,
therefore we have $(\kappa\coloneqq\kappa_{12}=\kappa_{21})$. With these
assumptions, which can be physically realizable, we will obtain the ${\cal
PT}$-symmetric form, that can display ${\cal PT}$-transition for varying
coupling constant or loss/gain parameters. The eigenspectrum of the system can
be found by solving the characteristic equation for the below matrix obtained
from Equation 1:

$\displaystyle{\cal M}_{2}=\begin{pmatrix}\omega+i\gamma&\kappa\\\
\kappa&\omega-i\gamma\end{pmatrix}$ (2)

By solving $\det[{\cal M}_{2}-\mathds{1}_{2}\lambda]=0$ (where
$\mathds{1}_{2}$ is a $2\times 2$ identity matrix), we obtain the
eigenspectrum of the two-channel system:

$\displaystyle\lambda_{1,2}=\omega\pm\sqrt{\kappa^{2}-\gamma^{2}}$ (3)

Above equation displays the canonical eigenspectrum of a non-Hermitian ${\cal
PT}$-symmetric system. ${\cal PT}$ transition from ${\cal PT}$-symmetric phase
to ${\cal PT}$-broken phase occurs for $\gamma>\kappa$, where the eigenvalues
becomes complex. A significant result worth mentioning is that it is not
possible to obtain mixing of ${\cal PT}$-symmetric and ${\cal PT}$-broken
phases for the two-channel case (unless other degrees of freedom such as
polarization is exploited [10]). We will show in the following sections that,
at least four-channels are required to obtain mixing of ${\cal PT}$-symmetric
and ${\cal PT}$-broken phases.

##  III Four-channel case

In this Section, we will re-derive the result of Reference [15] but we will
take all the coupling constants to be equal which will be useful for the
following section when deriving the general result for the N-channel case. For
the four-channel case in which channels display equal loss/gain, similar to
the two-channel system, we can write four coupled differential equations,
which can be then written in the form of a matrix whose characteristic
equation gives the four eigenvalues of the system:

$\displaystyle{\cal
M}_{4}=\begin{pmatrix}\omega+i\gamma&\kappa&\kappa&\kappa\\\
\kappa&\omega-i\gamma&\kappa&\kappa\\\ \kappa&\kappa&\omega+i\gamma&\kappa\\\
\kappa&\kappa&\kappa&\omega-i\gamma\end{pmatrix}$ (4)

Calculating $\det[{\cal M}_{4}-\mathds{1}_{4}\lambda]=0$ (where
$\mathds{1}_{4}$ is a $4\times 4$ identity matrix), yields the eigenspectrum
of the four-channel system:

$\displaystyle\lambda_{1,2}$ $\displaystyle=$ $\displaystyle\omega-\kappa\pm
i\gamma$ $\displaystyle\lambda_{3,4}$ $\displaystyle=$
$\displaystyle\kappa+\omega\pm\sqrt{4\kappa^{2}-\gamma^{2}}$ (5)

Above equations display mixing of ${\cal PT}$-symmetric and ${\cal PT}$-broken
phases as we mentioned in the previous section. The first two eigenvalues
$(\lambda_{1,2})$ are in ${\cal PT}$-broken phase without any dependence on
loss/gain and coupling parameters, whereas the last two $(\lambda_{3,4})$ can
be either ${\cal PT}$-symmetric or ${\cal PT}$-broken depending on $\kappa$
and $\gamma$. Specifically, for $4\kappa^{2}>\gamma^{2}$ ${\cal PT}$-symmetric
phase is obtained, giving rise to an overall mixing of ${\cal PT}$-symmetric
and ${\cal PT}$-broken phases. On the other hand, when
$4\kappa^{2}<\gamma^{2}$ all eigenvalues are in the ${\cal PT}$-broken phase.

##  IV N-channel case

We are now in a position to study the N-channel system. We limit our interest
in the case where $N/2$ of the channels are loss channels whereas the other
half are gain channels with equal loss/gain parameters to achieve ${\cal
PT}$-symmetry conditions. Moreover, we will assume that all channels couple to
each other with the same coupling constant. Same frequencies in all channels
are also assumed, similar to the two-channel and four-channel cases. It is
important to mention that, in order to obtain ${\cal PT}$-symmetry conditions
$N$ must be even, which should be apparent from the equal loss/gain condition.

In our theoretical calculations, again we do not make any restrictions on the
physical nature of the modes, so the theoretical model is quite general and
can be applied to coupled waveguides, coupled RLC-oscillators, OEOs or any
system displaying mode-coupling. Therefore, the results we will obtain for the
eigenspectrum is going to be valid for any N-channel system respecting the
conditions mentioned above for satisfying ${\cal PT}$-symmetry conditions.

Figure 2: Schematic depiction of N-channel coupling. Each of the N modes
coming from the left are coupled to other N-1 modes.

Figure 2 displays the schematic demonstration of the N-channel system. We will
again write the coupled-mode equations in the matrix form. For the N-channel
case we will have an $N\times N$ matrix, which will yield $N$ eigenvalues from
its characteristic equation. The form of the matrix for the N-channel case is
given below:

$\displaystyle{\cal
M}_{N}=\begin{pmatrix}\omega+i\gamma&\kappa&\ldots&\ldots&\kappa\\\
\kappa&\omega-i\gamma&\kappa&\ldots&\vdots\\\
\vdots&\kappa&\ddots&\ddots&\vdots\\\ \vdots&\ldots&\ddots&\ddots&\kappa\\\
\kappa&\ldots&\ldots&\kappa&\omega-i\gamma\par\end{pmatrix}$ (6)

The eigenvalues can again be found by solving $\det[{\cal
M}_{N}-\mathds{1}_{N}\lambda]=0$ (where $\mathds{1}_{N}$ is an $N\times N$
identity matrix) which can be calculated analytically from the equation for
$\lambda$ given below:

$\displaystyle\Big{[}$
$\displaystyle\lambda-(\omega-\kappa+i\gamma)\Big{]}^{\left(\frac{N}{2}-1\right)}\Big{[}\lambda-(\omega-\kappa-i\gamma)\Big{]}^{\left(\frac{N}{2}-1\right)}$
$\displaystyle\times$ $\displaystyle\Bigg{[}$
$\displaystyle\lambda^{2}-2\left(\omega+\left[\frac{N}{2}-1\right]\kappa\right)\lambda+\left(\omega+\left[\frac{N}{2}-1\right]\kappa\right)^{2}+\gamma^{2}-\left(\frac{N}{2}\right)^{2}\kappa^{2}\Bigg{]}=0$
(7)

The above equation is an Nth order polynomial equation for $\lambda$ with its
solutions yielding the $N$ eigenvalues of the system:

$\displaystyle\lambda_{1}=\lambda_{2}=\ldots=\lambda_{(N/2-1)}=\omega-\kappa+i\gamma$
$\displaystyle\lambda_{(N/2)}=\lambda_{(N/2+1)}=\ldots=\lambda_{(N-2)}=\omega-\kappa-i\gamma$
$\displaystyle\lambda_{(N-1)}=\left[\frac{N}{2}-1\right]\kappa+\omega+\sqrt{\left[\frac{N\kappa}{2}\right]^{2}-\gamma^{2}}$
$\displaystyle\lambda_{(N)}=\left[\frac{N}{2}-1\right]\kappa+\omega-\sqrt{\left[\frac{N\kappa}{2}\right]^{2}-\gamma^{2}}$
(8)

Figure 3: Counter plot of the function $f(N,\gamma)=[N/2]^{2}-\gamma^{2}$ for
a wide range of the parameters $N$ and $\gamma$. The red line separates the
${\cal PT}$-symmetric and ${\cal PT}$-broken regions for the eigenvalues
$\lambda_{(N-1)}$ and $\lambda_{(N)}$ where the lighter region corresponds to
${\cal PT}$-symmetric phase and the darker region corresponds to ${\cal
PT}$-broken phase.

As it can be seen from Equation 8, $N-2$ of the eigenvalues, consisting of two
pairs with $(N-2)/2$-fold degeneracy are in ${\cal PT}$-broken phase. The
remaining two eigenvalues can be either both ${\cal PT}$-symmetric or ${\cal
PT}$-broken depending on the values of loss/gain and coupling parameters. By
observing Equation 7, we can also arrive to the same conclusion which we
mentioned in Section II that for mixed ${\cal PT}$-symmetric and ${\cal
PT}$-broken phases, we need at least four-channels unless other degrees of
freedom such as polarization are taken into account.

Let us now define a function for fixed coupling constant, say $\kappa=1$ to
illustrate the ${\cal PT}$-symmetric and ${\cal PT}$-broken phases of
$\lambda_{(N-1)}$ and $\lambda_{(N)}$ for different values of $N$ and
$\gamma$:

$\displaystyle f(N,\gamma)=[N/2]^{2}-\gamma^{2}.$ (9)

The function $f(N,\gamma)$ has two distinct regions. When $f(N,\gamma)<0$ the
two eigenvalues $\lambda_{(N-1)}$ and $\lambda_{(N)}$ are in the ${\cal
PT}$-broken phase, whereas for $f(N,\gamma)>0$ they are ${\cal PT}$-symmetric,
which demonstrates an overall mixing of the eigenspectrum with coexisting
${\cal PT}$-symmetric and ${\cal PT}$-broken phases. We plot $f(N,\gamma)$ for
a wide range of the parameters $N$ and $\gamma$ in Figure 3. It can be seen
that as the number of channels N increases, ${\cal PT}$-symmetric phase
becomes more and more robust for fixed loss/gain value.

##  V Conclusion

We calculated the eigenspectrum of an N-channel system consisting of even
number of channels with $N/2$ pairs of them having equal loss/gain and showed
that the spectrum consists of two pairs of ${\cal PT}$-broken eigenvalues each
with a degeneracy of $(N-2)/2$ and two other eigenvalues that can either be in
the ${\cal PT}$-symmetric or ${\cal PT}$-broken phase depending on the values
of the loss/gain parameters and the coupling constant. Another significant
result we obtained is the lower limit for obtaining mixed ${\cal
PT}$-symmetric and ${\cal PT}$-broken phases, which is at least four-channels,
assuming no other degrees of freedom such as polarization is exploited.
Moreover, our results showed that, for the increasing number of channels, the
${\cal PT}$-symmetric phase becomes more robust for fixed loss/gain values.

## Acknowledgment

The author acknowledges fruitful discussions with Gökhan Alkaç.

## References

  * [1] C. M. Bender and S Boettcher, Real Spectra in Non-Hermitian Hamiltonians Having PT Symmetry, Physical Review Letters, 80, 5243 (1998).
  * [2] A. Mostafazadeh, Pseudo-Hermiticity versus PT-symmetry: The necessary condition for the reality of the spectrum of a non-Hermitian Hamiltonian, Journal of Mathematical Physics 43, 205 (2002).
  * [3] A. Mostafazadeh, Pseudo-Hermiticity versus PT-symmetry. II. A complete characterization of non-Hermitian Hamiltonians with a real spectrum, Journal of Mathematical Physics 43, 2814 (2002).
  * [4] A. Mostafazadeh, Pseudo-Hermiticity versus PT-symmetry III: Equivalence of pseudo-Hermiticity and the presence of antilinear symmetries, Journal of Mathematical Physics 43, 3944 (2002).
  * [5] R. El-Ganainy, K. G. Makris, D. N. Christodoulides, and Ziad H. Musslimani, Theory of coupled optical PT-symmetric structures, Optics Letters 32, 2632 (2007).
  * [6] K. G. Makris R. El-Ganainy and D. N. Christodoulides, Beam Dynamics in PT Symmetric Optical Lattices, Physical Review Letters, 100, 103904 (2008).
  * [7] Christian E. Rüter, K. G. Makris, Ramy El-Ganainy, D. N. Christodoulides, Mordechai Segev and Detlef Kip, Observation of parity–time symmetry in optics, Nature Physics, 6, 192 (2010).
  * [8] Shachar Klaiman, Uwe Günther, and Nimrod Moiseyev, Visualization of Branch Points in PT-Symmetric Waveguides, Physical Review Letters, 101, 080402 (2008).
  * [9] L. Feng, Z. J. Wong, R.-M. Ma, Y. Wang, X. Zhang, Single-mode laser by parity-time symmetry breaking, Science, 346, 972 (2014).
  * [10] E. Özgün, A. E. Serebryannikov, E. Ozbay and C. M. Soukoulis, Broadband mixing of PT-symmetric and PT-broken phases in photonic heterostructures with a one-dimensional loss/gain bilayer, Scientific Reports, 7, 15504 (2017).
  * [11] E. Özgün, T. Hakioğlu and E. Ozbay, Scattering of spin-1/2 particles from a PT-symmetric complex potential, Europhysics Letters, 131, 11001 (2020).
  * [12] Y. Choi, C. Hahn, J. W. Yoon and S. H. Song, Observation of an anti-PT-symmetric exceptional point and energy-difference conserving dynamics in electrical circuit resonators, Nature Communications, 9, 2182 (2018).
  * [13] J. Zhang and J. Yao, Parity-time–symmetric optoelectronic oscillator, Science Advances, 4, eaar678 (2018).
  * [14] E. Özgün, F. Uyar, T. Kartaloglu, E. Ozbay and I. Ozdur, A parity-time-symmetric optoelectronic oscillator with polarization multiplexed channels, Journal of Optics, 24, 055802 (2022).
  * [15] E. Özgün, E. Ozbay and I. Ozdur, Four-channel parity-time symmetry, Europhysics Letters, 140, 10001 (2022).