# Engineered dissipation to mitigate barren plateaus

Antonio Sannia  Institute for Cross-Disciplinary Physics and Complex Systems
(IFISC) UIB-CSIC, Campus Universitat Illes Balears, 07122, Palma de Mallorca,
Spain.  Francesco Tacchino  IBM Quantum, IBM Research Europe – Zurich,
Säumerstrasse 4, 8803 Rüschlikon, Switzerland  Ivano Tavernelli  IBM Quantum,
IBM Research Europe – Zurich, Säumerstrasse 4, 8803 Rüschlikon, Switzerland
Gian Luca Giorgi  Institute for Cross-Disciplinary Physics and Complex Systems
(IFISC) UIB-CSIC, Campus Universitat Illes Balears, 07122, Palma de Mallorca,
Spain.  Roberta Zambrini  Institute for Cross-Disciplinary Physics and Complex
Systems (IFISC) UIB-CSIC, Campus Universitat Illes Balears, 07122, Palma de
Mallorca, Spain.

###### Abstract

Variational quantum algorithms represent a powerful approach for solving
optimization problems on noisy quantum computers, with a broad spectrum of
potential applications ranging from chemistry to machine learning. However,
their performances in practical implementations crucially depend on the
effectiveness of quantum circuit training, which can be severely limited by
phenomena such as barren plateaus. While, in general, dissipation is
detrimental for quantum algorithms, and noise itself can actually induce
barren plateaus, here we describe how the inclusion of properly engineered
Markovian losses after each unitary quantum circuit layer can restore the
trainability of quantum models. We identify the required form of the
dissipation processes and establish that their optimization is efficient. We
benchmark our proposal in both a synthetic and a practical quantum chemistry
example, demonstrating its effectiveness and potential impact across different
domains.

##  I Introduction

While dissipation is generally detrimental to quantum technologies and is, in
fact, a limiting factor for currently available noisy quantum devices [1],
there are important exceptions. On the one hand, going beyond unitary
operations is not only needed to account for quantum measurements [2] but also
enables several applications ranging from quantum state preparation and
control to error correction [3]. On the other hand, noisy quantum platforms
are well suited for computing [4] in a hybrid quantum-classical setting, for
instance using the popular class of variational quantum algorithms (VQAs) [5].
In this framework, a parametric quantum circuit is used to prepare complex
quantum states and to evaluate, through sampling or with suitable quantum
measurements, a cost function that would be otherwise expensive to compute
classically. The optimization of the circuit parameters, aimed at minimizing
such cost, is instead assigned to a classical optimizer. As an example,
variational quantum eigensolvers (VQEs) [6, 7] can be used to approximate the
ground state of a given Hamiltonian H through a trial quantum circuit, also
called an ansatz, by minimizing the cost function given by the expectation
value of H. More in general, VQAs can be applied to classical optimization
problems [8], linear systems of equations [9, 10, 11], quantum simulations
[12, 13, 14, 15], quantum data compression [16], quantum machine learning [17,
18, 19], generative models [20, 21, 22], quantum foundations [23], quantum
compiling [24, 25, 26, 27], and quantum error correction [28]. In many of
these cases, the VQA can also be formulated as a ground-state problem with a
problem-specific (rather than a physically motivated) Hamiltonian, the choice
of which is generally non-trivial [5].

Despite their potential advantages, VQAs are known to suffer from a number of
bottlenecks that still prevent their implementation at scale. One particularly
serious drawback, specific to quantum operations, is the phenomenon known as
barren plateaus [29]. These manifest as a vanishing gradient of the loss
function with respect to the model parameters and can lead to severe
limitations in the training efficiency. Barren plateaus build up as the system
size increases, hence directly hindering the scalability of VQAs.
Specifically, it was found that when a quantum circuit ansatz reaches the
2-design limit the probability of randomly inducing a non-negligible gradient
decreases exponentially with the number of qubits [29]. This condition is
relatively easy to satisfy [30, 31, 32] implying a generalized loss of all
possible quantum speedups, even if the optimization strategy does not include
gradient computation [33].

Recent findings indicate that barren plateaus can originate from several
factors, including high circuit expressibility [34, 35], concentration of the
cost function [36], noise [37], entanglement excess [38], and globality of the
cost function [39, 40]. All these barren-plateau sources can be unified by
means of a Lie algebra theory that applies to a general (yet unitary) setting
[41, 42]. Several techniques to mitigate barren plateaus have been reported,
such as initialization strategies [43, 44], transferability of smooth
solutions [45], entanglement limitations [46], correlations and restrictions
of circuit parameters [34, 47], classical shadows [48], pre-trainings [49,
50], and layerwise learning [51, 52]. These mitigation strategies generally
consist in appropriately constraining the unitary ansatz. In this work, we
change the focus acting on the problem Hamiltonian without adapting the
quantum circuit. Remarkably, the occurrence of barren plateaus is closely
related to the unitarity of the ansatz [41, 42], while assessing these
phenomena beyond a fully unitary framework presents a promising research
avenue. Our goal is therefore to demonstrate the potential of non-unitary,
open-system dynamics as a powerful strategy to overcome the trainability
barrier and to ensure efficient convergence. To achieve it, it is reasonable
to expect that a non-trivial engineering of the dissipation processes will be
required. Indeed, while it has been already suggested that non-unitary
operations can increase the accuracy of ground-state VQE calculations in
quantum chemistry [53, 54], generic noise is known to actually induce barren
plateaus [37, 55]. In parallel, dissipation implemented simply by discarding
qubit registers, as in the class of so-called dissipative quantum neural
network models [56, 57, 58], is also not sufficient to ensure trainability, in
general [56]. The proof that a suitable set of operations can in fact be
constructed constitutes the main technical contribution of our work.

Building on previous results obtained in the field of engineered quantum
dissipation [59, 60], we consider a Markov dissipation modeled by a Gorini-
Kossakowski-Sudarshan-Lindblad (GKLS) Master Equation [61, 62, 63]. Our
strategy to design the proper dissipation is motivated by the recent
observation that training a local cost function (a cost function made of local
observables) is much more efficient than training a global one [64, 65, 66,
67]. Rigorous analyses have proven that local cost functions, unlike global
ones, are immune to barren plateaus in the case of shallow circuits [40, 39].
While formulating variational algorithms locally is preferable, mapping a
global problem to a local one is generally a non-trivial task. Here, we
propose variational algorithms based on non-unitary ansatzes and we show that
this represents an effective strategy to tackle barren plateaus (Sec. III),
setting the theory framework (Sec. IV) and discussing both a synthetic and a
chemical example (Sec. V and Sec. VI).

##  II General framework

Figure 1: Cross sections of the warm-up example cost function landscapes of
Sec. III.1 for a 20-qubit system. (a) Fully unitary ansatz $C_{u}$, (b) noisy
landscape $C_{n}$ with depolarizing probability $p=0.5$, and (c) engineered-
dissipation ansatz $C_{ed}$ in the case of $\Delta t\simeq 2,33$, which
corresponds to the maximum of $\mathop{\rm Var}[\partial
C_{ed}/\partial\theta_{j}]$ (see main text).

In the most general case, VQAs consist of minimizing a cost function whose
minimum faithfully corresponds to the solution of a considered problem [5]. In
the following, we consider the case where the cost function corresponds to the
expectation value of an $n$-qubit Hermitian operator H and, consequently, the
solution is its ground-state energy. Given an initial condition $\rho_{in}$
and a quantum circuit ansatz $U(\bm{\theta})$, where $\bm{\theta}$ is a free
parameter vector, the cost function has the form

$\displaystyle
C(\bm{\theta})=\Tr\\{HU(\bm{\theta})\rho_{in}U^{\dagger}(\bm{\theta})\\}.$

The minimization strategy consists of evaluating $C(\bm{\theta})$ by a quantum
device. The free parameters are optimized by a classical procedure.

A necessary condition for the circuit trainability is that, for a random
initialization of $\bm{\theta}$, the probability of finding a non-negligible
value of the cost function derivative with respect to a generic parameter
$\theta_{k}$, denoted as $\partial_{k}C$, is itself not negligible. An upper
bound on this probability is known to be proportional to the variance of the
derivative, which we write as $\mathop{\rm Var}[\partial_{k}C]$, setting the
limits to training efficiency [29, 34, 35, 36, 37, 38, 39, 40, 41, 42]. By
definition, we say that the cost function landscape presents a barren plateau
if $\mathop{\rm Var}[\partial_{k}C]$ exponentially decreases with the number
of qubits $n$, i.e., $\mathop{\rm Var}[\partial_{k}C]\in\mathcal{O}(e^{-pn})$
where $p$ is a positive integer.

This phenomenon strongly depends on the locality of $H$ where locality, in
this context, refers to the number of qubits on which the Hamiltonian acts
non-trivially. We expand the Hamiltonian such that

$\displaystyle H=c_{0}\mathbb{I}+\sum_{i=0}^{N}c_{i}H_{i}$ (1)

