##  1 Overview

Polynomial identity testing (PIT) is the fundamental problem of deciding
whether a given multi-variate algebraic circuit formally computes the zero
polynomial. PIT has a simple efficient randomized algorithm that only needs
blackbox access to the circuit: Pick a random point and check whether the
circuit evaluates to zero on that particular point.

Despite the fundamental nature of PIT and the simplicity of the randomized
algorithm, no efficient deterministic algorithm is known—even in the white-box
setting, where the algorithm has access to the description of the circuit. The
existence of such an algorithm would imply long-sought circuit lower bounds
[HeintzSchnorr1980, Agrawal2005, KabanetsImpagliazzo2004]. Conversely,
sufficiently strong circuit lower bounds yield blackbox derandomization for
all of BPP, the class of decision problems admitting efficient randomized
algorithms with bounded error [NisanWigderson1994, ImpagliazzoWigderson1997].
Although the known results leave gaps between the two directions, they suggest
that PIT acts as a BPP-complete problem in the context of derandomization, and
that derandomization of BPP can be achieved in a blackbox fashion if at all.

Blackbox derandomization of PIT for a class of polynomials $\Class$ in the
variables $x_{1}$, …, $x_{n}$ is equivalent to the efficient construction of a
substitution $G$ that replaces each $x_{i}$ by a low-degree polynomial in a
small set of fresh variables such that, for every nonzero polynomial $p$ from
$\Class$, $p(G)$ remains nonzero [ShpilkaYehudayoff2010, Lemma 4.1]. We refer
to $G$ as a generator, the fresh variables are its seed, and say that $G$ hits
the class $\Class$. If there are $l$ seed variables, and if $p$ and $G$ have
degree at most $n^{O(1)}$, then the resulting deterministic PIT algorithm for
$\Class$ makes $n^{O(l)}$ blackbox queries.

Much progress on derandomizing PIT has been obtained by designing such
substitutions and analyzing their hitting properties for interesting classes
$\Class$. Shpilka and Volkovich [ShpilkaVolkovich2008] introduced a generator,
by now dubbed the Shpilka–Volkovich generator or “SV generator” for short, and
proved that it hits sums of a bounded number of read-once formulas for
$l=O(\log n)$, later improved to $l=O(1)$ [MinahanVolkovich2017]. The
generator for $l=O(\log n)$ has also been shown to hit multilinear depth-4
circuits with bounded top fan-in [KMSV2009], multilinear bounded-read formulas
[AvMV2011], commutative read-once oblivious algebraic branching programs
[FSS2014], $\Sigma\mathrm{m}\bigwedge\Sigma\Pi^{O(1)}$ formulas [Forbes2015],
circuits with locally-low algebraic rank in the sense of [KumarSaraf2016], and
orbits of simple polynomial classes under invertible linear transformations of
the variables [MediniShpilka2021]. The generator is an ingredient in other
hitting set constructions, as well, notably constructions using the technique
of low-support rank concentration [ASS2013, AGKS2014, GKST2015, GKS2016,
SahaThankey2021, BhargavaGhosh2021]. It also forms the core of a “succinct”
generator that hits a variety of classes including depth-2 circuits [FSV2017].

### Vanishing ideal.

In this paper, we initiate a systematic study of the power of a generator $G$
through the set of polynomials $p$ such that $p(G)$ vanishes. The set has the
algebraic structure of an ideal and is known as the _vanishing ideal_ of $G$;
we denote it by $\VanIdeal{G}$. Our technical contributions can be understood
as precisely characterizing the vanishing ideal of the SV generator.

Characterizing the vanishing ideal facilitates two objectives:

Derandomization.

A generator $G$ hits a class $\Class$ of polynomials if and only if $\Class$
and $\VanIdeal{G}$ have at most the zero polynomial in common. If the
characterization of $\VanIdeal{G}$ is incompatible with being computable
within some resource bound, then $G$ trivially hits the class $\Class$ of
polynomials that are computable within the bound. In other words,
derandomization of PIT for $\Class$ reduces to proving lower bounds for
$\VanIdeal{G}$. Characterizing $\VanIdeal{G}$ yields explicit structure that
makes the lower bounds more tractable.

More generally, given a characterization of $\VanIdeal{G}$, in order to
derandomize PIT for a class $\Class$ it suffices to design another generator
$G^{\prime}$ that hits merely the polynomials in $\Class\cap\VanIdeal{G}$. As
$G$ hits the remainder of $\Class$, combining $G$ with $G^{\prime}$ yields a
generator for all of $\Class$. In this way, one may assume for free additional
structure about the polynomials in $\Class$, namely that the polynomials
moreover belong to $\VanIdeal{G}$. Characterizing the vanishing ideal of $G$
makes this additional structure explicit.

Lower bounds.

If we happen to know that $G$ hits the class $\Class$ of polynomials
computable within some resource bound, then any expression for a nonzero
polynomial in $\VanIdeal{G}$ yields an explicit polynomial that falls outside
$\Class$. Such a statement is often referred to as hardness of representation,
and can be viewed as a lower bound in the model of computation underlying
$\Class$ (provided the polynomial can be computed in the model at all).
Characterizing $\VanIdeal{G}$ makes explicit the polynomials to which the
lower bound applies.

We will illustrate how to make progress on both objectives through our
characterizations of the SV generator’s vanishing ideal.

### Rational function evaluations.

Another contribution of our paper is the development of an alternate view of
the SV generator, namely as evaluations of univariate rational functions of
low degree. We would like to promote the perspective for its intrinsic appeal
and its applicability. Among other benefits, it facilitates the study of the
vanishing ideal.

The transition goes as follows. The SV generator takes as additional
parameters a positive integer $l$ and a choice of distinct field elements
$a_{i}$ for each of the original variables $x_{i}$, $i\in[n]$. We refer to the
elements $a_{i}$ as _abscissas_ and denote the generator for a given value of
$l$ by $\SVl$ (suppressing the choice of abscissas). When $l=1$, $\SVparam{1}$
takes as seed two fresh variables, $y$ and $z$, and can be described
succinctly in terms of the Lagrange interpolants $L_{i}$ for the set of
abscissas:

$x_{i}\leftarrow z\cdot L_{i}(y)\doteq
z\cdot\prod_{j\in[n]\setminus\\{i\\}}\frac{y-a_{j}}{a_{i}-a_{j}}.$ (1)

By rescaling, the denominators on the right-hand side of eq.intro.SV can be
cleared, resulting in the following somewhat simpler substitution:

$x_{i}\leftarrow z\cdot\prod_{j\in[n]\setminus\\{i\\}}(y-a_{j}).$ (2)

The vanishing ideals of eq.intro.SVnonnorm and $\SVparam{1}$ are the same up
to rescaling the variables to match the rescaling from eq.intro.SV to
eq.intro.SVnonnorm.

More importantly, we apply the change of variables $z\leftarrow
z^{\prime}/\prod_{j\in[n]}(y-a_{j})$, resulting in a substitution that now
uses rational functions of the seed:

$x_{i}\leftarrow\frac{z^{\prime}}{y-a_{i}}.$ (3)

The notion of vanishing ideal naturally extends to rational function
substitutions. The change of variables from eq.intro.SVnonnorm to eq.intro.RFE
establishes that any polynomial vanishing on eq.intro.SVnonnorm also vanishes
on eq.intro.RFE. The change of variables is invertible (the inverse is
$z^{\prime}\leftarrow z\prod_{j\in[n]}(y-a_{j})$), so any polynomial vanishing
on eq.intro.RFE also vanishes on eq.intro.SVnonnorm. We conclude that the
vanishing ideal of eq.intro.RFE is the same as that of $\SVparam{1}$ up to
rescaling the variables.

Note that, for fixed $y$ and $z^{\prime}$, eq.intro.RFE may be interpreted as
first forming a univariate rational function
$f(\alpha)=\frac{z^{\prime}}{y-\alpha}$ (depending on $y$ and $z^{\prime}$ but
independent of $i$) and then substituting $x_{i}\leftarrow f(a_{i})$. As $y$
and $z^{\prime}$ vary, $f$ ranges over all rational functions with numerator
degree zero and denominator degree one. We denote eq.intro.RFE by
$\RFEparam{0}{1}$, where $\RFE$ is a short-hand for Rational Function
Evaluation, $0$ bounds the numerator degree, and $1$ bounds the denominator
degree.

As a generator, $\RFEparam{0}{1}$ naturally generalizes to $\RFEkl$ for
arbitrary $k,l\in\NN$.

######  Definition 1.1 (RFE Generator).

Let $\FF$ be a field and $\\{x_{1},\ldots,x_{n}\\}$ a set of variables. The
_Rational Function Evaluation Generator (RFE)_ for $\FF[x_{1},\ldots,x_{n}]$
is parametrized by the following data:

  * •

For each $i\in[n]$, a distinct _abscissa_ $a_{i}\in\FF$.

  * •

A non-negative integer $k$, the _numerator degree_.

  * •

A non-negative integer $l$, the _denominator degree_.

The generator takes as seed a rational function $f\in\FF(\alpha)$ such that
$f$ can be written as $g/h$ for some $g,h\in\FF[\alpha]$ with $\deg(g)\leq k$,
$\deg(h)\leq l$, and $h(a_{i})\neq 0$ for all $i\in[n]$. From $f$, it
generates the substitution $x_{i}\leftarrow f(a_{i})$ for each $i\in[n]$.

There are multiple ways to parametrize the seed of $\RFEkl$ using scalars,
such as by specifying the coefficients, evaluations, or roots for each of the
numerator and denominator. The flexibility to choose is a source of
convenience. As is customary in the context of blackbox derandomization of
PIT, we assume that $\FF$ is sufficiently large, possibly by taking a field
extension. We refer to a.rfe-as-hsg for a discussion on different
parametrizations as well as on how to obtain deterministic blackbox PIT
algorithms from the generator and how large the underlying field $\FF$ must
be.

The connection between $\RFEparam{0}{1}$ and $\SVparam{1}$ extends as follows.
For higher values of $l$, $\SVl$ is defined as the sum of $l$ independent
instantiations of $\SVparam{1}$. The same transformations as above relate
$\SVl$ and the sum of $l$ independent instantiations of $\RFEparam{0}{1}$. The
latter in turn is equivalent to $\RFElml$ by partial fraction decomposition.
The conclusion is that $\SVparam{l}$ is equivalent in power to $\RFElml$, up
to variable rescaling. We refer to a.rfe-sv for a formal treatment.

For parameter values $k\neq l-1$, there is no SV generator that corresponds to
$\RFEkl$, but $\SVparam{\max(k+1,l)}$ encompasses $\RFEkl$ (up to rescaling)
and uses at most twice as many seed variables. Thus, the $\RFE$-generator and
the $\SV$-generator efficiently hit the same classes of polynomials. However,
$\RFE$’s simpler univariate dependence on the abscissas—as opposed to $\SV$’s
multi-variate dependence—enables our approach for determining the vanishing
ideal. The moral is that, even though polynomial substitutions are sufficient
for derandomizing PIT, it nevertheless helps to consider rational
substitutions. Their use may simplify analysis and arguably yield more elegant
constructions.

### Generating set.

Our first result describes a small and explicit generating set for the
vanishing ideal of $\RFE$. It consists of instantiations of a single
determinant expression.

######  Theorem 1.2 (generating set).

Let $k,l\in\NN$, $\\{x_{i}:i\in[n]\\}$ be a set of variables, and let $a_{i}$
for $i\in[n]$ be distinct field elements. The vanishing ideal of $\RFEkl$ over
the given set of variables for the given choice of abscissas
$(a_{i})_{i\in[n]}$ is generated by the following polynomials over all choices
of $k+l+2$ variable indices $i_{1},i_{2},\dots,i_{k+l+2}\in[n]$:

$\VanGenkl[i_{1},i_{2},\dots,i_{k+l+2}]\doteq\det\bmatrix
a_{i_{j}}^{l}x_{i_{j}}&a_{i_{j}}^{l-1}x_{i_{j}}\dots
x_{i_{j}}a_{i_{j}}^{k}a_{i_{j}}^{k-1}\dots 1_{j=1}^{k+l+2}.$ (4)

Moreover, for any fixed set $C\subseteq[n]$ of $k+1$ variable indices, when
$i_{1},i_{2},\ldots,i_{k+l+2}$ range over all the choices of $k+l+2$ indices
from $[n]$ such that $i_{1}<\cdots<i_{k+l+2}$ and
$C\subseteq\\{i_{1},\ldots,i_{k+l+2}\\}$, the polynomials
$\VanGenkl[i_{1},i_{2},\dots,i_{k+l+2}]$ form a generating set of minimum
size.

The name “$\VanGen$” is a shorthand for “Elementary Vandermonde Circulation”.
Later we discuss a representation of polynomials using alternating algebra,
with connections to notions from network flow. In this representation,
polynomials in the vanishing ideal coincide with circulations, and
instantiations of $\VanGen$ are the elementary circulations.

We refer to the set $C$ in thm.ideal.generators as a _core_. The core $C$
plays a similar role as in a combinatorial sunflower except that, unlike the
petals of a sunflower, the sets $\\{i_{1},\ldots,i_{k+l+2}\\}$ do not need to
be disjoint outside the core.

######  Example 1.3.

Consider the special case where $k=0$ and $l=1$. The generator for
$\RFEparam{0}{1}$ when $i_{1}=1$, $i_{2}=2$, and $i_{3}=3$ is given by

$\VanGenparam{0}{1}[1,2,3]\doteq\vmatrix*a_{1}x_{1}&x_{1}1\\\
a_{2}x_{2}x_{2}1\\\
a_{3}x_{3}x_{3}1=(a_{1}-a_{2})x_{1}x_{2}+(a_{2}-a_{3})x_{2}x_{3}+(a_{3}-a_{1})x_{3}x_{1}.$

For any fixed $i^{*}\in[n]$, the polynomials
$\VanGenparam{0}{1}[i_{1},i_{2},i_{3}]$ form a generating set of minimum size
when $\\{i_{1},i_{2},i_{3}\\}$ ranges over all subsets of $[n]$ that contain
$C=\\{i^{*}\\}$, and $i_{1}<i_{2}<i_{3}$.

In general, the generators $\VanGenkl$ are nonzero multilinear homogeneous
polynomials of degree $l+1$ containing all multilinear monomials of degree
$l+1$.

Each generating set of minimum size in thm.ideal.generators yields a Gröbner
basis with respect to every monomial order that prioritizes the variables
outside $C$. A Gröbner basis is a special basis that allows solving ideal-
membership queries more efficiently as well as solving systems of polynomial
equations [CLO2013]. Computing Gröbner bases for general ideals is
exponential-space complete [KuhnleM1996, Mayr1997]. thm.ideal.generators
represents a rare instance of a natural and interesting ideal for which we
know a small and moreover explicit Gröbner basis.

To gain some intuition about dependencies between the generators $\VanGenkl$,
note that permuting the order of the variables used in the construction of
$\VanGenkl$ yields the same polynomial or minus that polynomial, depending on
the sign of the permutation. This follows from the determinant structure of
$\VanGenkl$ and is the reason why we need to fix the order of the variables in
order to obtain a generating set of minimum size. More profoundly, the
following relationship holds for every choice of $k+l+3$ indices
$i_{1},i_{2},\dots,i_{k+l+3}\in[n]$ and every univariate polynomial $w$ of
degree at most $k$:

$\det\bmatrix w(a_{i_{j}})&a_{i_{j}}^{l}x_{i_{j}}a_{i_{j}}^{l-1}x_{i_{j}}\dots
x_{i_{j}}a_{i_{j}}^{k}a_{i_{j}}^{k-1}\dots 1_{j=1}^{k+l+3}=0.$ (5)

The determinant in eq.generator.dependencies vanishes because the first column
of the matrix is a linear combination of the last $k+1$. A Laplace expansion
across the first column allows us to write the determinant of the matrix as a
linear combination of minors, and each minor is an instantiation of
$\VanGenkl$. Since the determinant vanishes, eq.generator.dependencies
represents a linear dependency for every nonzero polynomial $w$ of degree at
most $k$. In fact, when $\\{i_{1},\ldots,i_{k+l+3}\\}$ varies over subsets of
$[n]$ containing a fixed core of size $k+1$, the equations
eq.generator.dependencies generate _all_ linear dependencies among instances
of $\VanGenkl$.

As corollaries to thm.ideal.generators we obtain the following tight bounds on
$\VanIdeal{\RFEkl}$. The bounds hold for every way to choose the parameters in
def.rfegen.

  * •

The minimum degree of a nonzero polynomial in $\VanIdeal{\RFEkl}$ equals
$l+1$. This proves a conjecture by Fournier and Korwar [FournierKorwar2018]
(additional partial results reported in [Korwar2021]) that there exists a
polynomial of degree $l+1$ in $n=2l+1$ variables that $\SVl$ fails to hit. The
conjecture follows because the generators for $\VanIdeal{\SVl}$ have degree
$l+1$ and use $2l+1$ variables.

As none of the generators contain a monomial of support $l$ or less, the same
holds for every nonzero polynomial in $\VanIdeal{\RFEkl}$. This extends the
known property that $\SVl$ hits every polynomial that contains a monomial of
support $l$ or less [ShpilkaVolkovich2008].

  * •

The minimum sparseness, i.e., number of monomials, of a nonzero polynomial in
$\VanIdeal{\RFEkl}$ equals $\binom{k+l+2}{l+1}$. The generators $\VanGenkl$
realize the bound as they exactly contain all multilinear monomials of degree
$l+1$ that can be formed out of their $k+l+2$ variables.

The claim that no nonzero polynomial in $\VanIdeal{\RFEkl}$ contains fewer
than $\binom{k+l+2}{l+1}$ monomials requires an additional combinatorial
argument. It is a (tight) quantitative strengthening of the known property
that $\SVl$ hits every polynomial with fewer than $2^{l}$ monomials [AvMV2011,
GKST2015, Forbes2015, FSV2017]. Note that for $k=l-1$ we have that
$\binom{k+l+2}{l+1}=\binom{2l+1}{l+1}=\Theta(2^{2l}/\sqrt{l})$.

  * •

The minimum partition class size of a nonzero set-multilinear polynomial of
degree $l+1$ in $\VanIdeal{\RFEkl}$ equals $k+2$. Set-multilinearity is a
common restriction in works on derandomizing PIT and algebraic circuit lower
bounds. A polynomial $p$ of degree $l+1$ in a set of variables
$\\{x_{1},\ldots,x_{n}\\}$ is said to be set-multilinear if $[n]$ can be
partitioned as $[n]=X_{1}\sqcup X_{2}\sqcup\dots\sqcup X_{l+1}$ such that
every monomial in $p$ is a product $x_{i_{1}}\cdot x_{i_{2}}\cdot\dots\cdot
x_{i_{l+1}}$, where $i_{j}\in X_{j}$. Note that set-multilinearity implies
multilinearity but not the other way around.

As the generators $\VanGenkl$ are not set-multilinear, it is not immediately
clear from thm.ideal.generators whether $\VanIdeal{\RFEkl}$ contains
nontrivial set-multilinear polynomials. However, a variation on the
construction of the generators $\VanGenkl$ yields explicit set-multilinear
homogeneous polynomials in $\VanIdeal{\RFEkl}$ of degree $l+1$ where each
$X_{j}$ has size $k+2$. We denote them by $\SMVanOnekl$, where $\SMVanOne$
stands for “Elementary Set-Multilinear Vandermonde Circulation”. $\SMVanOnekl$
contains all monomials of the form $x_{i_{1}}\cdot x_{i_{2}}\cdot\dots\cdot
x_{i_{l+1}}$ with $i_{j}\in X_{j}$. For any variable partition $X_{1}\sqcup
X_{2}\sqcup\dots\sqcup X_{l+1}$ with $|X_{1}|=\cdots=|X_{l+1}|=k+2$,
$\SMVanOnekl$ is the only set-multilinear polynomial in $\VanIdeal{\RFEkl}$
with that variable partition, up to a scalar multiple.

### Membership test.

Our second characterization of the vanishing ideal of $\RFE$ can be viewed as
a structured membership test. Given a polynomial $p$, there is a generic way
to test whether $p$ belongs to the vanishing ideal of any generator $G$,
namely by symbolically substituting $G$ into $p$ and verifying that the result
simplifies to zero. When $G$ is a polynomial substitution, the well-known
transformation of a generator into a deterministic blackbox PIT algorithm
yields another test: verify $p(G)=0$ for a sufficiently large set of
substitutions into the seed variables [Ore1922, DeMilloLipton1978, Zippel1979,
Schwartz1980]. By clearing denominators, the same goes for rational
substitutions like $\RFEkl$. While these tests work, their genericity with
respect to $G$ implies that they cannot provide any $G$-specific insight into
whether or why a given $p$ belongs to the vanishing ideal of $G$. As we will
argue, our structured membership test does.

Several prior papers demonstrated the utility of partial derivatives and zero
substitutions in the context of derandomizing PIT using the SV generator,
especially for syntactically multilinear models [ShpilkaVolkovich2008,
KMSV2009, AvMV2011]. Building on the generating set of thm.ideal.generators,
we state a more structured test for membership in $\VanIdeal{\RFEkl}$ in terms
of those operations.

######  Theorem 1.4 (membership test).

Let $k,l\in\NN$, $X=\\{x_{i}:i\in[n]\\}$ be a set of variables, $a_{i}$ for
$i\in[n]$ be distinct field elements, and $Z$ be a set of at least $n-k-l-1$
nonzero field elements. A multilinear polynomial $p\in\FF[X]$ belongs to
$\VanIdeal{\RFEkl}$ if and only if both of the following conditions hold:

  1. 1.