where $\mathbb{I}$ is the $n$-qubit identity operator, $H_{i}$ is a generic
Hermitian operator and $c_{i}$ is a real coefficient. By definition, we say
that $H$ is local when all $H_{i}$ terms act non-trivially on at most $K$
qubits (where $K$ does not scale with $n$), as is the case of nearest-neighbor
interaction Hamiltonians. On the other hand, we call $H$ global if the $H_{i}$
operators act non-trivially on all qubits. It has been rigorously proved that,
under quite general conditions, for alternating layered ansatzes of arbitrary
depth, a barren plateau is inevitable in the global case, while in the local
case shallow circuits can prevent it [39, 40]. More precisely, if the number
of layers (L) does not increase faster than a logarithm in the number of
qubits, i.e., $L=\mathcal{O}(\log(n))$ then:

$\displaystyle\mathop{\rm
Var}[\partial_{k}C]=\Omega\left(\frac{1}{\text{poly}(n)}\right)$

which means that the probability of finding a non-negligible gradient does not
decrease faster than a polynomial implying the absence of barren plateaus.

It has also been empirically found that when a problem can be encoded with
both a global and a local cost function, the training is more effective in the
local case and the mitigation techniques work much better [34, 64, 65, 66].
Moreover, it has been also analytically shown that $\mathop{\rm
Var}[\partial_{k}C]$, in general, increases with the Hamiltonian locality [41,
42]. In the following, we will focus on the case where the Hamiltonian of the
problem is originally global and an equivalent local one is unknown. We will
show how the addition of a proper non-unitary layer in the variational scheme
allows the problem to be approximated with a local one where barren plateaus
are absent or easier to face.

##  III Non-unitary Ansatz

We propose to generalize the unitary framework for VQA to Markovian maps,
designing the proper dissipation [59, 60] amenable to experimental
implementation, also in near-term quantum devices. Previous non-unitary ansatz
in VQA’s have been considered for in silico implementations (classical post-
processing) to boost precision in chemical problems [53, 54]. The non-unitary
ansatz acting on the quantum state is

$\displaystyle\Phi(\bm{\sigma},\bm{\theta})\rho=\mathcal{E}(\bm{\sigma})\circ
U(\bm{\theta})\rho U^{\dagger}(\bm{\theta})$

where $U(\bm{\theta})$ is, as before, a parametric quantum circuit, and
$\mathcal{E}(\bm{\sigma})$ is a non-unitary superoperator, with their
respective tunable parameters $\bm{\theta}$ and $\bm{\sigma}$. Under Markovian
assumption, the non-unitary part of the ansatz is defined by a parametric
Liouvillian $\mathcal{L}=\mathcal{L}(\bm{\sigma})$ such that

$\displaystyle\mathcal{E}(\bm{\sigma})=e^{\mathcal{L}(\bm{\sigma})\Delta t},$

where $\Delta t$ is the interaction time with the environment. The Liouvillian
$\mathcal{L}$ is the superoperator that generates the dynamics of the GKLS
Master Equation [61, 62, 63]:

$\dot{\rho}=\mathcal{L}\rho\equiv-i[\mathfrak{H},\rho]+\sum_{i}\gamma_{i}(L_{i}\rho
L_{i}^{\dagger}-\frac{1}{2}\\{L_{i}^{\dagger}L_{i},\rho\\})$ (2)

where $\mathfrak{H}$ is the Hamiltonian responsible of the unitary part of the
evolution, $\\{\gamma_{i}\\}$ are the damping rates and the operators
$\\{L_{i}\\}$, called jump operators, identify the environment action on the
qubits.

For our specific proposal, we set the following conditions on $\mathcal{L}$:

  1. 1.

$\mathcal{L}$ can be expressed as the sum of $Q$ superoperators
$\mathcal{L}_{q}(\bm{\sigma}_{q})$, each of which acts non-trivially on at
most $K$ qubits. This expansion takes the form:

$\displaystyle\mathcal{L}(\bm{\sigma})=\sum_{q=1}^{Q}\mathcal{L}_{q}(\bm{\sigma}_{q}).$
(3)

  2. 2.

All the generators $\mathcal{L}_{q}$ commute with each other.

  3. 3.

All the generators $\mathcal{L}_{q}$ have exactly one stationary state,
$\rho_{ss,q}$.

  4. 4.

The generators $\mathcal{L}_{q}$ converge to their respective stationary
states at the same rate, which we refer to as the mixing time.

According to the previous definitions and under the assumptions 1. and 2., the
cost function reads:

$\displaystyle C(\bm{\theta},\bm{\sigma})$
$\displaystyle=\Tr{He^{\mathcal{L}(\bm{\sigma})\Delta
t}U(\bm{\theta})\rho_{in}U^{\dagger}(\bm{\theta})}$
$\displaystyle=\Tr{H\prod_{q=1}^{Q}e^{\mathcal{L}_{q}(\bm{\sigma}_{q})\Delta
t}U(\bm{\theta})\rho_{in}U^{\dagger}(\bm{\theta})}.$ (4)

###  III.1 Illustrative example

Before delving into the formal characterization, we present a simple example
to illustrate the strategy of tackling the barren plateaus with non-unitary
VQAs. Following [39], we consider a state preparation task formulated through
the optimization of a global Hamiltonian. In particular, we are interested in
minimizing $H=\mathbb{I}-\ket{\bm{0}}\bra{\bm{0}}$ supposing that all the
qubits are initialized in state $\ket{0}$. By applying the unitary ansatz
$U(\bm{\theta})=\bigotimes_{j=1}^{n}e^{-i\theta_{j}\sigma_{x}/2}$, we can
easily detect the occurrence of a barren plateau (see Fig. 1 a). In fact, the
cost function takes the form
$C_{u}=1-\prod_{j=1}^{n}\cos^{2}(\frac{\theta_{j}}{2})$ and a direct
calculation shows that $\mathop{\rm Var}[\frac{\partial
C_{u}}{\partial\theta_{j}}]=\frac{1}{8}(\frac{3}{8})^{n-1}$. In addition, the
derivative values are unbiased, i.e. $\left\langle\frac{\partial
C_{u}}{\partial\theta_{j}}\right\rangle=0$, implying an exponential
suppression of the probability to sample a non-negligible gradient as a
consequence of Chebyshev’s inequality.

We now emphasize that, as shown in [37], general noisy non-unitary
interactions are not able to mitigate barren plateaus and can actually induce
them [37]. For example, if we model the noise effect with a depolarizing
channel of probability $p$, the resulting cost function will be
$C_{n}=p(1-\frac{1}{2^{n}})+(1-p)C_{u}$ (see Fig. 1 b). This non-unitary
model, in addition to clearly having the same trainability problems of the
cost function $C_{u}$, also precludes the possibility of finding the correct
ground state energy.

We now introduce an engineered dissipation through a non-unitary layer,
ensuring that each $\mathcal{L}_{q}$ operator acts dissipatively on a single
qubit as per Eq. (2). We assume that, in general, the $\mathcal{L}_{q}$
structure lacks a Hamiltonian part and instead is determined by a single jump
operator $L_{q}=\ket{0}_{q}\bra{1}$ with a corresponding damping rate
normalized to $1$. This situation is advantageous because the stationary state
of the whole Liouvillian is the ground state of the sought-after problem. The
cost function in the presence of such an engineered dissipation is now
$C_{ed}=1-\prod_{j=1}^{n}[1-\sin^{2}(\frac{\theta_{j}}{2})e^{-\Delta t}]$ for
which the unbiased condition on the gradient still holds. However, in this
case, it is possible to efficiently avoid the barren plateau (see broader and
still deep minimum in Fig. 1 c). In fact, $\mathop{\rm Var}[\frac{\partial
C_{ed}}{\partial\theta_{j}}]=\frac{1}{8}e^{-\Delta t}(1+\frac{3}{8}e^{-2\Delta
t}-e^{-\Delta t})^{n-1}$ and if $\Delta t\sim\mathcal{O}(\log(n))$ the desired
polynomial scaling is achieved.

We observe that the $C_{ed}$ landscape takes a similar shape as the one in
[39], where the unitary ansatz was the same as we have considered here, but
the global Hamiltonian was replaced by an equivalent local one. In the
following, we will show how this landscape similarity is not a coincidence
because the considered non-unitary ansatz allows, in general, localizing a
problem that has been originally formulated as a global one. Moreover, we will
also prove that the logarithmic $\Delta t$ scaling is a general feature of the
procedure.

##  IV Theoretical results

Now we will show how a non-unitary layer that respects the constraints
introduced in Sec. III is able to make the Hamiltonian of the cost function
local. First, we note that the diagonalizability of the $\mathcal{L}_{q}$
superoperators is a mild condition that can be assumed to hold. In fact, if
$\mathcal{L}_{q}$ has a holomorphic dependence on the parameters
$\bm{\sigma}_{q}$ and if it is diagonalizable for a subspace of them, then it
will be diagonalizable in general, apart from the possible appearance of
countable exceptional points, which can be ignored in the discussion [68].
Anyway, in all the cases discussed in the following, diagonalizability will
always be ensured. Consequently, it is useful to expand the $\mathcal{L}_{q}$
superoperators through their dual basis. To keep the notation simple, we will
assume that each $\mathcal{L}_{q}$ acts non-trivially on exactly $K$ qubits
111The following calculations can be immediately generalized to cases where
each $\mathcal{L}_{k}$ operator acts on a generic number of qubits less than
or equal to $K$.. Introducing the Liouville notation for a generic matrix,
$\rho\rightarrow|\rho\rangle\\!\rangle$, and considering the dot product
$\langle\\!\langle\tau|\rho\rangle\\!\rangle=\Tr{\tau^{\dagger}\rho}$, then