condition]thm.ideal.membership.cond1 There are no monomials of degree $l$ or
less, nor of degree $n-k$ or more, in $p$.

  2. 2.

condition]thm.ideal.membership.cond2 For all disjoint subsets
$K,L\subseteq[n]$ with $|K|=k$ and $|L|=l$, and every $z\in Z$,
$\PartialZero{L}{K}{p}$ evaluates to zero upon the following substitution for
each $i\in\overline{K\cup L}$

$x_{i}\;\leftarrow\;z\cdot\frac{\prod_{j\in K}(a_{i}-a_{j})}{\prod_{j\in
L}(a_{i}-a_{j})}.$ (6)

A few technical comments regarding the statement: The first part of
thm.ideal.membership.cond1 in thm.ideal.membership reflects the known property
that $\SVl$ hits every multilinear polynomial that contains a monomial of
degree $l$ or less [ShpilkaVolkovich2008]. The second part expresses an
analogous property, and together they imply that all multilinear polynomials
on $n\leq k+l+1$ variables are hit by $\RFEkl$. In thm.ideal.membership.cond2,
$\PartialZero{L}{K}{p}$ denotes the polynomial obtained by taking the partial
derivative of $p$ with respect to every variable in $L$ and setting all the
variables in $K$ to zero. The order of the operations does not matter, and the
resulting polynomial depends only on variables in $\overline{K\cup L}$. The
set $Z$ serves only to separate the homogeneous components of $p$ (cf.
thm.ideal.membership.homogcase and fact.rfe-homog in s.test). One can just as
well use a fresh variable $z$ instead of using substitions from $Z$, in which
case thm.ideal.membership.cond2 becomes a single substitution for each $K$ and
$L$. The substitution eq.thm.ideal.membership.sub can be viewed as
$x_{i}\leftarrow z\cdot f_{K,L}(a_{i})$, where

$f_{K,L}(\alpha)\doteq\frac{\prod_{j\in K}(\alpha-a_{j})}{\prod_{j\in
L}(\alpha-a_{j})}$

is a valid seed for $\RFEkl$ when substituted into the variables in
$\overline{K\cup L}$.

All together, thm.ideal.membership can be understood as stating that a
multilinear polynomial $p$ is hit by $\RFEkl$ if and only if $p$ has a
monomial supported on few or all-but-few variables, or else there is a small
set of zero substitutions, $K$, and a small set of partial derivatives, $L$,
whose application to $p$ leaves a polynomial that is nonzero after
substituting $x_{i}\leftarrow z\cdot f_{K,L}(a_{i})$. As we mentioned above,
several prior papers demonstrated the utility of partial derivatives and zero
substitutions in the context of derandomizing PIT using the SV generator. By
judiciously choosing variables for those operations, these papers managed to
simplify $p$ and reduce PIT for $p$ to PIT for simpler instances, resulting in
efficient recursive algorithms. In s.test, we develop a general framework for
such algorithms and prove correctness directly from thm.ideal.membership.
Moreover, because thm.ideal.membership is a precise characterization, any
argument that $\SV$ or $\RFE$ hits a class of multilinear polynomials can be
converted into one within our framework, i.e., into an argument based on zero
substitutions and partial derivatives. Thus, thm.ideal.membership shows that
these tools harness the complete power of $\SV$ and $\RFE$ for multilinear
polynomials.

### Applications.

We illustrate the utility of our characterizations of the vanishing ideal of
$\RFE$ in the two directions mentioned before.

#### Derandomization.

To start, we demonstrate how thm.ideal.membership yields an alternate proof of
the result from [MinahanVolkovich2017] that $\SVparam{1}$—equivalently,
$\RFEparam{0}{1}$—hits every nonzero read-once formula $F$. Whereas the
original proof hinges on a clever ad-hoc argument, our proof (described in
s.test) is entirely systematic and amounts to a couple straightforward
observations in order to apply thm.ideal.membership.

As a proof of concept of the additional power of our characterization for
derandomization, we make progress in the model of read-once oblivious
algebraic branching programs (ROABPs).

######  Theorem 1.5 (ROABP hitting property).

For integer $l\geq 1$, $\SVl$ hits the class of polynomials computed by read-
once oblivious algebraic branching programs of width less than $(l/3)+1$ that
contain a monomial of degree at most $l+1$.

To the best of our knowledge, thm.roabp is incomparable to the known results
for ROABPs [RazShpilka2005, JQS2009, JQS2010, ForbesShpilka2013, FSS2014,
AGKS2014, AFSSV2015, GKST2015, GKS2016, GuoGurjar2020, SahaThankey2021,
BhargavaGhosh2021]. Without the restriction that the polynomial has a monomial
of degree at most $l+1$, thm.roabp would imply a fully blackbox polynomial-
time identity test for the class of constant-width ROABPs. No such test has
been proven to exist at this time; prior work requires either quasipolynomial
time or requires opening the blackbox, such as by knowing the order in which
the variables are read.

With the restriction, hitting the class in thm.roabp with a generator $G$
amounts to proving a lower bound on the width of ROABPs that compute nonzero
degree-$(l+1)$ polynomials in the vanishing ideal of $G$. For $G=\SVl$, such a
lower bound is interesting because there are nonzero degree-$(l+1)$
polynomials in the vanishing ideal, and, moreover, the vanishing ideal is
generated by such polynomials (thm.ideal.generators). That thm.roabp holds
suggests that the lower bound could hold for polynomials in $\VanIdeal{\SVl}$
of all degrees, in which case one would fully derandomize PIT for constant-
width ROABPs.

This stands in contrast to the case where $G=\SVparam{l+1}$. It is well-known
that $\SVparam{l+1}$ hits every polynomial containing a monomial of support
$l+1$ or less, and thus it hits the class in thm.roabp, irrespective of the
restriction on ROABP width. Indeed, there are _no_ nonzero polynomials of
degree $l+1$ or less in the vanishing ideal of $\SVparam{l+1}$. While the
lower bound necessary for $\SVparam{l+1}$ to hit the class in thm.roabp holds,
it does so only vacuously—a phenomenon that has no hope of extending to higher
degrees.

The method of proof of thm.roabp diverges significantly from prior uses of the
SV generator and therefore may be of independent interest. We elaborate on the
method more when we discuss the techniques of this paper, but for now, we
point out that most prior uses of the SV generator rely on combinatorial
arguments, i.e., arguments that depend only on which monomials are present in
the polynomials to hit. thm.roabp necessarily goes beyond this because there
is a polynomial in $\VanIdeal{\SVl}$ of degree $l+1$ that has the same
monomials as a polynomial computed by an ROABP of width 2, which by thm.roabp
is not in $\VanIdeal{\SVl}$ for $l\geq 4$. Namely, any instance of
$\SMVanOneparam{l-1}{l}$ contains exactly all the monomials of the form
$x_{i_{1}}\cdot x_{i_{2}}\cdot\dots\cdot x_{i_{l+1}}$ with
$(i_{1},\ldots,i_{l+1})\in X_{1}\times\cdots\times X_{l+1}$ for some disjoint
sets $X_{j}$; the same goes for $\prod_{j}\sum_{i_{j}\in X_{j}}x_{i_{j}}$,
which is computed by an ROABP of width 2.

#### Lower bounds.

Our result for ROABPs also illustrates this direction. Our derandomization
result for the class in thm.roabp is equivalent to a lower bound of at least
$(l/3)+1$ on the width of any ROABP computing a nonzero degree-$(l+1)$
polynomial in the vanishing ideal of $\SVl$, and in particular implies the
lower bound for $\VanGenparam{l-1}{l}$ and $\SMVanOneparam{l-1}{l}$. Other
hardness of representation results follow in a similar manner from prior
hitting properties of $\SV$ in the literature. The following lower bounds
apply to computing both $\VanGenparam{l-1}{l}$ and $\SMVanOneparam{l-1}{l}$:

  * •

Any syntactically multilinear formula must have at least
$\Omega(\log(l)/\log\log(l))$ reads of some variable [AvMV2011, Theorem 6.3].

  * •

Any sum of read-once formulas must have at least $\Omega(l)$ terms
[MinahanVolkovich2017, Corollary 5.2].

  * •

There exists an order of the variables such that any ROABP with that order
must have width at least $2^{\Omega(l)}$ [FSS2014, Corollary 4.3].

  * •

Any $\Sigma\mathrm{m}\bigwedge\Sigma\Pi^{O(1)}$ formula must have top fan-in
at least $2^{\Omega(l)}$ [Forbes2015]. See also [FSV2017arxiv, Lemma 5.12].

  * •

Lower bounds over characteristic zero for circuits with locally-low algebraic
rank [KumarSaraf2016, Lemma 5.2].

### Techniques.

Many of our results ultimately require showing that, under suitable
conditions, $\RFE$ hits a polynomial $p$. A recurring analysis fulfills this
role in the proofs of thm.ideal.generators,thm.ideal.membership,thm.roabp, as
well as several of the other results. We take intuition from the analytic
setting (e.g., $\FF=\RR$) and study the behavior of $p(\RFE)$ as a function of
the seed when the seed’s zeroes and poles are near the abscissas of chosen
variables of $p$. The behavior is dominated by the contributions of the
monomials of $p$ for which the variables with abscissas near zeros have
minimal degree and the variables with abscissas near poles have maximal
degree. Thus we may analyze a first approximation to $p(\RFE)$ by “zooming in”
on the contributions of the monomials in which the chosen variables have
extremal degrees. If the first approximation is nonzero, then we can conclude
that $\RFE$ hits $p$. We capture the technique in our Zoom Lemma
(Lemmalem.multi-zsub-deriv). Formal Laurent series can express the analytic
intuition purely algebraically. We provide a proof from first principles that
does not require any background in Laurent series and works over all fields.

thm.ideal.generators states the equality of two ideals,
$\langle\VanGenkl\rangle=\VanIdeal{\RFEkl}$, where $\langle\VanGenkl\rangle$
denotes the ideal generated by all instantiations of $\VanGenkl$, and
$\VanIdeal{\RFEkl}$ the vanishing ideal of $\RFEkl$.

  * •

The inclusion $\subseteq$ follows from linearizing the defining equations of
$\RFEkl$. The technique mirrors the use of resultants to compute implicit
equations for rational plane curves. This is where the univariate dependence
on the abscissas comes into play.

  * •

To establish the inclusion $\supseteq$ we first show that every equivalence
class of polynomials modulo $\langle\VanGenkl\rangle$ contains a
representative $p$ whose monomials exhibit the combinatorial structure of a
core. The structure allows us to apply lem.multi-zsub-deriv such that the
zoomed-in contributions of $p$ have only a single monomial. A single monomial
is nonzero at the evaluation in lem.multi-zsub-deriv, and we conclude that
$\RFEkl$ hits $p$.

The proof of thm.ideal.membership also relies on lem.multi-zsub-deriv.
Membership to the ideal is equivalent to the vanishing of all coefficients of
the expansion. The proof can be viewed as determining a small number of
coefficients sufficient to guarantee that their vanishing implies all
coefficients vanish. The restriction to multilinear polynomials $p$ allows us
to express the zoomed-in contributions of $p$ as the result of applying
partial derivatives and zero-substitutions.

thm.roabp makes use of the characterization of the minimum width of a read-
once oblivious algebraic branching program computing a polynomial $p$ as the
maximum rank of the monomial coefficient matrices of $p$ for various variable
partitions [Nisan1991]. We reduce to the case where $p$ is homogeneous of
degree $l+1$, whence the monomial coefficient matrices have a block-diagonal
structure. An application of lem.multi-zsub-deriv in the contrapositive yields
linear equations between elements of consecutive blocks under the assumption
that $\SVl$ fails to hit $p$. When some block is zero, the equations yield a
Cauchy system of equations on the rows or columns of its neighboring blocks;
since Cauchy systems have full rank, we deduce severe constraints on the row-
space/column-space of the neighboring blocks. A careful analysis turns this
observation into a rank lower bound of at least $(l/3)+1$ for a well-chosen
partition of the variables.

We point out that, in the preceding application, lem.multi-zsub-deriv is
instantiated several times in parallel to form a large system of equations on
the coefficients of $p$, and the whole system is necessary for the analysis.
This stands in contrast to most prior work using $\SV$, which uses knowledge
of how $p$ is computed to guide a search for a _single_ fruitful instantiation
of lem.multi-zsub-deriv.

### Alternating algebra representation.

The inspiration for several of our results stems from expressing the
polynomials $\VanGenkl$ using concepts from alternating algebra (also known as
exterior algebra or Grassmann algebra). In fact, the relationship between
thm.ideal.generators,thm.ideal.membership is based on the relationship
$\bdry^{2}=0$ from alternating algebra. Our original statement and proof of
the theorem made use of that framework, but we managed to eliminate the
alternating algebra afterwards. Still, as we find the perspective insightful
and potentially helpful for future developments, we describe the connection
briefly here and in more detail in s.simplicial-repr. We explain the intuition
behind thm.ideal.membership for the simple case where the degree of the
polynomial $p$ equals $l+1$. In that setting, belonging to the ideal generated
by the polynomials $\VanGenkl$ is equivalent to being in their linear span.

The alternating algebra $A$ of a vector space $V$ over a field $\FF$ consists
of the closure of $V$ under an additional binary operation, referred to as
“wedge” and denoted $\wedge$, which is bilinear, associative, and satisfies

$v\wedge v=0$ (7)

for every $v\in V$. This determines a well-defined algebra. When the
characteristic of $\FF$ is not 2, this can equivalently be understood as

$v_{1}\wedge v_{2}=-(v_{2}\wedge v_{1})$ (8)

for every $v_{1},v_{2}\in V$. In any case, for any $v_{1},v_{2},\dots,v_{k}\in
V$,

$v_{1}\wedge v_{2}\wedge\dots\wedge v_{k}$ (9)

is nonzero iff the $v_{i}$’s are linearly independent, and any permutation of
the order of the vectors in eq.wedge yields the same element of $A$ up to a
sign. The sign equals the sign of the permutation, whence the name
“alternating algebra.” If $V$ has a basis $X$ of size $n$, then a basis for
$A$ can be formed by all $2^{n}$ expressions of the form eq.wedge where the
$v_{i}$’s range over all subsets of $X$ and are taken in some fixed order.
Considering the elements of $X$ as vertices, the basis elements of $A$ can be
thought of as the oriented simplices of all dimensions that can be built from
$X$.

Anti-commutativity, the relation eq.anti-commute, arises naturally in the
context of network flow, where $X$ denotes the vertices of the underlying
graph, and a wedge $v_{1}\wedge v_{2}$ of level $k=2$ represents one unit of
flow from $v_{1}$ to $v_{2}$. Equation eq.anti-commute reflects the fact that
one unit of flow from $v_{1}$ to $v_{2}$ cancels with one unit of flow from
$v_{2}$ to $v_{1}$. The adjacent levels $k=1$ and $k=3$ also have natural
interpretations in the flow setting: $v_{1}$ (the element of $A$ of the form
eq.wedge with $k=1$) represents one unit of surplus flow at $v_{1}$ (the
vertex of the graph), and $v_{1}\wedge v_{2}\wedge v_{3}$ abstracts an
elementary circulation of one unit along the directed cycle $v_{1}\to v_{2}\to
v_{3}\to v_{1}$.

The different levels are related by so-called boundary maps. Boundary maps are
linear transformations that map a simplex to a linear combination of its
subsimplices of one dimension less. The maps are parametrized by a weight
function $w:X\to\FF$, and defined by

$\bdry_{w}:v_{1}\wedge v_{2}\wedge\dots\wedge
v_{m}\mapsto\sum_{i=1}^{m}(-1)^{i+1}w(v_{i})\;v_{1}\wedge\dots\wedge
v_{i-1}\wedge v_{i+1}\wedge\dots\wedge v_{m},$ (10)

an expression resembling the Laplace expansion of a determinant along a column
$[w(v_{i})]_{i=1}^{m}$. In the flow setting, using $w\equiv 1$, applying
$\bdry_{1}$ to $v_{1}\wedge v_{2}$ yields $v_{2}-v_{1}$, the superposition of
demand at $v_{1}$ and surplus at $v_{2}$ corresponding to one unit of flow
from $v_{1}$ to $v_{2}$. Likewise, $\bdry_{1}$ sends the abstract elementary
cycle $v_{1}\wedge v_{2}\wedge v_{3}$ to the superposition of the three edge
flows that make up the cycle. A linear combination $p$ of terms eq.wedge with
$k=2$ represents a valid circulation iff it satisfies conservation of flow at
every vertex, which can be expressed as $\bdry_{1}(p)=0$, i.e., $p$ is in the
kernel of $\bdry_{1}$. An equivalent criterion is for $p$ to be the
superposition of elementary circulations, which can be expressed as $p$ being
in the image of $\bdry_{1}$. The relationship between the image and the kernel
of boundary maps holds in general: For any linearly independent
$w_{0},\ldots,w_{m}$, it holds that

$\im\left(\bdry_{w_{m}}\circ\bdry_{w_{m-1}}\circ\dots\circ\bdry_{w_{0}}\right)\;=\;\bigcap_{i=0}^{m}\ker\left(\bdry_{w_{i}}\right).$
(11)

(When $w_{0},\ldots,w_{m}$ are linearly dependent,
$\bdry_{w_{m}}\circ\dots\circ\bdry_{w_{0}}$ is the zero map.)

In the context of $\RFE$, the set $X$ creates a vertex for each variable, and
simplices correspond to multilinear monomials. The anti-commutativity of
$\wedge$ coincides with the fact that swapping two arguments to $\VanGenkl$
means swapping two rows in eq.def.gen, which changes the sign of the
determinant. Using the above boundary maps, the right-hand side of eq.def.gen
can be viewed as $\bdry_{\omega}(v_{i_{1}}\wedge v_{i_{2}}\wedge\dots\wedge
v_{i_{k+l+2}})$, where
$\bdry_{\omega}\doteq\bdry_{w_{k}}\circ\bdry_{w_{k-1}}\circ\dots\circ\bdry_{w_{0}}$
and $w_{d}(v_{i})\doteq(a_{i})^{d}$. By eq.image.kernel, this means that
$\VanGenkl$ is in the kernel of $\bdry_{w_{d}}$ for each
$d\in\\{0,1,\dots,k\\}$, or equivalently, in the kernel of
$\bdry_{\widetilde{w}}$ for each $\widetilde{w}:X\to\FF$ of the form
$\widetilde{w}(v_{i})=w(a_{i})$ where $w$ is a polynomial of degree at most
$k$. This is precisely the condition eq.generator.dependencies. In fact,
eq.image.kernel implies that the linear span of the generators $\VanGenkl$
consists exactly of the polynomials of degree $l+1$ in this kernel. The latter
condition is precisely what thm.ideal.membership.cond2 in thm.ideal.membership
expresses.

### Further research.

In this paper we propose to investigate the power of generators through
characterizations of their vanishing ideals, develop such characterizations
for the generators $\SV$ and $\RFE$, and initiate a systematic study of the
repercussions. Completing the study entails going over classes $\Class$ of
polynomials of interest, checking whether they intersect nontrivially with the
vanishing ideal, and, if so, checking whether the intersection can be hit by
other generators. One specific target is the elimination of the degree
restriction from thm.roabp.

In addition to Shpilka–Volkovich, other specific generators for which we
suggest to characterize their vanishing ideals include Klivans–Spielman and
Gabizon–Raz, as well as generators based on isolating weight assignments find
recurring use in the PIT literature. Other possibilities may include the
$\RFE$ generator with pseudorandom abscissas, or higher-level work that
relates the vanishing ideal of a combination of generators to the vanishing
ideals of the constituent generators.

Along a related line of thought, the generator $\SVl$ is the canonical example
of an $l$-wise independent generator in the algebraic setting. Understanding
the power such generators more broadly should lead to useful insights for
derandomizing PIT. This work demonstrates explicit polynomials like
$\VanGenparam{l-1}{l}$ and $\SMVanOneparam{l-1}{l}$ that are not automatically
hit by $l$-wise independence; indeed they are not hit by $\SVl$. Is there a
deeper underlying reason related to $l$-wise independence?

Lastly, we remark that our derivation of $\VanGen$ from $\RFE$ in s.gens can
be abstracted to construct a polynomial $p$ in $\VanIdeal{G}$ for any
generator $G$, provided that $G$ substitutes rational functions for which the
numerators and denominators belong to linear spaces whose dimensions sum to
less than $n$. As the resulting polynomial $p$ is the determinant of an
$n\times n$ or smaller matrix, such a generator $G$ cannot derandomize PIT for
models that compute such determinants. One may view this phenomenon as a
possible barrier for derandomization. Our derivation works for the $\RFE$
generator, but does not directly apply to the SV generator, even though the
two are equivalent in derandomization power. It would be useful to have
sufficient conditions under which a generator and _all_ its reparametrizations
avoid the barrier.

### Organization.

We construct the generating set for the vanishing ideal (thm.ideal.generators)
in s.gens, followed by lem.multi-zsub-deriv in s.zoom. The ideal membership
test (thm.ideal.membership) is developed in s.test. We present the results on
sparseness in s.sparse, and the ones on set-multilinearity in s.set-
multiaffine. Background on ROABPs and our result on derandomizing PIT for
ROABPs (thm.roabp) are covered in s.constant-width-roabp. We end our paper in
s.simplicial-repr with a further discussion of the alternating algebra
representation. The appendix contains some details and a formal treatment of
the relationship between $\RFE$ and $\SV$.


# From integrability to chaos in quantum Liouvillians