$\displaystyle e^{\mathcal{L}_{q}(\bm{\sigma}_{q})\Delta
t}=\sum_{i=0}^{4^{K}-1}e^{\lambda_{q,i}\Delta
t}\frac{|r_{q,i}\rangle\\!\rangle\langle\\!\langle l_{q,i}|}{\langle\\!\langle
l_{q,i}|r_{q,i}\rangle\\!\rangle},$ (5)

where $\\{\lambda_{q,i}\\}_{i}$ is the set of the eigenvalues of
$\mathcal{L}_{q}$ and $\\{|r_{q,i}\rangle\\!\rangle\\}_{i}$ and
$\\{|l_{q,i}\rangle\\!\rangle\\}_{i}$ are the corresponding set of orthogonal
and normalized right and left eigenvectors 222we have omitted the
$\bm{\sigma}_{q}$ dependence of all these terms to make the notation easier.
Ordering the indexes as a function of the real parts of the eigenvalues such
that
$\Re{\lambda_{q,0}}>\Re{\lambda_{q,1}}\geq\dots\geq\Re{\lambda_{q,4^{K}-1}}$,
and assuming the uniqueness of the steady state
$|\rho_{ss,q}\rangle\\!\rangle$ (this is also a mild condition to satisfy), we
identify $\lambda_{0}=0$,
$|r_{q,0}\rangle\\!\rangle=|\rho_{ss,q}\rangle\\!\rangle$ and
$|l_{q,0}\rangle\\!\rangle=|\mathbb{I}\rangle\\!\rangle$.

Equation 5 can be used to calculate the order of magnitude of the mixing times
$\Delta t_{mix,q}$. Indeed, excluding the occurrence of singular phenomena
such as the skin effect [71], we have $\Delta
t_{mix,q}\sim\mathcal{O}(1/|\Re{\lambda_{q,1}}|)$, where the quantity
$|\Re{\lambda_{q,1}}|$ represents the spectral gap. Our constraint of a common
mixing time for all $\mathcal{L}_{q}$ is, then, satisfied if the spectral gaps
are equal and do not vary with $\bm{\sigma}_{q}$. Then, we deal with a unique
mixing time, which we will refer to as $\Delta t_{mix}$.

Defining
$|\rho(\bm{\theta})\rangle\\!\rangle\equiv|U(\bm{\theta})\rho_{in}U^{\dagger}(\bm{\theta})\rangle\\!\rangle$,
we can rewrite Eq. (4) with the introduced notation:

$\displaystyle C(\bm{\theta},\bm{\sigma})=\langle\\!\langle
H|e^{\mathcal{L}(\bm{\sigma})\Delta t}|\rho(\bm{\theta})\rangle\\!\rangle$
$\displaystyle=\langle\\!\langle
H|\sum_{i_{1},\dots,i_{Q}=0}^{4^{K}-1}e^{(\sum_{q=1}^{Q}\lambda_{q,i_{q}})\Delta
t}\bigotimes_{q=1}^{Q}\frac{|r_{q,i_{q}}\rangle\\!\rangle\langle\\!\langle
l_{q,i_{q}}|}{\langle\\!\langle
l_{q,i_{q}}|r_{q,i_{q}}\rangle\\!\rangle}|\rho(\bm{\theta})\rangle\\!\rangle.$

For a value of the time such that $\Delta t\apprle\log(Q)\cdot\Delta t_{mix}$
an approximation of the non-unitary superoperator, in which only the slowest
decay terms are taken, is justified (Theorem 6 in [72]):