Álvaro Rubio-García<EMAIL_ADDRESS>Instituto de Estructura de la
Materia, IEM-CSIC, Serrano 123, Madrid, E-28006, Spain  Rafael A. Molina
<EMAIL_ADDRESS>Instituto de Estructura de la Materia, IEM-CSIC, Serrano
123, Madrid, E-28006, Spain  Jorge Dukelsky<EMAIL_ADDRESS>Instituto de
Estructura de la Materia, IEM-CSIC, Serrano 123, Madrid, E-28006, Spain

###### Abstract

The dynamics of open quantum systems can be described by a Liouvillian, which
in the Markovian approximation fulfills the Lindblad master equation. We
present a family of integrable many-body Liouvillians based on Richardson-
Gaudin models with a complex structure of the jump operators. Making use of
this new region of integrability, we study the transition to chaos in terms of
a two-parameter Liouvillian. The transition is characterized by the spectral
statistics of the complex eigenvalues of the Liouvillian operators using the
nearest neighbor spacing distribution and by the ratios between eigenvalue
distances.

##  I Introduction

Quantum chaos has gained a renewed interest in the last few years due to its
relevance to important current topics of research like quantum thermalization,
non-equilibrium dynamics of quantum many-body systems and quantum information
[1, 2, 3, 4, 5]. In the standard Hermitian quantum mechanics valid for closed
systems, the dynamical regime of the system can be characterized through
spectral statistics [6]. Integrable systems present level clustering and the
nearest-neighbor spacing distribution follows the one-dimensional Poisson
distribution $P(s)=e^{-s}$ [7], while chaotic systems present level repulsion
with the $P(s)$ close to the Wigner surmise of Random Matrix Theory (RMT)
depending on their symmetry class, $P(s)\propto s^{\beta}$ for small $s$, with
$\beta=1,2,4$ for orthogonal, unitary, and symplectic symmetries respectively,
which is the content of the famous Bohigas-Giannoni-Schmit (BGS) conjecture
[8]. The BGS conjecture is well founded now in the semiclassical theory, valid
for systems with a proper classical limit [9, 10, 11] and supported by
overwhelming numerical and experimental evidence in many different quantum
systems [12, 13, 14]. The situation in many-body quantum systems is much less
clear, although some theoretical advances have been recently made [15, 16,
17]. Due to the symmetry under exchange of fermionic or bosonic particles, the
classical limit cannot be properly defined. Usually the BGS conjecture is
assumed to hold also for many-body quantum systems, mainly based on numerical
results, but a rigorous derivation is still lacking. The transition between
the integrable and chaotic universal limits is non-universal, depending on the
specifics of the particular system under study, although it has been explored
in much detail for different systems [18, 19]. In the transition between the
integrable and chaotic orthogonal cases, for example, some systems present
fractional level repulsion with $P(s)\propto s^{\beta}$ with the value of
$\beta$ varying continuously between the integrable case $\beta=0$ and the
corresponding RMT ensemble value $\beta=1$, while others present full level
repulsion but only for a fraction of levels [20]. Many systems, specially in
the many-body case, show the former behavior. However, the semiclassical
theory of Berry and Robnik for the transition predicts the latter one [19]. In
this case $P(0)=F$, with $F$ given by the fraction of regular orbits in phase
space of the classical limit for the model under consideration.

In open quantum systems, the theory is much less developed, even if the first
results came shortly after the proposal of the BGS conjecture [21]. Open
quantum systems, can be described by the Liouville equation, which
characterizes the time evolution of the density matrix operator. In the
Markovian approximation, the Liouvillian is a linear non-Hermitian operator,
and the Liouville equation can be written as a Lindblad master equation [22].
The Liouvillian, then, has complex eigenvalues instead of the real energies of
standard Hermitian quantum mechanics. The initial approach to the problem was
to study integrable or chaotic Hamiltonians with a weak coupling to the
environment. When the Hamiltonian was integrable, Grobe _et al._ studied the
spectral statistics in the complex plane and found good agreement with the
two-dimensional Poisson distribution [21]. In the chaotic limit there appears
universal cubic repulsion $P(s)\propto s^{3}$ for small values of $s$ as in
the Ginibre ensemble of non-Hermitian random matrices [23], although the
details of the full $P(s)$ distribution depend on the symmetries of the non-
Hermitian matrix [24, 25]. For an open quantum spin chain, the level spacing
distribution in the transition from integrability to chaos has been fitted by
a static two-dimensional Coulomb gas with harmonic confinement where level
repulsion is given by the inverse temperature showing fractional level
repulsion in the transition [26].

Very recently, the need for a different approach in the study of the
integrable and chaotic properties of open quantum systems has been triggered
as a consequence of the discovery of new families of integrable many-body
Liouvillians [27, 28, 29]. Extending the class of exactly solvable and quantum
integrable Liouvillians is an important step towards improving our
understanding of open quantum many-body systems. The statistical properties of
the complex spectra of random chaotic Liouvillians have been studied in a few
recent works [30, 31]. However, the transition between the exactly solvable
integrable limit and the chaotic limit in physical many-body Liouvillians
remains mostly unexplored.

In this letter we will extend the model of Ref. [28] based on the SU(2) spin 1
Richardson model to a line of integrability within the rational Richardson-
Gaudin (RG) class of integrable models. This new family of integrable
Liouvillians has a rich and complex structure of the jump operators and allows
for a simple parametrization along the line of integrability. We then define
in terms of a single parameter a Liouvillian that interpolates between
integrability and a fully chaotic limit. With these model Liouvillians we
characterize these transitions through the spectral statistics of the complex
eigenvalues of the Liouvillian operators and by the average properties of the
complex spacing ratios.

##  II Integrable Richardson-Gaudin Liouvillians

We consider a system described by a Hamiltonian $H$ weakly coupled to a
Markovian environment. The evolution of its density matrix can be expressed in
terms of the Lindblad master equation [22]

$\begin{split}\partial_{t}\rho=\mathcal{L}(\rho)=&-i\left[H,\rho\right]\\\
&+\sum_{j}\left(L_{j}\rho L_{j}^{\dagger}-\frac{1}{2}\rho
L_{j}^{\dagger}L_{j}-\frac{1}{2}L_{j}^{\dagger}L_{j}\rho\right),\end{split}$
(1)

where $\mathcal{L}$ is the Liouvillian superoperator acting on the space of
density matrices. The first term on the right hand side is the Hamiltonian
commutator describing the unitary evolution and the second term takes into
account the interaction with the environment through the Lindblad jump
operators $L_{j}$. The time evolution of a Liouvillian eigenstate
$\mathcal{L}(\rho_{i})=\lambda_{i}\rho_{i}$ is given by
$e^{\mathcal{L}t}\rho_{i}=e^{\lambda_{i}t}\rho_{i}$ and it can be shown that
the real part of the eigenvalues is always non positive,
$\textrm{Re}\left\\{\lambda_{i}\right\\}\leq 0$ [22, 32, 33]. Furthermore, the
Liouvillian spectrum consists of complex conjugate pairs
$(\rho_{i},\lambda_{i}),\,(\rho_{i}^{\dagger},\,\lambda_{i}^{*})$, implying
that a real eigenvalue corresponds to an hermitian density matrix. There is at
least one state with zero Liouvillian eigenvalue, the  _steady state_
$\rho_{\textrm{SS}}$, to which any initial state decays in the long time limit
$\lim_{t\rightarrow\infty}e^{\mathcal{L}t}\rho=\rho_{\textrm{SS}}$. As a
consequence, the trace of any eigenstate with negative real part must be zero.
The steady state can also be degenerate, in which case the system will decay
to a linear combination of states within the degenerate manifold depending on
the initial state.

Similar to Ref. [28], we study a chain of $L$ spins subject to a non uniform
magnetic field $h_{i}$

$H=\sum_{i=1}^{L}h_{i}S_{i}^{z},$ (2)

whose interaction with the environment is mediated by quantum jumps of the
form

$L_{z}^{a}=\sqrt{2\Gamma_{0}}\sum_{i=1}^{L}z_{i}^{a}S_{i}^{z},\quad
L_{\pm}=\sqrt{\Gamma_{\pm}}\sum_{i=1}^{L}x_{i}S_{i}^{\pm},$ (3)

where $a$ labels a set of different jumps and $z_{i}^{a}$ and $x_{i}$ are
arbitrary jump amplitudes. As we will see, the equivalence of the gain $L_{+}$
and loss $L_{-}$ operator strengths, $\Gamma=\Gamma_{+}=\Gamma_{-}$, is a
necessary condition for quantum integrability. In the particular case of a
single collective spin Hamiltonian, this condition could be relaxed [29]. Note
that for this choice the jumps become hermitian,
$\left(L_{z}^{a}\right)^{\dagger}=L_{z}^{a},\,\left(L_{+}\right)^{\dagger}=L_{-}$,
which immediately leads to a trivial steady state
$\rho_{\text{SS}}=\mathbb{I}$. The noisy spin model of Ref. [28] is a
particular case of the general RG integrable models corresponding to the limit
$x_{i}=1,\,z_{i}^{a}=\delta_{a,1}$ (see Eq. (12) below).

The vector representation of the Lindblad master equation [33] doubles the
Hilbert space $\mathcal{H}$ of dimension $\mathcal{N}$ by mapping the density
matrix $\rho$ into a vector $|\rho\rangle$ in the space $\mathcal{H\otimes H}$
of dimension $\mathcal{N}^{2}$. The mapping transforms every spin operator $S$
acting to the right or to the left of the density matrix into a superoperator
of the form

$\begin{split}S\rho\rightarrow&\ S\otimes\mathbb{I}\,|\rho\rangle\ \hat{=}\
S\,|\rho\rangle\\\ \rho S\rightarrow&\ \mathbb{I}\otimes S^{T}\,|\rho\rangle\
\hat{=}\ J\,|\rho\rangle\end{split}$ (4)

Thus, the operators acting to the right of the density matrix are mapped to a
new set of spin operators $J_{i}$ acting on a dual space with the same
dimension $\mathcal{N}$ of the original space. This mapping allows the
Lindblad master equation to be written in a matrix-vector form
$\mathcal{L}|\rho\rangle$ with a Liouvillian matrix $\mathcal{L}$ of dimension
$\mathcal{N}^{2}$ whose eigenvalues can be computed by an exact
diagonalization for small size systems or by an appropriate many-body
approximation [34, 35, 36, 37].

The complete Liouvillian in the Lindblad form is written, after a $\pi$
rotation around the $y$ axis in the dual space, as

$\begin{split}\mathcal{L}=&-i\sum_{i}h_{i}\left(S_{i}^{z}+J_{i}^{z}\right)\\\
&-\Gamma_{0}\sum_{a=1}^{n_{j}}\sum_{ij}z_{i}^{a}z_{i}^{a}\left(S_{i}^{z}+J_{i}^{z}\right)\left(S_{j}^{z}+J_{j}^{z}\right)\\\
&-\Gamma\sum_{ij}x_{i}x_{j}\frac{1}{2}\left[\left(S_{i}^{+}+J_{i}^{+}\right)\left(S_{j}^{-}+J_{j}^{-}\right)+h.c.\right]\end{split}$
(5)

The fact that the spin $S_{i}$ and the dual spin $J_{i}$ appear always as a
sum is a consequence of the equality of the gain and loss
$\Gamma_{+}\equiv\Gamma_{-}$ jumps operators. We can, therefore, define the
total spin operators $K_{i}=S_{i}+J_{i}$ for each lattice site $i$. In terms
of these new total spin operators the Liouvillian acquires the simpler
expression

$\begin{split}\mathcal{L}=&-i\sum_{i}h_{i}K_{i}^{z}-\Gamma\sum_{i}x_{i}^{2}\left(K_{i}\right)^{2}\\\
&+\sum_{i}\left[\Gamma
x_{i}^{2}-\Gamma_{0}\sum_{a}\left(z_{i}^{a}\right)^{2}\right]\left(K_{i}^{z}\right)^{2}\\\
&-\Gamma_{0}\sum_{a,i\neq j}z_{i}^{a}z_{i}^{a}K_{i}^{z}K_{j}^{z}\\\
&-\Gamma\sum_{i\neq
j}x_{i}x_{j}\frac{1}{2}\left(K_{i}^{+}K_{j}^{-}+K_{i}^{-}K_{j}^{+}\right).\end{split}$
(6)

We stress that the magnitude of the spins at each site of the chain $s_{i}$ is
completely arbitrary, with the only constraint imposed by the mapping that it
equals the magnitude of the dual spin $s_{i}=j_{i}$. Hence, the total spin at
site $i$ can take the following values $k_{i}=0,1,\cdots,2s_{i}$.

The Liouvillian (6) presents a set of $L+1$ weak symmetries [32], given by
operators that commute with the total Liouvillian
$\left[\mathcal{L},O\right]=0$. These weak symmetries are the total
magnetization operator $K^{z}=\sum_{i}K_{i}^{z}$ and the $L$ SU(2) Casimir
operators at each site

$K_{i}^{2}=\left(K_{i}^{z}\right)^{2}+\frac{1}{2}\left(K_{i}^{+}K_{i}^{-}+K_{i}^{-}K_{i}^{+}\right)=k_{i}\left(k_{i}+1\right)$
(7)

Since all symmetry operators commute $[K^{z},K^{2}_{i}]=0$, the Liouvillian
matrix is separated into blocks labelled by the set of $L$ quantum numbers
$k_{i}$ and each one, in turn, divided into sub-blocks labelled by the total
$k_{z}$. Due to the symmetry of the Liouvillian, the occurrence of a singlet
total spin $k_{i}=0$ at site $i$ effectively removes this site from the chain,
since the action of any SU(2) generator on this state annihilates it.
Therefore, the state with all singlet spins is the steady state of the
Liouvillian with 0 eigenvalue.

For the particular case of our Liouvillian (6) with an inhomogeneous magnetic
field along the $z$ axis, the eigenstates at opposite symmetry sectors $k_{z}$
and $-k_{z}$ are related by an anti-unitary operator $F=\mathcal{C}e^{i\pi
K^{x}}$, where $\mathcal{C}$ is the complex conjugation operator and $e^{i\pi
K^{x}}$ performs a $\pi$ rotation around the $x$ axis reversing the spins
along the $y$ and $z$ directions. The operator $F$ commutes with the
Liouvillian and with the $L$ Casimir operators $K^{2}_{i}$ but it does not
commute with the $K^{z}$ symmetry. In spite of the fact that $F$ cannot be
added to the chain of weak symmetries, it is useful to understand the
structure of the spectrum in each $k_{z}$ block. For every Liouvillian
eigenstate $|\rho_{i}\rangle$ (in vector form) with complex eigenvalue
$\lambda_{i}$, the eigenstate $|F\rho_{i}\rangle$ is also a Liouvillian
eigenstate with conjugate eigenvalue

$\mathcal{L}|F\rho_{i}\rangle=\mathcal{L}F|\rho_{i}\rangle=F\mathcal{L}|\rho_{i}\rangle=F\lambda_{i}|\rho_{i}\rangle=\lambda_{i}^{*}|F\rho_{i}\rangle.$
(8)

Moreover, the $|F\rho_{i}\rangle$ state has opposite total angular momentum
$K^{z}$ to $|\rho_{i}\rangle$

$\langle
F\rho_{i}|K^{z}|F\rho_{i}\rangle=\langle\rho_{i}|F^{\dagger}K^{z}F|\rho_{i}\rangle=\langle\rho_{i}|\left(-K^{z}\right)|\rho_{i}\rangle.$
(9)

Thus, each Liouvillian block with momentum $k_{z}$ contains the complex
conjugate eigenvalues of the corresponding block with momentum $-k_{z}$.
Blocks with $k_{z}=0$ are exceptional, since they contain all complex
conjugate pairs. This could be explained by the fact that in the subspaces of
$k_{z}=0$ the operators $F$ and $K^{z}$ do commute.

The Liouvillian (6) can be derived from the integrals of motion (IOM) of the
Richardson-Gaudin models [38, 39]

$R_{i}=K_{i}^{z}-G\sum_{j\left(\neq
i\right)=1}^{L}\frac{X_{ij}}{2}\left(K_{i}^{+}K_{j}^{-}+K_{j}^{+}K_{i}^{-}\right)+Z_{ij}K_{i}^{z}K_{j}^{z}.$
(10)

The local spins $k_{i}$ are in principle arbitrary. A possible choice for the
matrices $X$ and $Z$ such that the $L$ IOM commute among themselves,
$\left[R_{i},\,R_{j}\right]=0$, is [39]

$\begin{split}X_{ij}=&\
\frac{\sqrt{(1-\alpha)+\alpha\,\eta_{i}^{2}}\sqrt{(1-\alpha)+\alpha\,\eta_{j}^{2}}}{\eta_{i}-\eta_{j}}\\\
Z_{ij}=&\
\frac{(1-\alpha)+\alpha\,\eta_{i}\eta_{j}}{\eta_{i}-\eta_{j}},\end{split}$
(11)

with $\alpha\in[0,1]$ an interpolation variable within the region of
integrability. Note that the limits $\alpha=0,1$ assure the equality of the
two matrices $X$ and $Z$ and pertain to the rational or XXX RG family [40].
The IOM with $\alpha=0$ lead to the constant pairing Hamiltonian originally
solved by Richardson [41], while $\alpha=1$ leads to the recently studied
separable pairing Hamiltonian [42]. For intermediate values of $\alpha$ the
matrices $X$ and $Z$ are not equal, corresponding to the XXZ or hyperbolic RG
family. Moreover, the IOM have a set of $L+1$ parameters $(G,\eta_{i})$ that
are arbitrary complex numbers. In what follows, we consider the
$\eta^{\prime}s$ as real parameters and $G$ as pure imaginary, $G=ig$. While
this choice makes the IOM no longer Hermitian, the integrability properties of
these models are kept intact. In terms of these IOM we can obtain the set of
integrable Liouvillians (6) with the following linear combination

$\begin{split}\mathcal{L}_{\textrm{int}}=-i\sum_{i=1}^{L}\eta_{i}R_{i}=&-i\sum_{i}^{L}\eta_{i}K_{i}^{z}\\\
&-\frac{g}{2}\sum_{i\neq
j}^{L}\left[(1-\alpha)+\alpha\,\eta_{i}\eta_{j}\right]K_{i}^{z}K_{j}^{z}\\\
&-\frac{g}{2}\sum_{i\neq
j}^{L}x_{i}x_{j}\frac{1}{2}\left(K_{i}^{+}K_{j}^{-}+K_{i}^{-}K_{j}^{+}\right)\end{split}$
(12)

where $x_{i}=\sqrt{(1-\alpha)+\alpha\,\eta_{i}^{2}}$.

This family of integrable Liouvillians can be derived from am open quantum
system with Hamiltonian $H=\sum_{i}\eta_{i}S_{i}^{z}$, gain and loss jumps (3)

$L_{\pm}=\sqrt{\frac{g}{2}}\sum_{i=1}^{L}\sqrt{(1-\alpha)+\alpha\,\eta_{i}^{2}}\,S_{i}^{\pm}$
(13)

and two dephasing jumps

$L_{z}^{1}=\sqrt{(1-\alpha)\,\frac{g}{2}}\sum_{i=1}^{L}\,S_{i}^{z},\quad
L_{z}^{2}=\sqrt{\alpha\,\frac{g}{2}}\sum_{i=1}^{L}\,\eta_{i}S_{i}^{z}$ (14)

Figure 1: (a-c) Integrable and chaotic Liouvillian spectra for
$L=6,\,K^{z}=1,\,k_{i}=1$ in the rational constant (a), rational separable (b)
and chaotic (c) Liouvillian limits. (d-f) Level spacing distribution $P(s)$ of
the integrable and chaotic Liouvillians for different sizes $L=7,9,11$ and
$k_{i}=1,k_{z}=1$ with the same ordering as in (a-c). The dashed and dash-
dotted lines show the $\gamma=0$ and $\gamma=1$ limits of the interpolation
distribution $P(s;\gamma)$ (19), respectively, while the dotted line shows the
level statistics of the $AI^{\dagger}$ universality class.

In order to break integrability we will, in addition, consider a set of
$n_{j}$ random jumps that lead to a chaotic (non-integrable) Liouvillian.
These jumps are given by

$L_{\pm}^{a}=\sqrt{\Gamma}\sum_{i=1}^{L}w_{i}^{a}S_{i}^{\pm},\quad
a=1,\dots,n_{j},$ (15)

with $(w_{1}^{a},w_{2}^{a},\dots)$ a set of random orthonormal vectors that
are also orthonormal to the integrable gain and loss jump vector
$(x_{1},x_{2},\dots)$. Our results show that this choice of orthonormal
$w^{a}$ vectors displays a more definite chaotic behavior than non orthogonal
ones. The chaotic Liouvillian is then

$\begin{split}\mathcal{L}_{\textrm{chaotic}}=&-i\sum_{i}\eta_{i}K_{i}^{z}+\Gamma\sum_{a=1}^{n_{j}}\sum_{i}\left(w_{i}^{a}\right)^{2}\left[\left(K_{i}^{z}\right)^{2}-\left(K_{i}\right)^{2}\right]\\\
&-\Gamma\sum_{a=1}^{n_{j}}\sum_{i\neq
j}w_{i}^{a}w_{j}^{a}\frac{1}{2}\left(K_{i}^{+}K_{j}^{-}+K_{i}^{-}K_{j}^{+}\right),\end{split}$
(16)

When choosing the number of random orthogonal jumps one has to take into
account that when $n_{j}$ is equal to the total number of spins $L$, the full
Liouvillian becomes diagonal and, thus, integrable. So, the most chaotic
spectral statistics occurs for $n_{j}\approx L/2$ and this has been our choice
in Eq. 16.

The transition from integrability to chaos will be characterized by two
parameters $(\alpha,\beta)$, that can interpolate between the two rational
integrable liouvillians and the fully chaotic limit.