$\displaystyle e^{\mathcal{L}(\bm{\sigma})\Delta
t}=\bigotimes_{q=1}^{Q}\frac{|\rho_{ss,q}\rangle\\!\rangle\langle\\!\langle\mathbb{I}|}{\langle\\!\langle\mathbb{I}|\rho_{ss,q}\rangle\\!\rangle}$
$\displaystyle+\sum_{q=1}^{Q}e^{\lambda_{q,1}\Delta
t}\frac{|r_{q,1}\rangle\\!\rangle\langle\\!\langle l_{q,1}|}{\langle\\!\langle
l_{q,1}|r_{q,1}\rangle\\!\rangle}\bigotimes_{j\neq
q}^{Q}\frac{|\rho_{ss,j}\rangle\\!\rangle\langle\\!\langle\mathbb{I}|}{\langle\\!\langle\mathbb{I}|\rho_{ss,j}\rangle\\!\rangle}$
$\displaystyle+\mathchoice{{\scriptstyle\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{{\scriptscriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptscriptstyle\mathcal{O}$}}\left(e^{-\Delta
t/(\log(Q)\cdot\Delta t_{mix})}\right).$ (6)

Moving to the Heisenberg picture, the variational problem, defined by the cost
function of Eq. (4), can be formulated as a unitary optimization on an
effective Hamiltonian $H^{{}^{\prime}}$ such that:

$\displaystyle H^{{}^{\prime}}=e^{\mathcal{L^{\dagger}}(\bm{\sigma})\Delta
t}H.$ (7)

Writing $H$ as a linear combination of the $\mathcal{L}$ left eigenvectors
(which form a complete set)

$\displaystyle
H=\sum_{j_{1},\dots,j_{Q}}c_{j_{1},\dots,j_{Q}}\bigotimes_{q=1}^{Q}|l_{q,i_{q}}\rangle\\!\rangle$

and assuming that higher order terms in Eq. (IV) can be disregarded, we arrive
at the following approximate expression for $H^{{}^{\prime}}$:

$\displaystyle H^{{}^{\prime}}\simeq
c^{\prime}_{0}(\bm{\sigma})\mathbb{I}+\sum_{q=1}^{Q}c^{\prime}_{q}(\bm{\sigma})|\mathbb{I}_{q}\rangle\\!\rangle\otimes|l_{q,1}\rangle\\!\rangle,$
(8)

where $\mathbb{I}_{q}$ refers to the identity operator in all qubit Hilbert
spaces except the $q$-th subspace. At this point, we observe that
$H^{{}^{\prime}}$ takes the local form of Eq. (1) and is Hermitian, because
the time evolution of $H$, according to (7), can always be written as a sum of
time-decaying Hermitian matrices 333It is a direct consequence of the
following Liovillian property: if $\rho$ is a left eigenvector of
$\mathcal{L}$ (or a right one of $\mathcal{L}^{\dagger}$) with eigenvalue
equal to $\lambda$ then $\rho^{\dagger}$ will be in turn an eigenvector with
$\lambda^{*}$ as an eigenvalue. This statement can be proved easily by
applying the proof of Lemma 3 in [68] considering the adjoint of a generic
Liouvillian.. Therefore, neglecting terms in the expansion does not preclude
the hermiticity of the operator. Consequently, all the results concerning the
training of the hyperparameter vector $\bm{\theta}$ for local Hamiltonians can
be applied. We emphasize that the condition imposed on the mixing times
guarantees that each $\mathcal{L}_{q}$ operator contributes equally to the
expansion of the slowest terms of Eq. IV, resulting in the definition of
$H^{\prime}$. It is also easy to verify that the emergence of an effective
Hamiltonian with local character holds even when faster modes are included.

Our strategy is effective only if the $c^{\prime}$ coefficients are not
negligible. This property is achieved when the slowest decay right
eigenvectors of $\mathcal{L}$ have a significant overlap with $H$ and,
consequently, the choice of the Liouvillan has to be designed according to the
known Hamiltonian properties.

In addition to avoiding the problem of the barren plateau, we have to make
sure that the minimum energy of $H^{{}^{\prime}}$ and the minimum energy of
$H$ are close enough. Such a closeness will depend, in general, on the value
of the parameters $\bm{\sigma}$, and, consequently, the trainability of the
non-unitary layer is a property to take into account. To this end, we study
how the derivative of the cost function varies as a function of a generic
component of $\bm{\sigma}$, which we call $\sigma_{k}$. As in the previous
case of the unitary parameters, a key quantity to calculate is the variance of
$\sigma_{k}$, which is related to the probability of sampling a non-negligible
value of such parameter. If we call $d\mu(\bm{\theta})$ and
$d\mu(\bm{\sigma})$ the distribution volume elements of the free parameters,
we obtain

$\displaystyle\mathop{\rm Var}\left[{\frac{\partial
C}{\partial\sigma_{k}}}\right]$ $\displaystyle=\mathop{\rm
Var}\left[\Tr{\frac{\partial
H^{{}^{\prime}}}{\partial\sigma_{k}}\rho(\bm{\theta})}\right]$
$\displaystyle=\mathop{\rm
Var}\left[\Tr{\tilde{H}(\bm{\sigma})\rho(\bm{\theta})}\right]$
$\displaystyle=\int_{\bm{\sigma}}d\mu(\bm{\sigma})\int_{\bm{\theta}}d\mu(\bm{\theta})\left(\tilde{C}-\left\langle\tilde{C}\right\rangle\right)^{2}$
(9)

where $\tilde{H}$ is a local Hamiltonian that can be written in the form of
Eq. (8) and $\tilde{C}$ is its relative cost function computed with respect to
the unitary ansatz considered.

From Eq. (9) we learn that an exponential suppression of $\mathop{\rm
Var}[{\partial C/\partial\sigma_{k}}]$), according to [36], can occur if and
only if the unitary ansatz presents a barren plateau in the optimization of
$\tilde{H}$. Since $\tilde{H}$ and $H^{\prime}$ have the same local character,
considering a quantum circuit that is not affected by barren plateaus in the
case of local Hamiltonians (see Sec. II) directly excludes the presence of
barren plateaus in the non-unitary layer.

Moreover, we want to remark that the value of $\Delta t$ that ensures the
validity of the approximation in Eq. (IV) scales efficiently with the number
of generators that compose the model, that is, $\mathcal{O}(\log(Q))$,
independent on the specific Hamiltonian $H$ of the problem. Furthermore, the
condition $Q=\mathcal{O}(\text{poly}(n))$ ensures an efficient logarithmic
scaling even in the number of qubits, as required.

Let us also emphasize that the condition of Eq. (IV) can be alleviated by
including faster decay terms while still preserving a local structure for
$H^{{}^{\prime}}$. This means that for practical purposes we can choose a
value of $\Delta t$ significantly smaller than the theoretical threshold
discussed above.

Finally, our analysis extends seamlessly to considering a convex combination
of non-unitary maps generated by a Liouvillian that satisfies the constraints
introduced in Sec. III. Indeed, a non-unitary layer given by the convex
combination $\mathcal{E}=\sum_{j}\beta_{j}e^{\mathcal{L}_{j}\Delta t}$, in
addition to being a suitable quantum channel, will still transform a global
Hamiltonian to a local one. Importantly, if we compute the variance of the
derivative of the cost function with respect to either the $\beta_{j}$
coefficients or the $\mathcal{L}_{j}$ free parameters, we find that the result
reduces to the form shown in Eq. 9. This observation generalizes the
conclusions drawn from our previous analysis.

##  V Random Hamiltonian example

We will now provide a numerical demonstration of the scaling calculated
analytically in Sec. IV in a synthetic example. We consider a random
Hamiltonian whose ground state in each realization can be either in the
neighborhood of the state $\ket{\bm{0}}$ or of $\ket{\bm{1}}$. Assuming lack
of knowledge of the ground state, the initial guess for the ansatz is random.
For the sake of definiteness, we have fixed the minimum energy to the value of
-1.1 (see Appendix A for more details). Moreover, as in Refs. [29, 34], we
will employ a hardware-efficient, layered quantum circuit as the unitary
component of our ansatz (its form is shown in Appendix B). Considering the
known properties of the Hamiltonian, we define the non-unitary layer as
follows:

$\displaystyle\mathcal{E}(\sigma)=s(\sigma)e^{\mathcal{L}^{(1)}\Delta
t}+\left[1-s(\sigma)\right]e^{\mathcal{L}^{(2)}\Delta t},$

where $\sigma$ is a real free parameter, $s(\sigma)$ is the sigmoid function
that ensures the sum is a convex combination, and $\mathcal{L}^{(1)}$ and
$\mathcal{L}^{(2)}$ are Liouvillians made up of single-qubit dissipators such
that their corresponding stationary states are $\ket{\bm{0}}\bra{\bm{0}}$ and
$\ket{\bm{1}}\bra{\bm{1}}$, respectively. These two Liouvillians can be easily
derived considering the superoperators family identified in Appendix C in such
a way that they will respect the needed constraints presented in Sec. III.
This particular choice ensures a non-negligible overlap with the ground state
of the target Hamiltonian. Moreover, by training the parameter $\sigma$, we
can transform the outputs of the quantum circuit into states close to the
desired ground state. This convex combination can be implemented by
considering a stochastic Liouvillian, sampling $\mathcal{L}^{(1)}$ with
probability $s(\sigma)$ and $\mathcal{L}^{(2)}$ with complementary
probability.

Figure 2: Results of the random Hamiltonian example. (a) Scaling of the
partial derivative variances with respect to $\theta_{1}^{1}$ and $\sigma$.
Each point was determined from a statistics of $1000$ different random
samples. It includes different rotations directions, Hamiltonian and free
parameter realizations. Moreover, the $x$ parameter was selected at random in
the range $[-5,5]$ while all the angles were also selected at random in the
compact set $[0,2\pi]$. For the non-unitary ansatz, we varied $\Delta t$ from
$0.1$ to $3$ in steps of $0.1$, selecting the biggest variance for each value
of $n$. The variances with respect to $\theta_{1}^{1}$ are depicted by dashed
lines for the fully unitary ansatz whereas, for the non-unitary ansatz case,
they are represented by continuous lines. The partial derivative variance with
respect to $x$ is depicted using the black dotted lines in the upper part of
the plot. We only plot the cases where the non-unitary ansatz is beneficial
for training (for very short chains no improvement is observed compared to the
unitary ansatz ). (b) Variance averages with respect to $\theta_{1}^{1}$ as a
function of $\Delta t$ in the case of a 5-layer quantum circuit, for 5 and 8
qubits. (c)-(d) Evolution of the cost function evolution under the gradient
descent algorithm with a learning rate of $0.1$ and for $10$ random initial
conditions. The light lines represent single realizations, while the relative
averages are reported in a darker color. (c) The performances of the two
different cases (unitary vs non-unitary ansatz) are compared. (d) Cost
function calculated through a hybrid approach: the first 500 iterations
involve non-unitary ansatzes, while in the remaining 500 we impose the
condition $\Delta t=0$.

In our numerical experiment, we first consider computing the derivative
variance scaling with respect to the unitary ansatz parameters, with and
without the addition of the non-unitary layer. In the case of the non-unitary
ansatz, we have spanned $\Delta t$ from $0.1$ to $3$ with a resolution of
$0.1$ to find its optimal values. In Fig. 2(a) we can clearly observe that,
for a fixed number of layers of the quantum circuit, dissipation enables the
prevention of an exponential scaling of the gradient variance, which indicates
the absence of the barren plateau, as we predicted. We also remark that the
optimal values of $\Delta t$ found are always below the characteristic time
threshold indicated in Sec. IV. To give an idea of the role of the dissipation
strength, we plot the trend of the variance as a function of $\Delta t$ in
Fig. 2(b). We find a non-monotonic trend that is the result of the competition
between the constructive role of dissipation, leading to an increase of
locality, but also the drawbacks of a contractive dynamics. The latter at long
time would lead to a full Hamiltonian erasure (large $\Delta t$ in Fig. 2 2
(b)), while in the initial transient can rescale the Hamiltonian with a
consequent reduction of the variance (more visible for 5-qubits). For the non-
unitary layer to provide an advantage, the number of qubits has to be large
enough to obtain a maximum that overcomes the initial value, which corresponds
to the fully unitary case. In Fig. 2(a), we also show the behavior of the
derivative of the variance with respect to the free parameter of the unitary
layer as the number of qubits is changed. We can clearly exclude the presence
of an exponential scaling, in agreement with the theoretical results of Sec.
IV.

Finally, it is also important to test the accuracy of the ground-state energy
estimation for both the unitary and the non-unitary ansatz. To this end, we
applied $1000$ iterations of the gradient descent algorithm, with a learning
rate of $0.1$, to $10$ random initial conditions in the ideal case where the
gradient can be exactly estimated. The results, presented in Fig. 2(c), show
that the non-unitary circuit allows a faster convergence, with a relative
average error of $2.2\%$, while the fully unitary circuit has a slower
convergence but a higher accuracy of $1.5\%$. This is expected since the fully
unitary circuit preserves the original Hamiltonian and hence its ground state.

To take advantage of both properties, we propose a hybrid approach where the
non-unitary ansatz is applied to the initial state for the first $500$
iterations, and then only the fully unitary circuit is used for the last $500$
iterations by setting $\Delta t$ to zero. This method serves as a novel
initialization strategy, and from Fig. 2 (d), we observe a significant
improvement in convergence time compared to the fully unitary case.
Additionally, the average accuracy is the best found with a value of $1.2\%$.

##  VI Quantum chemistry example

Figure 3: Results of the gradient descent algorithm for the example of
$H_{2}$. (a) First 150 iterations of the algorithm. The blue line refers to
the non-unitary ansatz optimization with $\Delta t=0.5$ and a learning rate
equal to 1, the red dotted one to the unitary-ansatz optimized with a learning
rate equal to 1, while, for the red continuous one, the unitary ansatz takes a
learning rate equal to 0.1. The vertical black line indicates the ground-state
energy computed by a numerical diagonalization of the Hamiltonian. (b) Last
150 iterations. The green line shows the hybrid ansatz optimization and the
grey region refers to energy values falling within the considered error
threshold. The hybrid ansatz learning rate is set to 0.1.

Moving to a more practical application, we show the results of the beneficial
effects of our non-unitary model in the context of quantum chemistry. In
particular, we have studied the problem of finding the electronic ground state
of the $H_{2}$ molecule. In our numerical experiments, we studied the molecule
in its equilibrium configuration where the two hydrogen atoms are separated by
$0.74$ Å. The qubit Hamiltonian $H$ was obtained from a Jordan-Wigner
transformation [74], while the electronic orbitals were generated from the
STO-3G basis set [75]. Under these conditions, $H$ spans a four-qubit Hilbert
space. Following the previous example, the quantum circuit that defines the
unitary part of the ansatz is the one presented in Appendix 4 with a number of
layers fixed to $20$. As for the non-unitary layer, the nature of the problem
suggests choosing a Liouvillian whose corresponding stationary state is the
Hartree-Fock state [76]. Interestingly, the class of Liouvillians presented in
Appendix C can always accomplish this task by simply dissipating along the $z$
direction, because the Hartree-Fock state is generally a computational basis
state.

As in the previous example of Sec. V, we have applied the gradient descent
algorithm to ten random initial conditions both in the case of a fully unitary
ansatz, which uniquely consists of a quantum circuit, and in the case of a
non-unitary one. We found that a value of $\Delta t=0.5$ is sufficient to
significantly improve the convergence time and the ground state quality as
shown in Fig. 3. In particular, we have found that in the non-unitary case,
the algorithm is able to converge by fixing a learning rate of $1$, while in
the fully unitary case, we need to set the learning rate to $0.1$ to achieve
faster convergence time and a lower value of the final ground state. In all
the considered examples, we have fixed the number of iterations to 300. As in
the random Hamiltonian example, we observe that the introduction of the non-
unitary layer makes it possible to significantly reduce the algorithm time of
almost two orders of magnitude, at the cost of convergence to a less accurate
ground state.

As done in the previous section, we propose a hybrid approach in which, we
first apply the non-unitary ansatz, and after half of the total iteration
steps, we remove the dissipation and apply only the unitary ansatz to complete
the experiment. As shown in Fig. 3, this strategy is successful both in
reducing the convergence time and in improving the quality of the final ground
state. We have set an error threshold with respect to the ground state energy
calculated using a numerical diagonalization equal to 0.00159 Hartree, in line
with the standard threshold for defining chemical accuracy. Notably, only the
hybrid approach can estimate the ground state energy within this energy
interval.

##  VII Discussion

Extensive research in the field of variational quantum algorithms is devoted
to exploring strategies and techniques aimed at circumventing or minimizing
the occurrence of barren plateaus. This is essential for improving
scalability, algorithmic efficiency, and applicability in a wide range of
contexts. Here we have proposed and demonstrated the effectiveness of a
strategy based on engineered dissipation, that can be understood as a
localization of a global cost function argument. This dissipation-based
analysis goes beyond the general scaling laws recently presented in Refs. [41,
42], which are limited to the use of a unitary ansatz. Indeed, the presented
mitigation strategy has a time complexity that scales logarithmically as the
number of qubits is increased, which makes it possible to estimate the ground
state energy of systems that would be intractable with previous methods.

To assess the effectiveness of our approach, we have first examined the
disappearance of the barren plateau resulting from the presence of engineered
dissipation in a simple illustrative model that admits an analytical solution.
This toy model also serves as a useful point of comparison for the scenario
with general noise, which worsens the problem of attenuation of the barren
plateau. Interestingly, handling Hamiltonians such as
$H=\mathbb{I}-\ket{\psi}\bra{\psi}$, where $\ket{\psi}$ represents a generic
pure state, can be challenging with alternative proposals, such as
perturbative gadgets [77]. While both these approaches map a global
Hamiltonian into a local one, our engineered dissipation strategy does not
require to increase the Hilbert space dimension and the time complexity is
entirely unrelated to the decomposition of the Hamiltonian using Pauli
matrices.

The formal theoretical framework of this proposal stands on the analytical
proofs in Sec. IV, where we first demonstrate that the presence of a tailored
non-unitary layer makes it possible to map the cost function into an effective
one that corresponds to a local Hamiltonian. Moreover, we also prove that this
type of non-unitary layer can be efficiently trained. This theoretical
proposal opens a new avenue for implementations in noisy quantum processors.
Actually, the form of engineered dissipation here presented can be efficiently
implemented on experimental platforms, for instance through collision models
[78], as discussed in Appendix D.

We then tested the use of nonunitary ansatzes in two relevant numerical
examples. The first one addressed a random synthetic Hamiltonian, while the
second tackled the realistic task of determining the lower energy of a
Hydrogen molecule. Both problems display the effectiveness of our method in
speeding up the convergence and also the precision reached in the estimation
of the ground state energy. Indeed, the non-unitary approach can also be seen
as an efficient initialization protocol for a unitary ansatz suggesting a
powerful hybrid strategy where a unitary ansatz follows a non-unitary one.

Our results are consistent with other quantum machine learning protocols,
where the presence of losses can enhance performance. In a recent work, a
dissipative optimization algorithm, able to efficiently find Hamiltonian local
minima, has been reported [79]. Another relevant example can be found in the
realm of quantum reservoir computing [80], where dissipation can be
transformed into a constructive resource [81, 82, 83]. To stay within the
context of VQAs, it was demonstrated that incorporating stochastic noise can
prevent the occurrence of saddle points, which are detrimental to efficient
optimization [84]. Finally, we proposed a simple dissipation model that proved
very effective in mitigating barren plateaus. This indicates the potential for
devising more comprehensive strategies that fully explore the potential of
non-unitary architectures in variational quantum circuits but also in the
broader arena of quantum neural networks.

## Acknowledgements

We acknowledge the Spanish State Research Agency, through the María de Maeztu
project CEX2021-001164-M funded by the MCIN/AEI/10.13039/501100011033 and
through the QUARESC project (PID2019-109094GB-C21/AEI/10.13039/501100011033),
MINECO through the QUANTUM SPAIN project, and EU through the RTRP -
NextGenerationEU within the framework of the Digital Spain 2025 Agenda. We
also acknowledge funding by CAIB through the QUAREC project (PRD2018/47). The
CSIC Interdisciplinary Thematic Platform (PTI) on Quantum Technologies in
Spain is also acknowledged. GLG is funded by the Spanish Ministerio de
Educación y Formación Profesional/Ministerio de Universidades and co-funded by
the University of the Balearic Islands through the Beatriz Galindo program
(BG20/00085). The project that gave rise to these results received the support
of a fellowship from the ”la Caixa” Foundation (ID 100010434). The fellowship
code is LCF/BQ/DI23/11990081. IBM, the IBM logo, and ibm.com are trademarks of
International Business Machines Corp., registered in many jurisdictions
worldwide. Other product and service names might be trademarks of IBM or other
companies. The current list of IBM trademarks is available at
https://www.ibm.com/legal/copytrade.

##  Appendix A Random Hamiltonian structure

In Sec. V, we have proposed the optimization of a random Hamiltonian $H$ from
which we have only partial information about a possible ground state. We will
now show the criteria used to randomly generate $H$ in such a way that it is
constrained to this incomplete ground state knowledge. We point out that the
non-unitary optimization algorithm shown in Sec. V only has access to this
Hamiltonian property and, as a consequence, the key idea can be applied in a
general context.

Going now into the details of $H$, for the sake of definiteness, we will
impose the condition that all its eigenvalues $\lambda_{i}$ satisfy
$|\lambda_{i}|<1$, except for the maximum eigenvalue which is fixed at 1.1,
and the minimum eigenvalue which is set at -1.1. Additionally, we have
constrained the eigenvectors associated with these last two values to be one
of two states (up to a normalization factor):
$\ket{\psi_{1}}=\ket{\bm{0}}+0.1\ket{\phi_{1}}$ and
$\ket{\psi_{2}}=\ket{\bm{1}}+0.1\ket{\phi_{2}}$, selected randomly and
exclusively, where $\ket{\phi_{1}}$ and $\ket{\phi_{2}}$ are Haar-distributed
states. We stress that for $\ket{\psi_{1}}$ and $\ket{\psi_{2}}$ to be, in
general, viable eigenvectors, they must undergo a Gram-Schmidt
orthonormalization.

Finally, we underline that the algorithm of Sec. V is built from two possible
ground-state neighborhood guessing which are referred to two opposite cases:
the maximum and the minimum eigenvalue. Consequently, the numerical results
obtained indicate that the algorithm was able to successfully discriminate
them.

##  Appendix B Quantum circuit

Figure 4: Quantum circuit employed in the numerical simulations of Sections V
and VI.

We now show in detail the form of the quantum circuit considered in Sections V
and VI. Following the refs. [29, 34], it has the structure of a hardware-
efficient, layered quantum circuit structured as follows:

$U(\bm{\theta},\bm{d})=\prod_{l=1}^{D}WU_{l}(\bm{\theta}_{l},\bm{d}_{l}).$

This comprises $D$ layers of rotations and entangling gates, with the
entangler $W$ consisting of a series of controlled-$Z$ gates that correlate
adjacent qubits:

$W=\prod_{i=1}^{n-1}CZ_{i,i+1}.$

A $U_{l}$ gate, on the other hand, is formed by single-qubit rotations in
random directions:

$U_{l}(\bm{\theta}_{l},\bm{d}_{l})=\prod_{i=1}^{n}R_{d_{l}^{i}}(\theta_{l}^{i}),$

where $R_{d_{l}^{i}}(\theta_{l}^{i})$ is a rotation gate that applies an angle
of $\theta_{l}^{i}$ around the direction $d_{l}^{i}$ to the $i$-th qubit, and
$d_{l}^{i}$ is randomly chosen from the $x$, $y$, or $z$ axis, while
$\theta_{l}^{i}$ takes values in the set $[0,2\pi]$. To avoid any preferential
direction, the qubits are initially prepared in the state
$\rho_{in}=\ket{\psi_{0}}\bra{\psi_{0}}^{\otimes n}$, where
$\ket{\psi_{0}}=R_{y}(\frac{\pi}{4})\ket{0}$. The all circuit structure is
summarized in fig. 4.

##  Appendix C Liouvillians based on single-qubit dissipators

We now present a simple class of Liouvillians that satisfy the constraints
required for our purposes (see Sec. III) from which we have built the non-
unitary layer considered in sections V and VI. We consider a generator
$\mathcal{L}$ that takes the form of Eq. (3), where each operator
$\mathcal{L}_{q}$ implements single-qubit dissipation along a desired
direction. More precisely, their action is described by the following
relation:

$\displaystyle\mathcal{L}_{q}\rho=d_{q}\rho
d_{q}^{\dagger}-\frac{1}{2}\\{d_{q}^{\dagger}d_{q},\rho\\},$

where $d_{q}$ is a jump operator acting onto the q-th qubit Hilbert space:

$\displaystyle
d_{q}=\ket{\psi_{-}(\alpha_{q},\phi_{q})}_{q}\bra{\psi_{+}(\alpha_{q},\phi_{q})},$

with

$\displaystyle\ket{\psi_{+}(\alpha_{q},\phi_{q})}$
$\displaystyle=\cos(\frac{\alpha_{q}}{2})\ket{0}+e^{i\phi_{q}}\sin(\frac{\alpha_{q}}{2})\ket{1},$
$\displaystyle\ket{\psi_{-}(\alpha_{q},\phi_{q})}$
$\displaystyle=\sin(\frac{\alpha_{q}}{2})\ket{0}-e^{i\phi_{q}}\cos(\frac{\alpha_{q}}{2})\ket{1}.$

Here the parameters $\alpha_{q}$ and $\phi_{q}$ identify the direction of the
dissipation. This choice satisfies the required in Sec. III because the
spectral gap of each $\mathcal{L}_{q}$ is identically equal to $1/2$.
Moreover, each $\mathcal{L}_{q}$ has as a unique stationary state
$\ket{\psi_{-}(\alpha_{q},\phi_{q})}_{q}\bra{\psi_{-}(\alpha_{q},\phi_{q})}$,
and this information can be used to study the overlap with the target
Hamiltonian to ensure the effectiveness of the approach.

##  Appendix D Non-unitary layer implementation

Nonunitary operations can be effectively realised on digital quantum computing
architectures in several ways. This task is in fact closely related to the
general problem of implementing quantum simulation algorithms for open systems
[85, 86, 87, 88, 89]. The most direct approaches employ, for instance,
dilation theorems allowing one to recast dissipative dynamics into a unitary
evolution on an extended system [90, 91, 92]. Alternatively, or in combination
with the latter, one could in principle make use of engineered dissipations
leveraging intrinsic qubit noise [59, 93, 94, 95], randomized schemes [96, 97]
or controlled classical environments [98].

To give a concrete example, we now present an efficient implementation of the
nonunitary layer by means of the so-called collision models (CMs) [99]. The
general idea of CMs is to approximate the Markovian dynamics given by Eq. (2)
by letting the system qubits sequentially interact in a unitary way with a set
of ancillary qubits. In particular, the following approximation is considered:

$\displaystyle e^{\mathcal{L}\Delta t}\approx(\phi_{\Delta t/M})^{M}$

which consists of the alternating M steps where a quantum channel
$\phi_{\Delta t/M}$ approximates the evolution for a finite time equal to
$\Delta t/M$. By definition,

$\displaystyle\phi_{\Delta t/M}[\rho]=\Tr_{a}\\{V(\Delta
t/M)\rho\otimes\rho_{a}V^{\dagger}(\Delta t/M)\\}$ (10)

where $\Tr_{a}$ indicates the partial trace over the ancillary qubits space,
$\rho_{a}$ is their state which has to be properly prepared at each step and
$V(\Delta t/M)$ is a unitary operator which, depending on the particular
dynamics of interest, non-trivially acts on the system and ancillary qubits
space.

Importantly, for local Liouvillians as required by the criteria of Sec. III,
the number of resources needed by the collision-model strategy always scales
efficiently [92]. More quantitatively, let us call $\epsilon$ the difference
between the exact Liouvillan dynamics and the $M$-step collision model:

$\|e^{\mathcal{L}\Delta t}-(\phi_{\Delta t/M})^{M}\|_{1\rightarrow
1}=\epsilon,$

where $\left\|\cdot\right\|_{1\rightarrow 1}$ indicates the trace norm. Then,
it can be shown that the number of gates necessary to have an error $\epsilon$
is always a polynomial function of the number of qubits $n$.

According to Ref. [78], the class of models defined in Appendix Sec. C can be
efficiently implemented by coupling each qubit with an ancillary one
initialized in the state $\ket{0}$. Going more into details, in the definition
of $\phi_{\Delta t/M}$, according to Eq. 10, we will have
$\rho_{a}=\ket{0}\bra{0}^{\otimes n}$ and $V(\Delta
t/M)=\prod_{q=1}^{n}V_{q}(\Delta t/M)$ with

$\displaystyle V_{q}(\Delta t/M)=\exp(-i\sqrt{\Delta
t/M}(d_{q}\sigma^{+}_{q^{(a)}}+h.c.))$

where the index $q^{(a)}$ refers to the q-th ancillary qubit. In view of a
possible experimental implementation, we emphasize as a positive note that
each unitary $V_{q}(\Delta t/M)$ can be easily realized through a single-qubit
gate and a CNOT gate and that the qubit reset operation, currently available
on IBM Quantum devices [100], allows the number of ancillary qubits to be
fixed to $n$ for the entire execution of the collision-model algorithm.

## References

  * Preskill [2018] J. Preskill, Quantum computing in the nisq era and beyond, Quantum 2, 79 (2018).
  * Wiseman and Milburn [2009] H. M. Wiseman and G. J. Milburn,  _Quantum Measurement and Control_ (Cambridge University Press, 2009).
  * Harrington  _et al._ [2022] P. Harrington, E. Mueller, and K. Murch, Engineered dissipation for quantum information science, Nature Review Physics 4, 660 (2022).
  * Bharti  _et al._ [2022] K. Bharti, A. Cervera-Lierta, T. H. Kyaw, T. Haug, S. Alperin-Lea, A. Anand, M. Degroote, H. Heimonen, J. S. Kottmann, T. Menke, _et al._ , Noisy intermediate-scale quantum algorithms, Reviews of Modern Physics 94, 015004 (2022).
  * Cerezo  _et al._ [2021a] M. Cerezo, A. Arrasmith, R. Babbush, S. C. Benjamin, S. Endo, K. Fujii, J. R. McClean, K. Mitarai, X. Yuan, L. Cincio, _et al._ , Variational quantum algorithms, Nature Reviews Physics 3, 625 (2021a).
  * 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).
  * 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, _et al._ , The variational quantum eigensolver: a review of methods and best practices, Physics Reports 986, 1 (2022).
  * Moll  _et al._ [2018] N. Moll, P. Barkoutsos, L. S. Bishop, J. M. Chow, A. Cross, D. J. Egger, S. Filipp, A. Fuhrer, J. M. Gambetta, M. Ganzhorn, A. Kandala, A. Mezzacapo, P. Müller, W. Riess, G. Salis, J. Smolin, I. Tavernelli, and K. Temme, Quantum optimization using variational algorithms on near-term quantum devices, Quantum Science and Technology 3, 030503 (2018).
  * Bravo-Prieto  _et al._ [2019a] C. Bravo-Prieto, R. LaRose, M. Cerezo, Y. Subasi, L. Cincio, and P. J. Coles, Variational quantum linear solver (2019a).
  * 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 (2021).
  * Huang  _et al._ [2021] H.-Y. Huang, K. Bharti, and P. Rebentrost, Near-term quantum algorithms for linear systems of equations with regression loss functions, New Journal of Physics 23, 113021 (2021).
  * Li and Benjamin [2017] Y. Li and S. C. Benjamin, Efficient variational quantum simulator incorporating active error minimization, Physical Review X 7, 10.1103/physrevx.7.021050 (2017).
  * Heya  _et al._ [2019] K. Heya, K. M. Nakanishi, K. Mitarai, and K. Fujii, Subspace variational quantum simulator (2019).
  * Cîrstoiu  _et al._ [2020] C. Cîrstoiu, Z. Holmes, J. Iosue, L. Cincio, P. J. Coles, and A. Sornborger, Variational fast forwarding for quantum simulation beyond the coherence time, npj Quantum Information 6, 10.1038/s41534-020-00302-0 (2020).
  * Otten  _et al._ [2019] M. Otten, C. L. Cortes, and S. K. Gray, Noise-resilient quantum dynamics using symmetry-preserving ansatzes (2019).
  * Romero  _et al._ [2017] J. Romero, J. P. Olson, and A. Aspuru-Guzik, Quantum autoencoders for efficient compression of quantum data, Quantum Science and Technology 2, 045001 (2017).
  * Havlíček  _et al._ [2019] V. Havlíček, A. D. Córcoles, K. Temme, A. W. Harrow, A. Kandala, J. M. Chow, and J. M. Gambetta, Supervised learning with quantum-enhanced feature spaces, Nature 567, 209 (2019).
  * Benedetti  _et al._ [2019a] M. Benedetti, E. Lloyd, S. Sack, and M. Fiorentini, Parameterized quantum circuits as machine learning models, Quantum Science and Technology 4, 043001 (2019a).
  * Mangini  _et al._ [2021] S. Mangini, F. Tacchino, D. Gerace, D. Bajoni, and C. Macchiavello, Quantum computing models for artificial neural networks, Europhysics Letters 134, 10002 (2021).
  * Verdon  _et al._ [2017] G. Verdon, M. Broughton, and J. Biamonte, A quantum algorithm to train neural networks using low-depth circuits (2017).
  * Benedetti  _et al._ [2019b] M. Benedetti, D. Garcia-Pintos, O. Perdomo, V. Leyton-Ortega, Y. Nam, and A. Perdomo-Ortiz, A generative modeling approach for benchmarking and training shallow quantum circuits, npj Quantum Information 5, 10.1038/s41534-019-0157-8 (2019b).
  * Du  _et al._ [2020] Y. Du, M.-H. Hsieh, T. Liu, and D. Tao, Expressive power of parametrized quantum circuits, Physical Review Research 2, 10.1103/physrevresearch.2.033125 (2020).
  * Arrasmith  _et al._ [2019] A. Arrasmith, L. Cincio, A. T. Sornborger, W. H. Zurek, and P. J. Coles, Variational consistent histories as a hybrid algorithm for quantum foundations, Nature Communications 10, 10.1038/s41467-019-11417-0 (2019).
  * Khatri  _et al._ [2019a] S. Khatri, R. LaRose, A. Poremba, L. Cincio, A. T. Sornborger, and P. J. Coles, Quantum-assisted quantum compiling, Quantum 3, 140 (2019a).
  * Heya  _et al._ [2018] K. Heya, Y. Suzuki, Y. Nakamura, and K. Fujii, Variational quantum gate optimization (2018).
  * 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).
  * Carolan  _et al._ [2020] J. Carolan, M. Mohseni, J. P. Olson, M. Prabhu, C. Chen, D. Bunandar, M. Y. Niu, N. C. Harris, F. N. C. Wong, M. Hochberg, S. Lloyd, and D. Englund, Author correction: Variational quantum unsampling on a quantum photonic processor, Nature Physics 16, 367 (2020).
  * Johnson  _et al._ [2017] P. D. Johnson, J. Romero, J. Olson, Y. Cao, and A. Aspuru-Guzik, Qvector: an algorithm for device-tailored quantum error correction (2017).
  * 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, 1 (2018).
  * Dankert  _et al._ [2009] C. Dankert, R. Cleve, J. Emerson, and E. Livine, Exact and approximate unitary 2-designs and their application to fidelity estimation, Physical Review A 80, 10.1103/physreva.80.012304 (2009).
  * Brandão  _et al._ [2016] F. G. S. L. Brandão, A. W. Harrow, and M. Horodecki, Local random quantum circuits are approximate polynomial-designs, Communications in Mathematical Physics 346, 397 (2016).
  * Harrow and Mehraban [2018] A. W. Harrow and S. A. Mehraban, Approximate unitary $t$-designs by short random quantum circuits using nearest-neighbor and long-range gates, arXiv: Quantum Physics (2018).
  * Arrasmith  _et al._ [2021] A. Arrasmith, M. Cerezo, P. Czarnik, L. Cincio, and P. J. Coles, Effect of barren plateaus on gradient-free optimization, Quantum 5, 558 (2021).
  * Holmes  _et al._ [2022] Z. Holmes, K. Sharma, M. Cerezo, and P. J. Coles, Connecting ansatz expressibility to gradient magnitudes and barren plateaus, PRX Quantum 3, 010313 (2022).
  * Larocca  _et al._ [2022] M. Larocca, P. Czarnik, K. Sharma, G. Muraleedharan, P. J. Coles, and M. Cerezo, Diagnosing barren plateaus with tools from quantum optimal control, Quantum 6, 824 (2022).
  * Arrasmith  _et al._ [2022] A. Arrasmith, Z. Holmes, M. Cerezo, and P. J. Coles, Equivalence of quantum barren plateaus to cost concentration and narrow gorges, Quantum Science and Technology 7, 045015 (2022).
  * Wang  _et al._ [2021] 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 (2021).
  * Ortiz Marrero  _et al._ [2021] C. Ortiz Marrero, M. Kieferová, and N. Wiebe, Entanglement-induced barren plateaus, PRX Quantum 2, 040316 (2021).
  * Cerezo  _et al._ [2021b] M. Cerezo, A. Sone, T. Volkoff, L. Cincio, and P. J. Coles, Cost function dependent barren plateaus in shallow parametrized quantum circuits, Nature communications 12, 1 (2021b).
  * 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).
  * Ragone  _et al._ [2023] M. Ragone, B. N. Bakalov, F. Sauvage, A. F. Kemper, C. O. Marrero, M. Larocca, and M. Cerezo, A unified theory of barren plateaus for deep parametrized quantum circuits, arXiv 10.48550/arxiv.2309.09342 (2023).
  * Fontana  _et al._ [2023] E. Fontana, D. Herman, S. Chakrabarti, N. Kumar, R. Yalovetzky, J. Heredge, S. H. Sureshbabu, and M. Pistoia, The adjoint is all you need: Characterizing barren plateaus in quantum ansätze, arXiv 10.48550/arxiv.2309.07902 (2023).
  * Grant  _et al._ [2019] E. Grant, L. Wossnig, M. Ostaszewski, and M. Benedetti, An initialization strategy for addressing barren plateaus in parametrized quantum circuits, Quantum 3, 214 (2019).
  * Friedrich and Maziero [2022] L. Friedrich and J. Maziero, Avoiding barren plateaus with classical deep neural networks, Phys. Rev. A 106, 042433 (2022).
  * Mele  _et al._ [2022] A. A. Mele, G. B. Mbeng, G. E. Santoro, M. Collura, and P. Torta, Avoiding barren plateaus via transferability of smooth solutions in a hamiltonian variational ansatz, Phys. Rev. A 106, L060401 (2022).
  * Patti  _et al._ [2021] T. L. Patti, K. Najafi, X. Gao, and S. F. Yelin, Entanglement devised barren plateau mitigation, Phys. Rev. Res. 3, 033090 (2021).
  * Volkoff and Coles [2021] T. Volkoff and P. J. Coles, Large gradients via correlation in random parameterized quantum circuits, Quantum Science and Technology 6, 025008 (2021).
  * Sack  _et al._ [2022] S. H. Sack, R. A. Medina, A. A. Michailidis, R. Kueng, and M. Serbyn, Avoiding barren plateaus using classical shadows, PRX Quantum 3, 020365 (2022).
  * Dborin  _et al._ [2022] J. Dborin, F. Barratt, V. Wimalaweera, L. Wright, and A. G. Green, Matrix product state pre-training for quantum machine learning, Quantum Science and Technology 7, 035014 (2022).
  * Kieferova  _et al._ [2021] M. Kieferova, O. M. Carlos, and N. Wiebe, Quantum generative training using rényi divergences, arXiv 10.48550/ARXIV.2106.09567 (2021).
  * Skolik  _et al._ [2021] A. Skolik, J. R. McClean, M. Mohseni, P. van der Smagt, and M. Leib, Layerwise learning for quantum neural networks, Quantum Machine Intelligence 3, 10.1007/s42484-020-00036-4 (2021).
  * Tacchino  _et al._ [2021] F. Tacchino, S. Mangini, P. K. Barkoutsos, C. Macchiavello, D. Gerace, I. Tavernelli, and D. Bajoni, Variational learning for quantum artificial neural networks, IEEE Transactions on Quantum Engineering 2, 1 (2021).
  * Mazzola  _et al._ [2019] G. Mazzola, P. J. Ollitrault, P. K. Barkoutsos, and I. Tavernelli, Nonunitary operations for ground-state calculations in near-term quantum computers, Physical Review Letters 123, 10.1103/physrevlett.123.130501 (2019).
  * Benfenati  _et al._ [2021] F. Benfenati, G. Mazzola, C. Capecci, P. K. Barkoutsos, P. J. Ollitrault, I. Tavernelli, and L. Guidoni, Improved accuracy on noisy devices by nonunitary variational quantum eigensolver for chemistry applications, Journal of Chemical Theory and Computation 17, 3946 (2021).
  * Melo  _et al._ [2023] A. Melo, N. Earnest-Noble, and F. Tacchino, Pulse-efficient quantum machine learning, Quantum 7, 1130 (2023).
  * Beer  _et al._ [2020] K. Beer, D. Bondarenko, T. Farrelly, T. J. Osborne, R. Salzmann, D. Scheiermann, and R. Wolf, Training deep quantum neural networks, Nature Communications 11, 10.1038/s41467-020-14454-2 (2020).
  * Bondarenko and Feldmann [2020] D. Bondarenko and P. Feldmann, Quantum autoencoders to denoise quantum data, Phys. Rev. Lett. 124, 130502 (2020).
  * Poland  _et al._ [2020] K. Poland, K. Beer, and T. J. Osborne, No free lunch for quantum machine learning (2020).
  * Verstraete  _et al._ [2009] F. Verstraete, M. M. Wolf, and J. Ignacio Cirac, Quantum computation and quantum-state engineering driven by dissipation, Nature physics 5, 633 (2009).
  * Pastawski  _et al._ [2011] F. Pastawski, L. Clemente, and J. I. Cirac, Quantum memories based on engineered dissipation, Physical Review A 83, 012304 (2011).
  * Breuer  _et al._ [2002] H.-P. Breuer, F. Petruccione, _et al._ ,  _The theory of open quantum systems_ (Oxford University Press on Demand, 2002).
  * Gorini  _et al._ [1976] V. Gorini, A. Kossakowski, and E. C. G. Sudarshan, Completely positive dynamical semigroups of n-level systems, Journal of Mathematical Physics 17, 821 (1976).
  * Lindblad [1976] G. Lindblad, On the generators of quantum dynamical semigroups, Communications in Mathematical Physics 48, 119 (1976).
  * Khatri  _et al._ [2019b] S. Khatri, R. LaRose, A. Poremba, L. Cincio, A. T. Sornborger, and P. J. Coles, Quantum-assisted quantum compiling, Quantum 3, 140 (2019b).
  * LaRose  _et al._ [2019] R. LaRose, A. Tikku, É. O’Neel-Judy, L. Cincio, and P. J. Coles, Variational quantum state diagonalization, npj Quantum Information 5, 10.1038/s41534-019-0167-6 (2019).
  * Bravo-Prieto  _et al._ [2019b] C. Bravo-Prieto, R. LaRose, M. Cerezo, Y. Subasi, L. Cincio, and P. J. Coles, Variational quantum linear solver (2019b).
  * Letcher  _et al._ [2023] A. Letcher, S. Woerner, and C. Zoufal, From tight gradient bounds for parameterized quantum circuits to the absence of barren plateaus in qgans (2023), arXiv:2309.12681 [quant-ph] .
  * Minganti  _et al._ [2018] F. Minganti, A. Biella, N. Bartolo, and C. Ciuti, Spectral theory of liouvillians for dissipative phase transitions, Phys. Rev. A 98, 042118 (2018).
  * Note [1] The following calculations can be immediately generalized to cases where each $\EuScript{L}_{k}$ operator acts on a generic number of qubits less than or equal to $K$.
  * Note [2] We have omitted the $\bm{\sigma}_{q}$ dependence of all these terms to make the notation easier.
  * Haga  _et al._ [2021] T. Haga, M. Nakagawa, R. Hamazaki, and M. Ueda, Liouvillian skin effect: Slowing down of relaxation processes without gap closing, Phys. Rev. Lett. 127, 070402 (2021).
  * Kastoryano  _et al._ [2012] M. J. Kastoryano, D. Reeb, and M. M. Wolf, A cutoff phenomenon for quantum markov chains, Journal of Physics A: Mathematical and Theoretical 45, 075307 (2012).
  * Note [3] It is a direct consequence of the following Liovillian property: if $\rho$ is a left eigenvector of $\EuScript{L}$ (or a right one of $\EuScript{L}^{\dagger}$) with eigenvalue equal to $\lambda$ then $\rho^{\dagger}$ will be in turn an eigenvector with $\lambda^{*}$ as an eigenvalue. This statement can be proved easily by applying the proof of Lemma 3 in [68] considering the adjoint of a generic Liouvillian.
  * Jordan and Wigner [1928] P. Jordan and E. Wigner, Über das paulische Äquivalenzverbot, Zeitschrift für Physik 47, 631 (1928).
  * Hehre  _et al._ [1969] W. J. Hehre, R. F. Stewart, and J. A. Pople, Self-consistent molecular-orbital methods. i. use of gaussian expansions of slater-type atomic orbitals, The Journal of Chemical Physics 51, 2657 (1969).
  * Slater [1951] J. C. Slater, A simplification of the hartree-fock method, Phys. Rev. 81, 385 (1951).
  * Cichy  _et al._ [2022] S. Cichy, P. K. Faehrmann, S. Khatri, and J. Eisert, A perturbative gadget for delaying the onset of barren plateaus in variational quantum algorithms, arXiv 10.48550/ARXIV.2210.03099 (2022).
  * Cattaneo  _et al._ [2023] M. Cattaneo, M. A. Rossi, G. García-Pérez, R. Zambrini, and S. Maniscalco, Quantum simulation of dissipative collective effects on noisy quantum computers, PRX Quantum 4, 10.1103/prxquantum.4.010324 (2023).
  * Chen  _et al._ [2023] C.-F. Chen, H.-Y. Huang, J. Preskill, and L. Zhou, Local minima in quantum systems (2023).
  * Mujal  _et al._ [2021] P. Mujal, R. Martínez-Peña, J. Nokkala, J. García-Beni, G. L. Giorgi, M. C. Soriano, and R. Zambrini, Opportunities in quantum reservoir computing and extreme learning machines, Advanced Quantum Technologies 4, 2100027 (2021).
  * Kubota  _et al._ [2023] T. Kubota, Y. Suzuki, S. Kobayashi, Q. H. Tran, N. Yamamoto, and K. Nakajima, Temporal information processing induced by quantum noise, Phys. Rev. Res. 5, 023057 (2023).
  * Sannia  _et al._ [2022] A. Sannia, R. Martínez-Peña, M. C. Soriano, G. L. Giorgi, and R. Zambrini, Dissipation as a resource for quantum reservoir computing, arXiv (2022), 2212.12078 [quant-ph] .
  * Domingo  _et al._ [2023] L. Domingo, G. Carlo, and F. Borondo, Taking advantage of noise in quantum reservoir computing, Scientific Reports 13, 8790 (2023).
  * Liu  _et al._ [2022] J. Liu, F. Wilde, A. A. Mele, L. Jiang, and J. Eisert, Noise can be helpful for variational quantum algorithms, arXiv preprint arXiv:2210.06723 (2022).
  * Miessen  _et al._ [2023] A. Miessen, P. J. Ollitrault, F. Tacchino, and I. Tavernelli, Quantum algorithms for quantum dynamics, Nature Computational Science 3, 25 (2023).
  * Kliesch  _et al._ [2011] M. Kliesch, T. Barthel, C. Gogolin, M. Kastoryano, and J. Eisert, Dissipative quantum church-turing theorem, Phys. Rev. Lett. 107, 120501 (2011).
  * 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).
  * Kamakari  _et al._ [2022] H. Kamakari, S.-N. Sun, M. Motta, and A. J. Minnich, Digital quantum simulation of open quantum systems using quantum imaginary–time evolution, PRX Quantum 3, 010320 (2022).
  * Ramusat and Savona [2021] N. Ramusat and V. Savona, A quantum algorithm for the direct estimation of the steady state of open quantum systems, Quantum 5, 399 (2021).
  * Cleve and Wang [2017] R. Cleve and C. Wang, Efficient Quantum Algorithms for Simulating Lindblad Evolution, in  _44th International Colloquium on Automata, Languages, and Programming (ICALP 2017)_, Leibniz International Proceedings in Informatics (LIPIcs), Vol. 80, edited by I. Chatzigiannakis, P. Indyk, F. Kuhn, and A. Muscholl (Schloss Dagstuhl–Leibniz-Zentrum fuer Informatik, Dagstuhl, Germany, 2017) pp. 17:1–17:14.
  * Head-Marsden  _et al._ [2021] K. Head-Marsden, S. Krastanov, D. A. Mazziotti, and P. Narang, Capturing non-markovian dynamics on near-term quantum computers, Phys. Rev. Research 3, 013182 (2021).
  * Cattaneo  _et al._ [2021] M. Cattaneo, G. D. Chiara, S. Maniscalco, R. Zambrini, and G. L. Giorgi, Collision models can efficiently simulate any multipartite markovian quantum dynamics, Physical Review Letters 126, 10.1103/physrevlett.126.130403 (2021).
  * Murch  _et al._ [2012] K. W. Murch, U. Vool, D. Zhou, S. J. Weber, S. M. Girvin, and I. Siddiqi, Cavity-assisted quantum bath engineering, Phys. Rev. Lett. 109, 183602 (2012).
  * Rost  _et al._ [2021] B. Rost, L. D. 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 preprint (2021), arXiv:2108.01183 .
  * Leppäkangas  _et al._ [2022] J. Leppäkangas, N. Vogt, K. R. Fratus, K. Bark, J. A. Vaitkus, P. Stadler, J.-M. Reiner, S. Zanker, and M. Marthaler, A quantum algorithm for solving open quantum system dynamics on quantum computers using noise, arXiv preprint arXiv:2210.12138 (2022).
  * Chenu  _et al._ [2017] A. Chenu, M. Beau, J. Cao, and A. del Campo, Quantum simulation of generic many-body open system dynamics using classical noise, Phys. Rev. Lett. 118, 140403 (2017).
  * 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 , 1 (2023).
  * Potočnik  _et al._ [2018] A. Potočnik, A. Bargerbos, F. A. Y. N. Schröder, S. A. Khan, M. C. Collodo, S. Gasparinetti, Y. Salathé, C. Creatore, C. Eichler, H. E. Türeci, A. W. Chin, and A. Wallraff, Studying light-harvesting models with superconducting circuits, Nature Communications 9, 10.1038/s41467-018-03312-x (2018).
  * Ciccarello  _et al._ [2022] F. Ciccarello, S. Lorenzo, V. Giovannetti, and G. M. Palma, Quantum collision models: Open system dynamics from repeated interactions, Physics Reports 954, 1 (2022).
  * Qiskit contributors [2023] Qiskit contributors, Qiskit: An open-source framework for quantum computing (2023).