$\mathcal{L}(\alpha,\beta)=(1-\beta)\mathcal{L}_{\text{int}}(\alpha)+\beta\,\mathcal{L}_{\text{chaotic}},$
(17)

with $(\alpha,\beta)=(0,0)$ the rational constant model, $(1,0)$ the separable
model and $(\alpha,1)$ the chaotic limit.

##  III Liouvillians Spectral statistics in the transition from integrability
to chaos

In this work we study the spectrum of Liouvillians of spin 1/2 chains. Because
the $L+1$ weak symmetries block diagonalize the Liouvillian space, we will
focus on the $k_{i}=1,k_{z}=1$ subspaces, as this represents the block with
maximal dimension, adequate for the study of spectral statistics. While the
subspace $k_{z}=0$ is bigger, in principle, than $k_{z}=1$, its spectrum is
divided in half by the $F$ symmetry and results into two blocks of smaller
dimension. To preserve the scaling of the Liouvillian spectral characteristics
with the system size we fix $\Gamma=\frac{1}{L}$ in the subsequent
computations.

Figs. 1(a,b) show some typical examples of the spectrum of the rational
constant and separable models, respectively, for Liouvillians with $L=6$,
while Fig. 1(c) shows the spectrum of the chaotic Liouvillian limit with the
same parameters and $n_{j}=L/2=3$ random jumps. The parameters $\eta_{i}$ are
chosen as $\eta_{i}=1/2+i/L+\xi_{i}$, where the values of $\xi_{i}$ are random
variables coming from an uniform distribution $\xi_{i}\in[-1/L,1/L]$. This is
a convenient choice for random values of $\eta_{i}$ that approximately keeps
the relative strength of the Hamiltonian and the jumps constant as we increase
size. For the chaotic jumps, the values $\omega_{i}^{a}$ are drawn from a
uniform distribution $[-1/2,1/2]$ and then orthonormalized using a Grand-
Schmidt procedure. The ordered patterns of the two integrable spectra, as we
will see below, have random local distribution of the eigenvalues, while the
chaotic spectrum presents level repulsion.

Figure 2: Characterizations of the integrable to chaotic transition using the
$\gamma$ 19 parameter (a) and the $\langle r\rangle$ (b) and
$-\langle\cos\theta\rangle$ (c) averages of the complex spacing ratios 20.
Each point represents a sample of 20 Liouvillians of matrix dimension $8350$.

To characterize the integrable to chaos Liouvillian transition we measure the
nearest neighbor level statistics. Before that, an unfolding procedure is done
in order to rescale the local level density to unity, which results in
adimensional quantities to be able to compare different Liouvillian systems to
the universal results of Random Matrix Theory [6]. The unfolding procedure can
be tricky and if it is not performed carefully may give rise to misleading
results [43]. We use a local unfolding procedure where the level distances are
rescaled by the local level density [21]. This local density is computed by
calculating the area of a circle containing the $n$ nearest neighbors of each
level. The nearest-neighbor spacings are then

$S_{i}=\sqrt{\frac{n}{\pi d_{i,n}^{2}}}d_{i,1},$ (18)

where $d_{i,n}$ is the distance in the complex plane between level $i$ and its
$n$-th neighbor. Then, the value of $S_{i}$ is rescaled by its mean to obtain
the final unfolded spacings $s_{i}$, $s_{i}=S_{i}/\left<S_{i}\right>$. The
unfolded spacings, by construction, have $\left<s\right>=1$ as required for a
proper comparison. The distribution of $s_{i}$, $P(s)$, is then contrasted
with the theoretical results as explained in the paragraph below. There is a
certain arbitrariness in choosing the number $n$, but it has to be large
enough to not depend on small scale fluctuations and small enough to capture
the local statistics. We have found $n=20$ to be a good compromise between
these two limits, the results being quite stable and robust in a wide range
around this value.

For integrable Liouvillians, the level statistics follow a Poisson
distribution in the plane
$P(s)=\frac{\pi}{2}s\exp\left(-\frac{\pi}{4}s^{2}\right)$, while the level
statistics of chaotic Liouvillians belong to the $AI^{\dagger}$ universality
class of non-Hermitian symmetric matrices [24, 25]. To observe the transition
from quantum integrability to chaos we fit the level spacing distribution
$P(s)$ to a distribution governed by a parameter $\gamma$

$P(s;\gamma)=A(\gamma)s^{2\gamma+1}e^{-B(\gamma)s^{2}},\quad\gamma\in\left[0,1\right]$
(19)

with $A(\gamma),B(\gamma)$ such that
$\int_{0}^{\infty}P(s)ds=\int_{0}^{\infty}sP(s)ds=1$. Similar to the Brody
distribution for the transition between integrability and chaos in standard
Hermitian Hamiltonians [18], the limits of this distribution go from the
integrable limit $\gamma=0$, which corresponds to the Poisson distribution in
the plane, to the chaotic limit $\gamma=1$, corresponding to the Wigner
surmise of the level spacing distribution of the Ginibre ensemble [23]. Figs.
1(d-f) shows the limit $\gamma=0$ (dashed line) and $\gamma=1$ (dash-dotted).
While the distribution in the limit $\gamma=1$ does not exactly correspond to
the level statistics of the $AI^{\dagger}$ universality class (dotted line in
Fig. 1), it shows the same level repulsion $P(s)\sim s^{3}$ and both
distributions are sufficiently similar so that $\gamma=1$ is a good
characterization of the chaotic limit.

The unfolded level statistics of the two rational and chaotic limits are shown
in Figs. 1(d-f) for Liouvillians with sizes $L=7,9,11$ and Hilbert space
dimension $357,\,2907$ and $24068$, respectively. For each of these sizes $L$
we compute a sample of several Liouvillian spectra so that there are in total
more than 700.000 eigenvalues to draw statistics from. The two integrable
limits show levels statistics that are very close to the 2d Poisson
distribution in the plane, indicating a random distribution of the local
eigenvalues, while the chaotic limit shows similar level statistics to the
$AI^{\dagger}$ universality class. For each Liouvillian limit there is also a
convergence to the respective 2d Poisson and $AI^{\dagger}$ distributions from
smaller to larger $L$’s.

Fig. 2(a) shows the transition from integrable to chaotic level statistics
characterized by the interpolation parameter $\gamma$, which we fit to the
level spacings using a Maximum Likelihood Estimation. Each value of the
parameters $(\alpha,\beta)$ represents the statistics of a sample of 20
Liouvillians of size $L=10$, with a Hilbert space dimension of $8350$
eigenvalues, and $n_{j}=5$ random jumps. For values of $\beta=0$ the
transition parameter $\gamma$ is very close to 0, showing integrable
statistics for both rational integrable Liouvillians, and for $\beta>0$ there
is a smooth transition in terms of $\gamma$ from integrable to chaotic
statistics up to the limit $\beta=1$, with the integrability breaking at
relatively low values of $\beta$. This transition is shown in more detail in
Fig. 3 for different system sizes $L$ and for $\alpha=0,0.5,1$. We have
rescaled the transition parameter $\beta$ by $L^{2}$. Our computations suggest
that the fully chaotic limit $\gamma=1$ is reached at a critical point
$\beta_{c}(\alpha)=C(\alpha)L^{-2}$, with $C(\alpha)$ a constant that only
depends on the properties of the integrable Liouvillian 12. This points to a
sharp transition from integrability to chaos in the thermodynamical limit as
soon as random Lindblad jumps 15 are added to the system at $\beta>0$. The
transition occurs earlier for $\alpha=0$ probably because there is a wider
distribution of Hamiltonian parameters in the integrable part of the model
with $\alpha=1$. However, the system sizes we are able to reach are not large
enough to make a definite claim in this regard.

Figure 3: Dependence of the fitting parameter $\gamma$ on the transition
parameter $\beta$ rescaled with $L^{2}$ for the three Liouvillian limits:
rational constant $\alpha=0$, general hyperbolic $\alpha=0.5$ and separable
$\alpha=1$, and for different system sizes $L=9,10,11$.

Apart from the interpolation distribution $P(s;\gamma)$ we study the
integrable to chaos transition using the complex spacing ratios recently
introduced by Sá  _et al._ [44]

$z_{k}=\frac{\lambda_{k}^{\textrm{NN}}-\lambda_{k}}{\lambda_{k}^{\textrm{NNN}}-\lambda_{k}}$
(20)

with $\lambda_{k}$ the $k$-th eigenvalue and
$\lambda_{k}^{\textrm{NN}},\lambda_{k}^{\textrm{NNN}}$ its first and second
neighbors, respectively. The advantage of the spacing ratios is that, because
they are adimensional quantities, the unfolding procedure is not needed, at
least if the local density of states does not vary in the scale of the typical
minimal distance between levels. The equivalent quantity for closed systems
introduced in Ref. [45] has become widely used for this reason. The complex
spacing ratios have only been computed for a few examples [31, 44, 46] yet so
it is particularly important to compare their results to the more standard
spacing distributions in the integrable to chaotic transition in Liouvillians.

Figs. 2(b,c) show the means of the absolute and phase angle values,
respectively, of the complex ratios $z_{k}=r_{k}e^{\theta_{k}}$, which were
determined by Sá  _et al._ to be good indicators of the transition. They find
in the limit $L\rightarrow\infty$, the values $(\langle
r\rangle,-\langle\cos\theta\rangle)=(2/3,0)$ for integrable statistics and
$\approx(0.72,0.2)$ for the $AI^{\dagger}$ universality class. Our results for
$\langle r\rangle$ [Fig. 2(b)] coincide with these limits for both integrable
$\beta=0$ and chaotic $\beta=1$ families of Liouvillians, while
$-\langle\cos\theta\rangle$ has a slow convergence to the infinite size limit,
also mentioned by Sá  _et al._ and its limits are a bit lower than those of
$L\rightarrow\infty$. Apart from the finite size effects, both measurements of
the complex spacing rations show a good agreement with the $P(s;\gamma)$
interpolation distribution.

##  IV Conclusions

To summarize, we have presented a new continuous family of integrable many-
body Liouvillians based on the rational and hyperbolic Richardson-Gaudin
models. In vector form, these Liouvillians can be written as a fully coupled
network of spins of the XXX or XXZ type, representing open quantum spin
systems with gain, loss, and dephasing jumps. This very general family of
integrable Liouvillians increases greatly the number of systems that can be
investigated with exact solutions, which can open new avenues of research for
the study of the properties of open quantum many-body systems, a topic that is
still in its infancy.

By adding chaotic jumps, we define two-parameter Liouvillians that describe
the transition between the different integrable models and chaotic non-
integrable Liouvillians with the same degrees of freedom. We have
characterized these transitions studying the spectral statistics of the
Liouvillians’s eigenvalues in the complex plane using two methods, the
nearest-neighbor spacing distribution and the average properties of the
complex spacing ratios. The complex ratios do not require the difficult
unfolding procedure, which is an important advantage as in the case of complex
spectra the unfolding is even more critical than for real spectra.

Within the integrable line, we find excellent agreement between the spacing
distribution and the distribution of a Poisson process in the plane, while in
the chaotic case we find a very good agreement with the $AI^{\dagger}$ random
matrix universality class. The agreement improves as we reach larger sizes. We
have derived a simple interpolation formula between the linear and cubic
repulsion cases as a function of one parameter. The formula captures nicely
the transition between integrable and chaotic dynamics in the spacing
distribution of the complex Liouvillian eigenvalues, allowing a quantitative
description. Our results show that the Liouvillian transitions to chaotic
dynamics with the addition of random Lindblad jumps occurs faster as the size
of the system is increased, pointing to a sharp transition in the
thermodynamic limit as soon as the random jumps are added, although we cannot
reach large enough sizes to make a definite commitment. In the case of the
ratios, the behavior of the average value of their modulus and phase shows
excellent agreement with the behavior of the fitted parameter value of the
interpolation formula as we scan the transition between integrability and
chaos in the parameter space of our model. These results give an extra support
for the use of the complex ratios for characterizing the onset of chaos in
open quantum systems. A similar analysis in other models could help to
elucidate the role of chaos in information scrambling and thermalization in
quantum many body open systems. We anticipate that the applications of these
tools will be as widespread as the equivalent spectral analysis in closed
quantum many body systems.

###### Acknowledgements.

We acknowledge financial support from Projects No. PGC2018-094180-B-I00
(MCIU/AEI/FEDER, EU) and and CAM/FEDER Project No. S2018/TCS-4342 (QUITEMAD-
CM). This research has been also supported by CSIC Research Platform on
Quantum Technologies PTI-001.

## References

  * Srednicki [1994] M. Srednicki, Chaos and quantum thermalization, Phys. Rev. E 50, 888 (1994).
  * Ueda [2020] M. Ueda, Quantum equilibration, thermalization and prethermalization in ultracold atoms, Nature Reviews Physics 2, 669 (2020).
  * Maldacena  _et al._ [2016] J. Maldacena, S. H. Shenker, and D. Stanford, A bound on chaos, Journal of High Energy Physics 2016, 106 (2016).
  * Lewis-Swan  _et al._ [2019] R. J. Lewis-Swan, A. Safavi-Naini, J. J. Bollinger, and A. M. Rey, Unifying scrambling, thermalization and entanglement through measurement of fidelity out-of-time-order correlators in the dicke model, Nature Communications 10, 1581 (2019).
  * Sieberer  _et al._ [2019] L. M. Sieberer, T. Olsacher, A. Elben, M. Heyl, P. Hauke, F. Haake, and P. Zoller, Digital quantum simulation, trotter errors, and quantum chaos of the kicked top, npj Quantum Information 5, 78 (2019).
  * Haake [1991] F. Haake,  _Quantum Coherence in Mesoscopic Systems_ (Springer, 1991) pp. 583–595.
  * Berry  _et al._ [1977] M. V. Berry, M. Tabor, and J. M. Ziman, Level clustering in the regular spectrum, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 356, 375 (1977).
  * Bohigas  _et al._ [1984] O. Bohigas, M. J. Giannoni, and C. Schmit, Characterization of chaotic quantum spectra and universality of level fluctuation laws, Phys. Rev. Lett. 52, 1 (1984).
  * Sieber and Richter [2001] M. Sieber and K. Richter, Correlations between periodic orbits and their role in spectral statistics, Physica Scripta T90, 128 (2001).
  * Müller  _et al._ [2004] S. Müller, S. Heusler, P. Braun, F. Haake, and A. Altland, Semiclassical foundation of universality in quantum chaos, Phys. Rev. Lett. 93, 014103 (2004).
  * Müller  _et al._ [2009] S. Müller, S. Heusler, A. Altland, P. Braun, and F. Haake, Periodic-orbit theory of universal level correlations in quantum chaos, New Journal of Physics 11, 103025 (2009).
  * Stöckmann [1999] H.-J. Stöckmann,  _Quantum chaos: an introduction_ (Cambridge University Press, Cambridge, 1999).
  * Gutzwiller [2007] M. Gutzwiller, Quantum chaos, Scholarpedia 2, 3146 (2007).
  * Gómez  _et al._ [2011] J. Gómez, K. Kar, V. Kota, R. Molina, A. Relaño, and J. Retamosa, Many-body quantum chaos: Recent developments and applications to nuclei, Physics Reports 499, 103 (2011).
  * Geiger  _et al._ [2017] B. Geiger, J.-D. Urbina, Q. Hummel, and K. Richter, Semiclassics in a system without classical limit: The few-body spectrum of two interacting bosons in one dimension, Phys. Rev. E 96, 022204 (2017).
  * Bertini  _et al._ [2018] B. Bertini, P. Kos, and T. c. v. Prosen, Exact spectral form factor in a minimal model of many-body quantum chaos, Phys. Rev. Lett. 121, 264101 (2018).
  * Kos  _et al._ [2018] P. Kos, M. Ljubotina, and T. c. v. Prosen, Many-body quantum chaos: Analytic connection to random matrix theory, Phys. Rev. X 8, 021062 (2018).
  * Brody [1973] T. A. Brody, A statistical measure for the repulsion of energy levels, Lettere al Nuovo Cimento (1971-1985) 7, 482 (1973).
  * Berry and Robnik [1984] M. V. Berry and M. Robnik, Semiclassical level spacings when regular and chaotic orbits coexist, Journal of Physics A: Mathematical and General 17, 2413 (1984).
  * Prosen and Robnik [1994] T. Prosen and M. Robnik, Semiclassical energy level statistics in the transition region between integrability and chaos: transition from brody-like to berry-robnik behaviour, Journal of Physics A: Mathematical and General 27, 8059 (1994).
  * Grobe  _et al._ [1988] R. Grobe, F. Haake, and H.-J. Sommers, Quantum distinction of regular and chaotic dissipative motion, Phys. Rev. Lett. 61, 1899 (1988).
  * Breuer and Petruccione [2007] H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, Oxford, 2007).
  * Ginibre [1965] J. Ginibre, Statistical ensembles of complex, quaternion, and real matrices, Journal of Mathematical Physics 6, 440 (1965), https://doi.org/10.1063/1.1704292 .
  * Jaiswal  _et al._ [2019] A. B. Jaiswal, A. Pandey, and R. Prakash, Universality classes of quantum chaotic dissipative systems, EPL (Europhysics Letters) 127, 30004 (2019).
  * Hamazaki  _et al._ [2020] R. Hamazaki, K. Kawabata, N. Kura, and M. Ueda, Universality classes of non-hermitian random matrices, Phys. Rev. Research 2, 023286 (2020).
  * Akemann  _et al._ [2019] G. Akemann, M. Kieburg, A. Mielke, and T. Prosen, Universal signature from integrability to chaos in dissipative open quantum systems, Phys. Rev. Lett. 123, 254101 (2019).
  * Medvedyeva  _et al._ [2016] M. V. Medvedyeva, F. H. L. Essler, and T. Prosen, Exact bethe ansatz spectrum of a tight-binding chain with dephasing noise, Phys. Rev. Lett. 117, 137202 (2016).
  * Rowlands and Lamacraft [2018] D. A. Rowlands and A. Lamacraft, Noisy spins and the richardson-gaudin model, Phys. Rev. Lett. 120, 090401 (2018).
  * Ribeiro and Prosen [2019] P. Ribeiro and T. Prosen, Integrable quantum dynamics of open collective spin models, Phys. Rev. Lett. 122, 010401 (2019).
  * Denisov  _et al._ [2019] S. Denisov, T. Laptyeva, W. Tarnowski, D. Chruściński, and K. Życzkowski, Universal spectra of random lindblad operators, Phys. Rev. Lett. 123, 140403 (2019).
  * Sá  _et al._ [2020] L. Sá, P. Ribeiro, and T. Prosen, Spectral and steady-state properties of random liouvillians, Journal of Physics A: Mathematical and Theoretical 53, 305303 (2020).
  * van Caspel and Gritsev [2018] M. van Caspel and V. Gritsev, Symmetry-protected coherent relaxation of open quantum systems, Phys. Rev. A 97, 052106 (2018).
  * 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).
  * Lee  _et al._ [2014] T. E. Lee, C.-K. Chan, and S. F. Yelin, Dissipative phase transitions: Independent versus collective decay and spin squeezing, Phys. Rev. A 90, 052109 (2014).
  * Dhar  _et al._ [2018] H. S. Dhar, M. Zens, D. O. Krimer, and S. Rotter, Variational renormalization group for dissipative spin-cavity systems: Periodic pulses of nonclassical photons from mesoscopic spin ensembles, Phys. Rev. Lett. 121, 133601 (2018).
  * Buča and Jaksch [2019] B. Buča and D. Jaksch, Dissipation induced nonstationarity in a quantum gas, Phys. Rev. Lett. 123, 260401 (2019).
  * Can  _et al._ [2019] T. Can, V. Oganesyan, D. Orgad, and S. Gopalakrishnan, Spectral gaps and midgap states in random quantum master equations, Phys. Rev. Lett. 123, 234103 (2019).
  * Dukelsky  _et al._ [2001] J. Dukelsky, C. Esebbag, and P. Schuck, Class of exactly solvable pairing models, Phys. Rev. Lett. 87, 066403 (2001).
  * Ortiz  _et al._ [2005] G. Ortiz, R. Somma, J. Dukelsky, and S. Rombouts, Exactly-solvable models derived from a generalized gaudin algebra, Nuclear Physics B 707, 421 (2005).
  * Dukelsky  _et al._ [2004] J. Dukelsky, S. Pittel, and G. Sierra, Colloquium: Exactly solvable richardson-gaudin models for many-body quantum systems, Rev. Mod. Phys. 76, 643 (2004).
  * Richardson and Sherman [1964] R. Richardson and N. Sherman, Exact eigenstates of the pairing-force hamiltonian, Nuclear Physics 52, 221 (1964).
  * Stouten  _et al._ [2019] E. Stouten, P. W. Claeys, J.-S. Caux, and V. Gritsev, Integrability and duality in spin chains, Phys. Rev. B 99, 075111 (2019).
  * Gómez  _et al._ [2002] J. M. G. Gómez, R. A. Molina, A. Relaño, and J. Retamosa, Misleading signatures of quantum chaos, Phys. Rev. E 66, 036209 (2002).
  * Sá  _et al._ [2020] L. Sá, P. Ribeiro, and T. Prosen, Complex spacing ratios: A signature of dissipative quantum chaos, Phys. Rev. X 10, 021019 (2020).
  * Oganesyan and Huse [2007] V. Oganesyan and D. A. Huse, Localization of interacting fermions at high temperature, Phys. Rev. B 75, 155111 (2007).
  * Peron  _et al._ [2020] T. Peron, B. M. F. de Resende, F. A. Rodrigues, L. d. F. Costa, and J. A. Méndez-Bermúdez, Spacing ratio characterization of the spectra of directed random networks, Phys. Rev. E 102, 062305 (2020).