aainstitutetext: Department of Theoretical Physics, Tata Institute of
Fundamental Research, 1 Homi Bhabha Road, Mumbai 400005, India

# Complexity Growth and the Krylov-Wigner function

Ritam<EMAIL_ADDRESS>a,2 Anirban Ganguly,222Visiting student
at TIFR<EMAIL_ADDRESS>a,3 Souparna
<EMAIL_ADDRESS>a,4 Onkar
<EMAIL_ADDRESS>
###### Abstract

For any state in a $D$-dimensional Hilbert space with a choice of basis, one
can define a discrete version of the Wigner function – a quasi-probability
distribution which represents the state on a discrete phase space. The Wigner
function can, in general, take on negative values, and the amount of
negativity in the Wigner function has an operational meaning as a resource for
quantum computation. In this note, we study the growth of Wigner negativity
for a generic initial state under time evolution with chaotic Hamiltonians. We
introduce the Krylov-Wigner function, i.e., the Wigner function defined with
respect to the Krylov basis (with appropriate phases), and show that this
choice of basis minimizes the early time growth of Wigner negativity in the
large $D$ limit. We take this as evidence that the Krylov basis (with
appropriate phases) is ideally suited for a dual, semi-classical description
of chaotic quantum dynamics at large $D$. We also numerically study the time
evolution of the Krylov-Wigner function and its negativity in random matrix
theory for an initial pure state. We observe that the negativity broadly shows
three phases: it rises gradually for a time of $O(\sqrt{D})$, then hits a
sharp ramp and finally saturates close to its upper bound of $\sqrt{D}$.

##  1 Introduction

In general, given a classical phase space, there is a more-or-less standard
procedure one can follow to quantize the theory. But given a quantum theory,
it is not always clear whether there exists an underlying classical
description which is “useful”. For instance, in the AdS/CFT correspondence,
the large $N$ boundary conformal field theory is a quantum mechanical system
at strong coupling, so the underlying classical description in terms of
classical boundary fields is far from being useful. However, it is a highly
non-trivial fact that a new useful classical description emerges in the form
of the bulk gravitational theory, with the bulk $\hbar\sim\frac{1}{N^{2}}$.
Given a strongly coupled quantum system, it would be nice to have a general
strategy to identify when it admits a dual, semi-classical description,
particularly one involving gravity.

A related question has been studied in quantum information theory – given a
quantum circuit, does it admit a “useful” classical description? A quantum
circuit consists of some initial state which is sequentially acted upon by a
set of quantum gates/operations, followed by a measurement of some observable.
The output of the quantum circuit (over several runs) is a probability
distribution for the observable of interest. We say that the quantum circuit
admits a useful classical description if the same probability distribution can
be sampled by an efficient classical algorithm. Here, by efficient we mean
that the running time of the classical algorithm should scale polynomially
with the spatial size and depth of the quantum circuit.

A partial answer to this question is given by the Gottesman-Knill theorem
gottesman1998heisenberg ; Aaronson:2004xuh  and its generalizations Mari_2012
; Veitch_2012 ; Pashayan_2015 ; Wang_2019 . In order to state the result, we
need to introduce the notion of the _Wigner function_ Wigner  (see also
Dhar:1992rs ; Dhar:1992hr ; Lin:2004nb ; Balasubramanian:2005mg  for some
applications of the Wigner function formalism in holography). The Wigner
function $W_{\psi}(q,p)$ for a state $\psi$ is a quasi-probability
distribution, which gives a representation of the quantum state in phase
space, parametrized here by $(q,p)$. This should be contrasted against the
wavefunction $\psi(q)$ of the state, whose absolute square gives a probability
distribution over a half-dimensional (Lagrangian) submanifold, parametrized
here by $q$. These ideas can also be extended to a quantum system with a
finite $D$-dimensional Hilbert space WOOTTERS19871 ; Leonard ; sphere ;
quantumcomp ; Galois ; Gross_2006 ; classicality ; in this case, one picks an
orthonormal basis for the Hilbert space and corresponding to this choice, the
Wigner function is defined on a discrete phase space lattice
$\mathbb{Z}_{D}\times\mathbb{Z}_{D}$.555The discrete phase space construction
works best when $D$ is a prime number, since in this case $\mathbb{Z}_{D}$ has
the structure of a field. The formalism can also be generalized to include
powers of primes. Importantly, the Wigner function is not a true probability
distribution because it can take on negative values. Nevertheless, it is a
natural quantum analog of a classical probability distribution in phase space;
for instance, it is real and unit normalized, and under time evolution follows
a quantum generalization of the classical Liouville’s equation for the time
evolution of phase space densities. When the Wigner function for a state is
everywhere non-negative, it acts as a genuine probability distribution, and
one might roughly expect such states to be “classical”. A Wigner function can
be associated not only to states, but also to unitary quantum gates (more
generally quantum channels) and to measurements. This brings us back to the
original question – when does a quantum circuit admit a useful classical
description? It turns out that any “local” quantum circuit which uses elements
(initial state, quantum gates, measurements etc.) with positive Wigner
functions admits an efficient classical implementation Mari_2012 ; Veitch_2012
; Pashayan_2015 ; Wang_2019 . Thus, the question of whether one can find an
equivalent classical description for the quantum circuit translates to whether
one can find a basis in terms of which the elements of a quantum circuit admit
positive Wigner function representations.

On the other hand, given a choice of the basis, classically simulable circuits
such as above evolve within a very small subset of states, called _stabilizer
states_. Stabilizer states and circuits are not universal – in order to
achieve universality, one needs non-stabilizer elements in the circuit. The
negativity of the Wigner function turns out to be an operationally meaningful
measure (i.e., a monotone under stabilizer operations) of the non-stabilizer
content of a state. Much like entanglement plays the role of a resource in the
theory of quantum communication, the negativity of the Wigner function plays
the role of a resource in the theory of quantum computation using stabilizer
circuits Veitch_2012 ; Veitch_2014 . Indeed, Wigner negativity (or “non-
stabilizerness”) is an important feature of various interesting physical
systems including conformal field theories White:2020zoz , condensed matter
systems Liu:2020yso ; Turkeshi:2023lqu ; Niroula:2023meg , topological field
theories Fliss:2020yrd , and chaotic quantum systems Leone2021quantumchaosis ;
OLIVIERO2021127721 ; Goto .

In this note, our aim is to study the time evolution of the negativity of the
Wigner function for generic initial states under chaotic Hamiltonian
evolution. In this context, we introduce the _Krylov-Wigner function_ – the
Wigner function with respect to a specific choice of basis called the Krylov
basis Balasubramanian:2022tpr .666See also Parker:2018yvk ; Jian:2020qpp ;
Rabinovici_2021 ; caputa2021geometry  for some previous work related to Krylov
basis in the context of operator spreading. Given an initial state $\psi_{0}$
and a chaotic Hamiltonian $H$, the Krylov basis is generated by repeatedly
acting with $H$ on $\psi_{0}$,

$\psi_{0},\,H\psi_{0},\;H^{2}\psi_{0}\cdots$

and then orthonormalizing using the Gram-Schmidt process:

$\displaystyle|0_{\mathcal{K}}\rangle$ $\displaystyle=$
$\displaystyle|\psi_{0}\rangle,$ $\displaystyle|1_{\mathcal{K}}\rangle$
$\displaystyle=$
$\displaystyle\frac{1}{\sqrt{N_{1}}}\left(H|\psi_{0}\rangle-\langle
0_{\mathcal{K}}|H|\psi_{0}\rangle\,|0_{\mathcal{K}}\rangle\right),$
$\displaystyle|2_{\mathcal{K}}\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{N_{2}}}\left(H^{2}|\psi_{0}\rangle-\langle
0_{\mathcal{K}}|H^{2}|\psi_{0}\rangle\,|0_{\mathcal{K}}\rangle-\langle
1_{\mathcal{K}}|H^{2}|\psi_{0}\rangle|1_{\mathcal{K}}\rangle\right),$
$\displaystyle\vdots$ (1)

where $N_{i}$ are normalization factors. The aspect of chaos that will be
important for us is ergodicity, i.e., that the time evolution of a generic
initial state should cover the entire (or a large subspace of the) Hilbert
space. One of our main results is that the Krylov basis (up to multiplication
by individual phases) minimizes the early time growth of the Wigner negativity
in the large-$D$ limit. We view this as evidence that the Krylov basis is a
potentially good basis for finding an efficient “close to classical”
description for the quantum time evolution of the initial state.777However,
note that we have not given an algorithm to explicitly demonstrate that the
time evolution of the state can be efficiently modeled on a classical computer
in this basis. Our result builds on the work of Balasubramanian:2022tpr ,
where it was shown that the Krylov basis minimizes the amount of early time
spreading of the wavefunction under time evolution, a quantity which was
referred to as the “spread complexity”. This already suggests that the
wavefunction for the time-evolving state is as peaked as possible in the
Krylov basis. One of the motivations of this paper is to show that there is a
sense in which the Krylov basis makes the time-evolving wavefunction as
classical as possible, and to further explore connections between the work of
Balasubramanian:2022tpr  and computational complexity, say along the lines of
nielsen2005geometric ; Nielsen_2006 ; Jefferson_2017 ; Chapman_2018 ;
Brown_2018 ; Balasubramanian:2019wgd ; Balasubramanian:2021mxo  (see also
Craps:2023ivc  for recent work in this direction).

Incidentally, the Krylov basis makes an interesting appearance in a certain
model of two dimensional quantum gravity called the double-scaled SYK (DSSYK)
model Cotler:2016fpe ; Berkooz:2018jqr ; Lin:2022rbf . The DSSYK model is a
model of $N$ Majorana fermions with all to all $p$ body interactions where the
couplings are drawn from a Gaussian distribution, and where both $p$ and $N$
are taken to infinity with $\lambda=\frac{2p^{2}}{N}$ fixed. The model is of
interest to quantum gravity because upon averaging over the couplings, the
“chord-diagram” expansion for the observables in this model, in a certain
triple scaling limit, turns into a calculation of observables in Jackiw-
Teitelboim (JT) gravity maldacena2016conformal . In this correspondence, the
Krylov basis in the DSSYK model turns into the “chord-length” basis in the
dual description Lin:2022rbf ; Rabinovici:2023yex . In the JT gravity limit,
the chord length precisely corresponds to the length of the bulk wormhole, and
furthermore, the growth of the spread complexity of Balasubramanian:2022tpr
corresponds to the growth of the length (what would have been volume in higher
dimensions) of the wormhole interior, in line with the “Complexity = Volume”
conjecture susskind2014computational ; Susskind:2014moa ; Stanford:2014jda .
Thus, at least in these models, the dual gravitational description seems to
utilise the Krylov basis in a natural way. That the Krylov basis minimizes the
early growth of the Wigner negativity in the large-$D$ limit suggests a
possibly general underlying principle: _the gravitational variables correspond
to a basis which makes the time evolution as classical as possible_.

In order to further explore this proposal for “gravity as efficient classical
implementation of quantum computation”,888See also Czech:2017ryf ;
Caputa:2018kdj  for similar/related ideas. it would be helpful to study Wigner
negativity and its growth in models such as DSSYK or the Saad-Shenker-Stanford
matrix model Saad:2019lba . A significant obstacle is that the Wigner function
and its negativity are difficult to compute analytically, but as a first step
towards this goal, here we numerically study the time evolution of the Krylov-
Wigner function and its negativity in Gaussian random matrix theory (GUE). We
observe that the negativity broadly shows three phases (see figure 1):
initially, the KW function spreads locally over phase space and its negativity
shows a gradual rise. At a time $t\sim\sqrt{D}$, as the Wigner function starts
to cover the entire phase space, the negativity hits a sharp ramp. After this,
the Wigner negativity hits a plateau and saturates close to its upper bound of
$\sqrt{D}$. On the other hand, in a generic basis, the Wigner negativity
encounters a sharp ramp at $t=0$ and saturates close to its upper bound value
within an $O(1)$ amount of time evolution (see figure 3). These results are
consistent with the claim that the Krylov basis minimizes the early-time
growth of Wigner negativity, and further suggest that this basis is useful for
times not scaling exponentially in $\log\,D$.

The plan for the paper is as follows: in section 2 we review some background
material about discrete Wigner functions and the interpretation of Wigner
negativity as a resource for quantum computation. In section 3, we introduce
the Krylov-Wigner function, i.e., the Wigner function defined with respect to
the Krylov basis, and show that this choice of computational basis minimizes
the early-time growth of the negativity in the large $D$ limit. In section 4,
we present some numerical results on the growth of the Krylov-Wigner
negativity in random matrix theory. Finally, we end with some discussion on
loose ends, open questions and future directions.

##  2 Preliminaries

###  2.1 The Wigner function

Consider a particle moving along on a line under some potential.
Traditionally, the classical phase space is the two-dimensional plane
coordinatized by the position of the particle $q$ and the momentum $p$. The
phase space comes equipped with a symplectic structure, i.e., a closed, non-
degenerate 2-form $\omega$, which in the above coordinates is given by

$\omega=dq\wedge dp.$ (2)

Equivalently, we can use the above symplectic form to define Poisson brackets
on the space of functions on phase space:

$\left\\{f,g\right\\}_{PB}=\omega^{IJ}\partial_{I}f\partial_{J}g,$ (3)

where $I=1,2$ runs over phase space indices and $\omega^{IJ}$ is the inverse
matrix corresponding to the 2-form $\omega$ thought of as a matrix. Time
evolution in phase space is governed by Hamilton’s equations, and
correspondingly, given a probability density $u(q,p)$ in phase space, its time
evolution is governed by the Louiville equation:

$\frac{du}{dt}=\left\\{H,u\right\\}_{PB}.$ (4)

Quantum mechanically, the particle is described by a Hilbert space. Functions
on phase space get promoted to operators on Hilbert space, while the Poisson
bracket gets promoted to the commutator. Given a Hamiltonian operator, the
time evolution of a density matrix on the Hilbert space is given by

$\frac{d\rho}{dt}=\frac{i}{\hbar}\left[H,\rho\right].$ (5)

Given a quantum theory, how does one recover the underlying classical phase
space? The Wigner function approach is one way to proceed, where one defines a
quasi-probability distribution called the Wigner function:

$W(q,p)=\frac{1}{2\pi\hbar}\int_{-\infty}^{\infty}dy\,\matrixelement{q-\frac{y}{2}}{\rho}{q+\frac{y}{2}}\exp(-\frac{i}{\hbar}py).$
(6)

Equivalently, one can write the Wigner function as

$W(q,p)=\frac{1}{2\pi\hbar}\text{Tr}\left(\rho\,A(q,p)\right),$ (7)

where the operators $A(q,p)$ are called _phase point operators_ :

$\langle
q^{\prime}|A(q,p)|q^{\prime\prime}\rangle=\delta(q-\frac{(q^{\prime}+q^{\prime\prime})}{2})e^{-\frac{ip}{\hbar}(q^{\prime}-q^{\prime\prime})}.$
(8)

The phase-point operators are naturally related to the _Heisenberg-Weyl_ (HW)
displacement operators:

$w(q,p)=e^{\frac{i}{\hbar}(p\,\mathbf{q}-q\,\mathbf{p})}=e^{-\frac{i}{2\hbar}pq}Z(p)X(q),$
(9)
$Z(p)=e^{\frac{i}{\hbar}p\,\mathbf{q}},\;\;X(q)=e^{-\frac{i}{\hbar}q\,\mathbf{p}}$

where $\mathbf{q}$ and $\mathbf{p}$ are canonically conjugate operators which
satisfy $\left[\mathbf{q},\mathbf{p}\right]=i\hbar$. The phase point operators
are symplectic Fourier transforms of the HW operators:

$A(q,p)=\int\frac{dp^{\prime}dq^{\prime}}{2\pi\hbar}e^{\frac{i}{\hbar}(p^{\prime}q-q^{\prime}p)}w(q^{\prime},p^{\prime}),$
(10)

and can be thought of, heuristically, as being delta-function localized at the
phase-space point $(q,p)$.

One can think of the Wigner function as a representation of the quantum state
in classical phase space, but note that the Wigner function is, generally
speaking, not everywhere positive, so it is not a classical probability
distribution. Some useful properties of the Wigner function are:

  1. 1.

The Wigner function is real and unit normalized:

$\int_{-\infty}^{\infty}dx\int_{-\infty}^{\infty}dp\,W(x,p)=1.$ (11)

  2. 2.

Integrating the Wigner function along one-direction gives the probability
density. For e.g.,

$\int_{-\infty}^{\infty}dp\,W(x,p)=\langle x|\rho|x\rangle.$ (12)

  3. 3.

The time evolution of the Wigner function is controlled by the Moyal equation,
which reduces to the Louiville equation in the $\hbar\to 0$ limit:

$\displaystyle\frac{dW}{dt}$ $\displaystyle=$
$\displaystyle\frac{2}{\hbar}\sin\left(\frac{\hbar}{2}(\partial_{x^{\prime}}\partial_{p}-\partial_{p^{\prime}}\partial_{x})\right)h(x^{\prime},p^{\prime})W(x,p)\Big{|}_{x^{\prime}=x,p^{\prime}=p}$
(13) $\displaystyle=$ $\displaystyle\left\\{h,W\right\\}_{PB}+O(\hbar),$

where $h(x,p)$ is defined as in equation (6), but with $\rho$ replaced with
the Hamiltonian $H$.

As mentioned above, the Wigner function is not always everywhere positive, and
so it is not quite a probability distribution. States for which the Wigner
function is positive may be regarded as _classical_ states Anatole . Hudson’s
theorem HUDSON1974249 ; Soto  shows that the only states for which the Wigner
function is everywhere positive are _Gaussian_ states:

$\psi(q)\sim e^{-\frac{1}{2}aq^{2}-bq},$ (14)

where $a,b\in\mathbb{C}$.999When $a=1$ and $b=0$, this is simply the harmonic-
oscillator vacuum, which has a angularly symmetric Wigner function centered at
the origin. Turning on an imaginary part for $a$ causes the Wigner function to
get squeezed in one direction. The real and imaginary parts of $b$ displace
the center of the Wigner function. Thus, states with everywhere positive
Wigner functions form a very restricted class of states in the Hilbert space.
But even if we were to start with a classical/Gaussian state, the positivity
of the Wigner function will not, in general, be preserved under time
evolution. This only happens for at most inhomogenous, quadratic (i.e., free)
Hamiltonians Soto2 . The space of such at most quadratic Hamiltonians
generates a group, which we may call the “Clifford” group. The Clifford group
is the subgroup of the unitary group which normalizes the HW group, or in
other words, maps HW operators to HW operators. Any quantum evolution which
starts from a Gaussian state and involves only Clifford gates will keep the
state within the restricted class of Gaussian states. This leads to a massive
reduction in the complexity of the problem: instead of having to solve for the
entire wavefunction $\psi(t,q)$ for every $q$ as a function of time, we now
only need to keep track of the finite set of parameters $(a(t),b(t))$. In this
sense, quantum dynamics gets mapped to generalized classical dynamics for such
systems, and time evolution can be efficiently simulated on a classical
computer Bartlett .

###  2.2 The discrete Wigner function

So far, the discussion has been about infinite dimensional Hilbert spaces,
where the phase space is non-compact. Our main focus in this paper will be on
quantum systems with finite dimensional Hilbert spaces. For finite dimensional
quantum systems, there is a natural analog of the Wigner function, and is
referred to as the _discrete_ Wigner function WOOTTERS19871 ; Leonard ; sphere
; quantumcomp ; Galois ; Gross_2006 ; classicality . In the finite dimensional
case, one defines a _discrete phase space_ $\mathcal{P}$ as a lattice of size
$D\times D$, where $D$ is the dimension of the Hilbert space. The discrete
Wigner function formalism works best when $D$ is a prime number.101010For $D$
either prime or a power of a prime, there exists a field of order $D$. When
$D$ is prime, this field is simply $\mathbb{Z}_{D}$. For simplicity, we will
take $D$ to be a prime number in this work.111111For a Hilbert space with non-
prime dimension, we can always embed it in a larger Hilbert space with prime
dimension to construct Wigner function representations. In this case, we think
of the discrete phase space $\mathcal{P}$ as
$\mathbb{Z}_{D}\times\mathbb{Z}_{D}$, with a _symplectic inner product_
between two points $\vec{\alpha}=(q,p)$ and
$\vec{\beta}=(q^{\prime},p^{\prime})$ given by:

$\left[\vec{\alpha},\vec{\beta}\right]=pq^{\prime}-qp^{\prime}.$ (15)

Let $\left\\{|k\rangle\right\\}_{k=0}^{D-1}$ be a particular orthonormal basis
for the Hilbert space, which is chosen to begin with. In terms of this basis,
one defines a discrete version of the Heisenberg-Weyl (HW) operators:

$w(q,p)=e^{-\frac{2\pi i}{D}\frac{D+1}{2}qp}Z(p)X(q),$ (16)

where $(q,p)\in\mathcal{P}$ and

$Z(p)|q\rangle=e^{\frac{2\pi
ipq}{D}}|q\rangle,\;\;\;X(q)|q^{\prime}\rangle=|q^{\prime}+q\rangle.$ (17)

HW operators are closed under multiplication:

$w(q,p)w(q^{\prime},p^{\prime})=e^{\frac{2\pi
i}{D}\frac{(pq^{\prime}-qp^{\prime})}{2}}w(q+q^{\prime},p+p^{\prime}).$ (18)

In terms of the HW operators, the discrete phase-point operators are defined
by a discrete version of the symplectic Fourier transform WOOTTERS19871 :

$A(q,p)=\sum_{k,\ell=0}^{D-1}\widehat{\delta}_{2q,k+l}e^{\frac{2\pi
i}{D}(k-\ell)p}\outerproduct{k}{\ell}.$ (19)

Note that the hatted Kronecker delta $\widehat{\delta}$ is the $\text{mod}\,D$
version, i.e., it is one when $(k+\ell)=2q\,\text{mod}\,D$, and zero
otherwise. We will often use the notation $\vec{\alpha}=(q,p)$ to denote a
phase-space point, and correspondingly denote the phase-point operator as
$A_{\vec{\alpha}}$. These phase-point operators satisfy the following
properties:

$\begin{split}&\Tr(A_{\vec{\alpha}})=1,\\\
&\Tr(A_{\vec{\alpha}}A_{\vec{\beta}})=D\,\delta_{\vec{\alpha},\vec{\beta}},\\\
&\frac{1}{D}\sum_{\vec{\alpha}}A_{\vec{\alpha}}=\mathbf{1}.\end{split}$ (20)

The phase-point operators defined above are the finite-dimensional analogs of
the phase-point operators in the continuous case (see equation (8)). In
analogy with the continuous case, the discrete Wigner function for a density
matrix $\rho$ is now defined in terms of discrete phase point operators as:

$W_{\rho}(\vec{\alpha})=\frac{1}{D}\mathrm{Tr}\left(\rho
A_{\vec{\alpha}}\right).$ (21)

Equivalently, any density matrix can be written in terms of its Wigner
function as $\rho=\sum_{\vec{\alpha}}W_{\vec{\alpha}}A_{\vec{\alpha}}$. The
discrete Wigner function satisfies the following properties:

  1. 1.

The discrete Wigner function is real and is unit normalized, i.e.,

$\Tr(\rho)=1\implies\sum_{{\vec{\alpha}}}W_{\rho}(\vec{\alpha})=1.$ (22)

  2. 2.

Summing over one of the directions, say $p$, reduces the Wigner function to
the probability density along the other direction:

$\sum_{p=0}^{D-1}W_{\rho}(x,p)=\langle x|\rho|x\rangle.$ (23)

  3. 3.

The time evolution of the discrete Wigner function follows a discrete version
of the Moyal equation:

$\frac{dW_{\rho}(\vec{\alpha})}{dt}=\frac{2}{D}\sum_{\vec{\beta},\vec{\gamma}}\sin\left(\frac{4\pi}{D}\mathcal{A}_{\vec{\alpha}\vec{\beta}\vec{\gamma}}\right)H(\vec{\beta})W_{\rho}(\vec{\gamma}),$
(24)

where

$\mathcal{A}_{\vec{\alpha}\vec{\beta}\vec{\gamma}}=\alpha_{2}(\gamma_{1}-\beta_{1})+\beta_{2}(\alpha_{1}-\gamma_{1})+\gamma_{2}(\beta_{1}-\alpha_{1}),$
(25)

and the Hamiltonian is expressed in the basis of the phase point operators as
$H=\displaystyle\sum_{\vec{\alpha}}H(\vec{\alpha})A_{\vec{\alpha}}$. To see
that the discrete Wigner function satisfies the above differential equation,
recall that the density matrix $\rho(t)$ follows the evolution equation (we
will henceforth set $\hbar=1$):

$\frac{\mathrm{d}\rho}{\mathrm{d}t}=i[\rho,H].$ (26)

Substituting the density matrix in terms of its Wigner function in here we get

$\frac{\mathrm{d}W_{\rho}(\vec{\alpha})}{\mathrm{d}t}=-\frac{i}{D}\sum_{\vec{\beta},\vec{\gamma}}\mathrm{Tr}\left(A_{\vec{\alpha}}\left[A_{\vec{\beta}},A_{\vec{\gamma}}\right]\right)H(\vec{\beta})W_{\rho}(\vec{\gamma}).$
(27)

This determines the time evolution of the Wigner function in the discrete
phase space. Using the explicit form of the phase point operators given in
equation (19) above, we then land on equation (24) WOOTTERS19871 .

So far we have discussed Wigner functions for states, but Wigner functions can
also be assigned to quantum gates, and more generally to quantum
operations/channels Mari_2012 ; Wang_2019 . For a quantum channel
$\mathcal{E}$, we consider the corresponding Choi-Jamiolkowski state:

$\sigma_{\mathcal{E}}=(\mathbf{1}\otimes\mathcal{E})(|\Omega\rangle\langle\Omega|),$
(28)

where

$|\Omega\rangle=\sum_{k=1}^{D}|k\rangle\otimes|k\rangle,$ (29)

is the maximally entangled state on two copies of the Hilbert space. Then, the
Wigner function for the quantum operation $\mathcal{E}$ is defined as

$W_{\mathcal{E}}(\vec{\beta}|\vec{\alpha})=\frac{1}{D}\mathrm{Tr}\left(A^{T}_{\vec{\alpha}}\otimes
A_{\vec{\beta}}\,\sigma_{\mathcal{E}}\right)=\frac{1}{D}\mathrm{Tr}\left(A_{\vec{\beta}}\,\mathcal{E}(A_{\vec{\alpha}})\right).$
(30)

Given an input state $\rho$ with the Wigner function $W_{\rho}(\vec{\alpha})$,
the Wigner function of the output state $\mathcal{E}(\rho)$ is given by

$W_{\mathcal{E}(\rho)}(\vec{\beta})=\sum_{\vec{\alpha}}W_{\mathcal{E}}(\vec{\beta}|\vec{\alpha})\,W_{\rho}(\vec{\alpha}).$
(31)

Similarly, for a positive operator valued measure $\\{M_{k}\\}$, the
corresponding Wigner function is defined as

$W_{k}(\vec{\alpha})=\mathrm{Tr}\left(M_{k}\,A_{\vec{\alpha}}\right).$ (32)

In this way, Wigner functions can be assigned to general elements of a quantum
circuit.

###  2.3 Negativity of the Wigner function

As in the continuous case, the discrete Wigner function is also a quasi-
probability distribution, in that it can take on negative values. It is
natural to regard the set of states with non-negative Wigner functions as
being _classical_ , but interestingly, the sense in which such states are
classical has a computational flavor. Given the choice of an orthonormal
basis, states with a non-negative Wigner function are again known to be
Gaussian Gross_2006 :

$\langle q|\psi\rangle\sim e^{\frac{2\pi i}{D}(aq^{2}+bq)},$ (33)

where $a,b\in\mathbb{Z}_{D}$. The corresponding discrete Wigner functions are
highly constrained. We say that a subspace $\mathcal{L}$ of phase space is an
_isotropic subspace_ if for any $\vec{\alpha},\vec{\beta}\in\mathcal{L}$, we
have

$\left[\vec{\alpha},\vec{\beta}\right]=0.$ (34)

It follows from equation (18) that the HW group elements corresponding to an
isotropic subspace form a subgroup. We say that $\mathcal{L}$ is maximally
isotropic if its cardinality is $D$. A maximally isotropic subspace is the
discrete analog of a Lagrangian submanifold in phase space. It turns out that
positivity is a rather restrictive criterion in the discrete case, and implies
that the Wigner function has to be of the form Gross_2006 :

$W_{\psi}(q,p)=\begin{cases}&\frac{1}{D}\;\;\cdots\;\;\text{if}\;(q,p)\in\mathcal{L}\\\
&0\;\;\cdots\;\;\text{if}\;(q,p)\notin\mathcal{L},\end{cases}$ (35)

where $\mathcal{L}$ is some maximally isotropic subspace. Thus, states with
non-negative Wigner functions are “translationally invariant” along some
“line” in phase space. For this reason, such states are called _stabilizer
states_ , i.e., they are stabilized by the maximally isotropic subgroup of the
HW group corresponding to $\mathcal{L}$.

Note that even though a stabilizer state has a positive Wigner function, in
general, unitary evolution will introduce negativity. However, there is a
special subgroup of unitaries which preserve the positivity of the Wigner
function. This is called the _Clifford group_ , and is defined as the subset
of unitary operators which map HW operators to HW operators:

$Uw(\vec{\alpha})U^{\dagger}=w(S(\vec{\alpha})),$ (36)

for some $S:\mathcal{P}\to\mathcal{P}$. In other words, the Clifford group is
the normalizer of the HW group. There is a very elegant characterization of
the Clifford group in terms of symplectic affine transformations on phase
space Gross_2006 . A _symplectic_ transformation $S:\mathcal{P}\to\mathcal{P}$
on phase space is a linear map (with coefficients in $\mathbb{Z}_{D}$) which
preserves the symplectic inner product:

$\left[S(\vec{\alpha}),S(\vec{\beta})\right]=\left[\vec{\alpha},\vec{\beta}\right].$
(37)

Corresponding to every symplectic map, there exists a unitary operator
$\mu(S)$ such that

$\mu(S)\,w(\vec{\alpha})\,\mu(S)^{\dagger}=w(S(\vec{\alpha})).$ (38)

Up to a phase, every Clifford unitary is of the form
$U=w(\vec{\alpha})\,\mu(S)$. Thus, the Clifford group is essentially a
representation of the group of symplectic affine transformations on phase
space. The Wigner function transforms covariantly under the Clifford group,
namely

$W_{U\psi}(\vec{\beta})=W_{\psi}(S(\vec{\beta})+\vec{\alpha}).$ (39)

Since the Clifford group normalizes the Heisenberg-Weyl group, Clifford
rotations preserve the non-negativity of the Wigner function, as can be seen
from equation (39). Consequently, Clifford rotations map stabilizer states to
stabilizer states. Thus, any quantum circuit which is initialized to a
stabilizer state and involves only Clifford operations evolves within a very
restricted subset of the Hilbert space. In fact, it is a well-known fact that
such quantum circuits can be simulated efficiently on a classical computer.
This is the content of the _Gottesman-Knill theorem_ gottesman1998heisenberg ;
Aaronson:2004xuh .

More generally, any circuit which involves initializing to a stabilizer state,
quantum operations with positive Wigner functions, and measurements defined by
positive operator valued measures (POVM) $\\{M_{k}\\}$ with positive Wigner
functions can be efficiently simulated on a classical computer Mari_2012 ;
Veitch_2012 . The essential idea is that in a circuit where all elements have
positive Wigner functions, the Wigner function of the initial state $\psi_{0}$
constitutes a genuine probability distribution on phase space, the Wigner
function of each quantum operation constitutes a stochastic matrix, and the
Wigner function of each POVM $M_{k}$ constitutes a probability over $k$ for
every point in phase space. The outcome of the quantum circuit:

$P_{k}=\mathrm{Tr}\left(M_{k}\,\mathcal{E}_{m}\cdots\mathcal{E}_{1}(|\psi_{0}\rangle\langle\psi_{0}|)\right),$
(40)

translates to classical stochastic evolution in terms of the Wigner function
representations:

$P_{k}=\sum_{\vec{\alpha}}\sum_{\vec{\beta_{1}}\cdots,\vec{\beta}_{m}}\sum_{\vec{\gamma}}W_{M_{k}}(\vec{\alpha})\,W_{\mathcal{E}_{m}}(\vec{\beta_{m}}|\vec{\beta}_{m-1})\cdots
W_{\mathcal{E}_{1}}(\vec{\beta_{1}}|\vec{\gamma})\,W_{\psi_{0}}(\vec{\gamma}).$
(41)

Thus, quantum circuits which use elements with positive Wigner functions
essentially reduce to probabilistic classical dynamics. As long as each of
these probability distributions can be efficiently sampled,121212This is
usually enforced by requiring local quantum operations, i.e., those which act
on a few qubits at a time. then we can also efficiently sample the output
probability distribution of the quantum circuit on a classical computer.
Interestingly, one can generalize this argument to the case where circuit
elements are not Wigner positive, and one finds that the more the Wigner
negativity in the circuit, the less efficient the classical simulation
Pashayan_2015 ; Wang_2019 .

###  2.4 Resource theory of quantum computation

Based on the above discussion, it is natural to expect that the negativity of
the Wigner function is a measure of the non-classicality – or the inherent
quantumness – of a given state, or of a quantum circuit, and indeed, there is
an operational sense in which this is the case. It is clear that circuits
which do not involve Wigner negativity only evolve within a very small subset
of states in the Hilbert space. Thus, such circuit models cannot be universal.
In order to achieve universal quantum computation, one needs an additional
resource, namely non-stabilizer states Bravyi_2005 . A state is said to be
_magic_ if it is not a stabilizer state. Much like entanglement plays the role
of a resource in the theory of quantum communication with local unitary
operations and classical communication being treated as “cheap”, _magic_ plays
the role of a resource in the theory of quantum computation, with _stabilizer
protocols_ , i.e., quantum operations using elements with positive Wigner
functions (such as stabilizer states, Clifford unitaries and measurements in
the computational basis) being treated as being cheap Veitch_2014 .

In order to quantify the notion of magic, one defines a _magic monotone_
$\mathcal{M}$ as any real-valued measure on density matrices which is, on
average, non-increasing under any stabilizer protocol Veitch_2014 . By on
average, we mean that if $\Lambda$ is a stabilizer protocol potentially
involving a measurement, and

$\Lambda(\rho)=\sigma_{i},$ (42)

with probability $p_{i}$, then

$\mathcal{M}(\rho)\geq\sum_{i}p_{i}\mathcal{M}(\sigma_{i}).$ (43)

For instance, it could happen that
$\mathcal{M}(\rho)<\mathcal{M}(\sigma_{i})$; we only require the monotone to
decrease on average.131313Note also that the monotone is not necessarily
required to be convex. It turns out that the _logarithmic negativity_ of the
Wigner function, sometimes also called _mana_ :

$\mathcal{M}(\psi)=\log\,\mathcal{N}(\psi),\;\;\;\mathcal{N}(\psi)=\sum_{x,p}|W_{\psi}(x,p)|,$
(44)

is a magic monotone. For e.g., the Wigner negativity remains invariant under a
Clifford unitary, and more generally either decreases under either a quantum
channel with positive Wigner function, or decreases on average under a POVM
with a positive Wigner function Wang_2019 . The logarithmic negativity has
several nice properties. It is evident from the definition that the log-
negativity is positive definite, and is zero if and only if the Wigner
function is nowhere negative. The log-negativity is also bounded above by
$\frac{1}{2}\log\,D$. To see this, we note from equation (20) that

$\sum_{\vec{\alpha}}W_{\vec{\alpha}}^{2}=\frac{1}{D}\Tr(\rho^{2}).$ (45)

Together with Jensen’s inequality and the fact that $\mathrm{Tr}\,\rho^{2}\leq
1$, this gives us the desired bound:

$\sum_{\vec{\alpha}}|W_{\vec{\alpha}}|\leq
D\sqrt{\sum_{\vec{\alpha}}W_{\vec{\alpha}}^{2}}=\sqrt{D}\sqrt{\mathrm{Tr}\,\rho^{2}}\leq\sqrt{D}.$
(46)

Most importantly, the log-negativity has a very direct operational meaning, in
terms of _magic distillation_. Any stabilizer protocol which consumes resource
states $\rho$ to produce $m$ copies of a target state $\sigma$ with some
finite probability, requires, on average, at least
$m\frac{\mathcal{M}(\sigma)}{\mathcal{M}(\rho)}$ copies of $\rho$ Veitch_2014
. This follows from the additivity of logarithmic negativity under tensor
factorization, together with the fact that the logarithmic negativity
decreases on average under stabilizer operations. Thus, the negativity of the
Wigner function constitutes an operationally meaningful notion of the inherent
“quantumness” of a state. In the rest of this work, we will mainly focus on
the negativity as opposed to the log-negativity:

$\mathcal{N}(\psi)=\sum_{x,p}|W_{\psi}(x,p)|,$ (47)

and regard it as a measure of the non-classicality of a state.

##  3 The Krylov-Wigner function

Having discussed the importance of Wigner negativity in the theory of quantum
computation, our aim now is to study the behavior of the Wigner function and
its negativity under time evolution with a chaotic Hamiltonian for generic
initial states. But before we go on to study the time evolution of the
negativity, we have to make a choice for the basis with respect to which we
define the Wigner function.

###  3.1 Krylov basis

Most physical Hamiltonians come with a natural locality structure, i.e., a
tensor product structure for the Hilbert space. For instance, the physical
system could be made out of qubits. The Hamiltonian is then written as a
linear combination of sufficiently local operators, i.e., products of finitely
many single qubit operators. The locality structure, of course, presents a
natural candidate for a computational basis. But here, we focus on another
potentially natural candidate called the _Krylov basis_ , which is defined as
follows: given a normalized initial state $\psi_{0}$, we generate a basis by
successively applying the Hamiltonian to $\psi_{0}$ and then iteratively
orthogonalizing using the Graham-Schmidt process. So, we denote $\psi_{0}$ to
be the state $|0\rangle$. Then, we consider the state $H|0\rangle$,
orthogonalize it with respect to $|0\rangle$ and normalize to get the state
$|1\rangle$:

$|1\rangle=\frac{1}{\sqrt{\langle 0|H^{2}|0\rangle-(\langle
0|H|0\rangle)^{2}}}\left(H|0\rangle-\langle 0|H|0\rangle|0\rangle\right),$
(48)

Then we take the state $H^{2}|0\rangle$, orthogonalize with respect to
$|0\rangle$ and $|1\rangle$, normalize and call it $|2\rangle$, and keep
going. For a chaotic Hamiltonian $H$ without any symmetries, and a
generic141414For instance, the initial state should not be an energy
eigenstate. Any sufficiently random superposition of energy eigenstates would
suffice. initial state $\psi_{0}$, we expect that this procedure generates a
basis that spans the entire Hilbert space.

The action of the Hamiltonian in the Krylov basis takes the following simple
form:

$H\ket{k}=a_{k}\ket{k}+b_{k+1}\ket{k+1}+b_{k}\ket{k-1},$ (49)

where

$a_{n}=\langle n|H|n\rangle,$ (50)

and the $b_{n}$ are real, positive coefficients (with the convention that
$b_{0}=b_{D}=0$). In matrix form, the Hamiltonian looks like:

$H=\begin{pmatrix}a_{0}&b_{1}&0&0&\cdots&0&0\\\
b_{1}&a_{1}&b_{2}&0&\cdots&0&0\\\ 0&b_{2}&a_{2}&b_{3}&\cdots&0&0\\\
0&0&b_{3}&a_{3}&\cdots&0&0\\\
\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\\
0&0&0&0&\cdots&a_{D-2}&b_{D-1}\\\ 0&0&0&0&\cdots&b_{D-1}&a_{D-1}\\\
\end{pmatrix}$ (51)

The Krylov basis has the nice property that it minimizes the “spreading” of
the wavefunction under time evolution Balasubramanian:2022tpr . This can be
quantified as follows: consider the time evolved state
$|\psi(t)\rangle=e^{-itH}|\psi_{0}\rangle$. For any ordered, orthonormal basis
$\mathcal{B}=\left\\{|n_{\mathcal{B}}\rangle\right\\}_{n=0}^{D-1}$, we can
defined a measure of the spread of the state $\psi(t)$:

$C_{\mathcal{B}}(t)=\sum_{n=0}^{D-1}c_{n}|\langle
n_{\mathcal{B}}|\psi(t)\rangle|^{2},$ (52)

where $c_{n}$ is any positive, monotonically increasing weight function; for
e.g., a natural choice is to take $c_{n}=n$. We may think of
$C_{\mathcal{B}}(t)$ as some notion of “spread complexity”, i.e., a measure of
how spread out the state is in the given choice of basis. Let us say we choose
$|0\rangle_{\mathcal{B}}=|\psi_{0}\rangle$, so that $C_{\mathcal{B}}(0)$ is
set to zero. This does not fix anything about the remaining $(D-1)$ basis
vectors of $\mathcal{B}$. But one could ask for more: one could look for a
basis where the early growth of $C_{\mathcal{B}}(t)$ with respect to time is
minimized. In other words, we would like a basis which successively minimizes
the time derivatives
$\\{C^{\prime}_{\mathcal{B}}(0),C^{\prime\prime}_{\mathcal{B}}(0),\cdots\\}$
and so on up to the $(D-1)$th derivative. It was shown in
Balasubramanian:2022tpr  that the Krylov basis minimizes all of these
derivatives. In this sense, the Krylov basis quantifies the minimal and
inevitable spread in the wavefunction that a state undergoes under quantum
evolution. To be clear, adding individual phases to the basis vectors does not
change $C_{\mathcal{B}}(t)$, so any generalized Krylov basis with individual
phases also minimizes the spreading.

Incidentally, the Krylov basis naturally shows up in an interesting way in the
context of the double-scaled SYK model Berkooz:2018jqr ; Lin:2022rbf ;
Rabinovici:2023yex . The DSSYK model is a model of $N$ Majorana fermions with
all-to-all $p$-body interactions. The couplings in the Hamiltonian are chosen
randomly from a Gaussian ensemble, and one takes a double-scaling limit where
$p\to\infty$, $N\to\infty$ with $\lambda=\frac{2p^{2}}{N}$ fixed. The model is
of particular interest because JT gravity arises in a certain triple scaling
limit of the DSSYK model. It turns out that calculations in the DSSYK model
can be organized nicely in a particular basis called the “chord-number” basis.
This basis is the natural generalization of the gravitational length basis
(i.e., a basis for the gravitational Hilbert space labelled by the length of
the bulk wormhole) in JT gravity. From the perspective of the boundary quantum
system, the chord-number basis turns out to be the Krylov basis, and the
length of the bulk wormhole maps onto the Krylov index which orders the basis.
Thus, the dual “gravitational” description of the DSSYK model naturally seems
to utilize the Krylov basis. It is natural to wonder what makes the Krylov
basis special, i.e., whether there is a deeper underlying reason why the
gravitational variables become manifest in the Krylov basis. We will take some
baby steps towards addressing this question in section 3.3.

###  3.2 Krylov phase space

The results of Balasubramanian:2022tpr  suggest that the Krylov basis is very
natural in the sense that the time-dependent wavefunction is as peaked as
possible in this basis. This suggests that the state looks as “classical” as
possible in the Krylov basis, somewhat like coherent states, except that the
latter are peaked in phase space. In order to show that there is a sense in
which the Krylov basis makes the time-evolving state as classical as possible,
it thus seems natural to study the time-evolution in phase space. Following
our discussion in section 2.2, we can define a discrete phase space
$\mathcal{P}=\mathbb{Z}_{D}\times\mathbb{Z}_{D}$ with respect to the Krylov
basis, and a corresponding Wigner function for any state $\psi$:

$W_{\psi}(q,p)=\frac{1}{D}\sum_{k,\ell=0}^{D-1}\widehat{\delta}_{2q,k+\ell}e^{\frac{2\pi
i}{D}(k-\ell)p}\langle k|\psi\rangle\langle\psi|\ell\rangle.$ (53)

Of particular interest to us will be the Wigner function of the time evolved
initial state $|\psi(t)\rangle=e^{-itH}|\psi_{0}\rangle$. We will call this
object the _Krylov-Wigner function_. The Krylov-Wigner function gives us a way
of visualizing how the state $|\psi\rangle\langle\psi|$ evolves over the
discrete Krylov phase space. In words, one can expand the evolved density
matrix $|\psi(t)\rangle\langle\psi(t)|$ in terms of the discrete phase-point
operators – which recall form an orthonormal basis for the space of operators
– and the KW-function is the corresponding wavefunction.151515Note that there
is a different thing one could have done, which is to construct an operator
version of the Krylov basis by constructing the operators
$\left\\{\rho,[H,\rho],[H,[H,\rho]],\cdots\right\\}$ where
$\rho=|\psi\rangle\langle\psi|$, and then orthonormalize these with respect to
the trace norm. However, this would not give a natural phase-space structure.
It’s easy to see that the KW function at time $t=0$ is given by

$W_{\psi_{0}}(q,p)=\frac{1}{D}\delta_{q,0}.$ (54)

In order to see how the Wigner function evolves with time, we can do a Taylor
expansion in $t$:

$\displaystyle W_{\psi(t)}(q,p)$ $\displaystyle\simeq$
$\displaystyle\frac{1}{D}\delta_{q,0}-\frac{it}{D}\sum_{k,\ell}\widehat{\delta}_{2q,k+\ell}e^{\frac{2\pi
i}{D}(k-\ell)p}\left(\langle k|H|0\rangle\langle 0|\ell\rangle-\langle
k|0\rangle\langle 0|H|\ell\rangle\right))+\cdots$ (55) $\displaystyle=$
$\displaystyle\frac{1}{D}\delta_{q,0}-\frac{2t}{D}\widehat{\delta}_{2q,1}\sin\left(\frac{2\pi
p}{D}\right)+\cdots.$

Note that the second term is only non-zero when $q=\frac{D+1}{2}$. So, it
would naively seem like the time evolution spreads the Wigner function non-
locally. However, this is only a matter of convention and can be easily cured
by representing the Wigner function as being defined on the lattice of half-
integer separated points
$\left\\{0,\frac{1}{2},1,\cdots,\frac{D-1}{2}\right\\}$, with the
understanding that the point $Q=\frac{1}{2}$ in this lattice corresponds to
$q=\frac{D+1}{2}$ in $\mathbb{Z}_{D}$.161616We will use $Q$ to denote half-
integer separated points. This is merely a re-ordering of the phase-point
operators, and a convenient one because with respect to this ordering, the
spread of the Wigner function has a more local appearance. For instance, at
first order in the expansion, the KW function spreads to the point
$Q=\frac{1}{2}$, at second order it spreads up to $Q=1$, and so on.

We can now express the Hamiltonian as a function on the Krylov phase space:

$\displaystyle H(Q,p)$ $\displaystyle=$ $\displaystyle\frac{1}{D}\Tr(HA(Q,p))$
(56) $\displaystyle=$
$\displaystyle\frac{1}{D}\sum_{k,\ell=0}^{D-1}\widehat{\delta}_{2Q,k+\ell}e^{\frac{2\pi
i}{D}(k-\ell)p}\langle\ell|H|k\rangle$ $\displaystyle=$
$\displaystyle\frac{1}{D}\sum_{k,\ell=0}^{D-1}\widehat{\delta}_{2Q,k+\ell}e^{\frac{2\pi
i}{D}(k-l)p}\Tr\bigg(a_{k}\delta_{k,\ell}+b_{k+1}\delta_{k+1,\ell}+b_{k}\delta_{k-1,\ell}\bigg{missing})$
$\displaystyle=$
$\displaystyle\frac{1}{D}\left(a_{Q}+b_{Q+\frac{1}{2}}\cos\frac{2\pi
p}{D}\right),$

where in the third line we have used equation (49) for how the Hamiltonian
acts on the Krylov basis vectors. In phase space, the time evolution of the KW
function is controlled by equation (24), with the Hamiltonian function given
above.

###  3.3 Minimizing the growth of negativity

In this section, we would like to show that there is a sense in which time
evolution in the Krylov basis looks as classical as possible. For this
purpose, we consider the generalized family of Krylov bases obtained by
multiplying individual basis elements by phases:

$\mathcal{K}=\left\\{e^{i\phi_{k}}|k\rangle\right\\}_{k=0}^{D-1},$ (57)

where $|k\rangle$ is the $k$th basis vector of _the_ Krylov basis, as defined
around equation (48) (with no extra phases). The claim is that for some
particular choice of the phases $\\{\phi_{k}\\}$, the generalized Krylov basis
minimizes the early time growth of Wigner negativity in the large-$D$ limit.
This result is quite analogous to the one of Balasubramanian:2022tpr , where
it was shown that the Krylov basis minimizes the growth of the spread
complexity, and our derivation below will follow a similar strategy as in
Balasubramanian:2022tpr .

It is worth taking a moment to explain what we mean by “minimizing the early
time growth of the Wigner negativity”. Let $\mathcal{B}$ be some basis, and
let $W_{\mathcal{B}}$ be the Wigner function defined with respect to that
basis. If we wish to minimize the Wigner negativity at $t=0$, then it is clear
that we must pick one of the basis vectors to be $\psi_{0}$, and without loss
of generality we call this vector $|0\rangle_{\mathcal{B}}=|\psi_{0}\rangle$.
With this choice, the log of the Wigner negativity at $t=0$ vanishes, and is
thus also minimized. Now, consider the functions:

$\overline{W}^{(m)}_{\mathcal{B}}(Q,p)=\sum_{n=0}^{m}W^{(n)}_{\mathcal{B}}(Q,p)\frac{t^{n}}{n!},$
(58)

where
$W^{(n)}_{\mathcal{B}}=\frac{d^{n}W_{\mathcal{B}}}{dt^{n}}\Big{|}_{t=0}$.
These functions $\overline{W}^{(m)}_{\mathcal{B}}$ are essentially
increasingly accurate Taylor series approximations to the Wigner function in
an infinitesimally small neighborhood of $t=0$. Let

$\mathcal{S}_{\mathcal{B}}=\left\\{\mathcal{N}_{\mathcal{B}}^{(0)},\mathcal{N}_{\mathcal{B}}^{(1)}\cdots,\mathcal{N}_{\mathcal{B}}^{(D-1)}\right\\},$
(59)

where

$\mathcal{N}^{(m)}_{\mathcal{B}}=\sum_{Q,p}|\overline{W}^{(m)}_{\mathcal{B}}(Q,p)|.$
(60)

is the negativity of the $m$th Taylor series approximation.

Definition: For two bases $\mathcal{B}_{1}$ and $\mathcal{B}_{2}$, we say that
$\mathcal{S}_{\mathcal{B}_{1}}<\mathcal{S}_{\mathcal{B}_{2}}$ if there exists
some $k<D$ such that
$\mathcal{N}^{(m)}_{\mathcal{B}_{1}}=\mathcal{N}^{(m)}_{\mathcal{B}_{2}}$ for
all $m<k$, and
$\mathcal{N}^{(k)}_{\mathcal{B}_{1}}<\mathcal{N}^{(k)}_{\mathcal{B}_{2}}$. We
say that a basis $\mathcal{B}$ minimizes the growth of the Wigner negativity
if for any other basis $\mathcal{B}^{\prime}$, we have
$\mathcal{S}_{\mathcal{B}}\leq\mathcal{S}_{\mathcal{B}^{\prime}}$.

Intuitively, if we imagine the time evolution as happening in discrete steps
with time intervals $\delta t$, then what we wish to do is to find a class of
bases which minimizes the negativity at the first time step, then within that
class find a sub-class which minimizes the negativity at the second time-step,
and so on. In order to minimize the negativity after the first time step, one
would naively want to include the vector $e^{-i\delta t\,H}|\psi_{0}\rangle$
as a basis element. However, one must first orthogonalize it with respect to
$\psi_{0}$. Similarly, one would naively want to keep including the vectors
$e^{-i\,n\,\delta t\,H}|\psi_{0}\rangle$ as basis vectors at successive time
steps, but after first orthogonalizing them with respect to the previous basis
vectors. In the $\delta t\to 0$ limit, this is precisely the Krylov basis. The
following result makes this intuition precise, but instead of doing it for
discrete time steps, we will try and implement this minimization at the level
of the above Taylor series approximations in time. Having explained what we
mean by minimizing the early time growth of the Wigner negativity, we can now
state our result:

Claim: _In the large- $D$ limit, the basis $\mathcal{B}$ which minimizes the
early time growth of Wigner negativity belongs to the class $\mathcal{K}$ of
generalized Krylov bases._

Proof: Our proof will be inductive in nature. We will assume that for the
basis to minimize the first $(d-1)$ elements of $\mathcal{S}_{\mathcal{B}}$,
the basis vectors $\\{|n_{\mathcal{B}}\rangle\\}_{n=0}^{d-1}$ must be of the
generalized Krylov form, with some particular choice of phases. Now our aim is
to show that the requirement that the function
$\mathcal{N}^{(d)}_{\mathcal{B}}$ be minimized forces us to pick the $d$th
basis vector to be the $d$th Krylov basis vector, up to a phase. To this end,
consider a basis $\mathcal{B}$ such that its first $d$ basis vectors agree
with that of the Krylov basis $\mathcal{K}$, up to phases which are determined
by the requirement that the first $(d-1)$ elements of
$\mathcal{S}_{\mathcal{B}}$ are minimized. The time-derivatives of the Wigner
function at $t=0$, with respect to the basis $\mathcal{B}$, are given by:

$W_{\mathcal{B}}^{(m)}(Q,p)=\frac{1}{D}\sum_{k,\ell=0}^{D-1}\widehat{\delta}_{2Q,k+\ell}e^{\frac{2\pi
i}{D}(k-\ell)p}\sum_{n=0}^{m}(-1)^{n}i^{m}\frac{m!}{n!(m-n)!}\langle\ell_{\mathcal{B}}|H^{n}|\psi(0)\rangle\langle\psi(0)|H^{m-n}|k_{\mathcal{B}}\rangle.$
(61)

For $m\leq(d-1)$, the derivatives of the Wigner function for $\mathcal{B}$
agree precisely with that of $\mathcal{K}$. This means that the Taylor
expansion of the Wigner function around $t=0$ agrees up to and including
$O(t^{d-1})$ terms. The first deviation happens at $m=d$:

$W_{\mathcal{B}}^{(d)}(Q,p)=\frac{1}{D}\sum_{k,\ell=0}^{D-1}\widehat{\delta}_{2Q,k+\ell}e^{\frac{2\pi
i}{D}(k-\ell)p}\sum_{n=0}^{d}(-1)^{n}i^{d}\frac{d!}{n!(d-n)!}\langle\ell_{\mathcal{B}}|H^{n}|\psi(0)\rangle\langle\psi(0)|H^{d-n}|k_{\mathcal{B}}\rangle.$
(62)

In particular, the terms which are different from the generalized Krylov basis
correspond to $n=0$ and $n=d$. Now,

$H^{d}|\psi(0)\rangle=c_{d}|d_{\mathcal{K}}\rangle+|\chi\rangle,$ (63)

where $|\chi\rangle$ is some linear combination of the Krylov basis vectors of
order lower than $d$, and

$c_{d}=e^{-i\phi_{d}}\prod_{n=1}^{d}b_{n},$

where the phase $\phi_{d}$ is the phase of the $d$th basis vector in
$\mathcal{K}$ which is arbitrary for now, but will be fixed momentarily. Thus,
we get

$\displaystyle W_{\mathcal{B}}^{(d)}(Q,p)$ $\displaystyle=$
$\displaystyle\frac{1}{D}\sum_{k,\ell=0}^{D-1}\widehat{\delta}_{2Q,k+\ell}e^{\frac{2\pi
i}{D}(k-\ell)p}\sum_{n=1}^{d-1}(-1)^{n}i^{d}\frac{d!}{n!(d-n)!}\langle\ell_{\mathcal{K}}|H^{n}|\psi(0)\rangle\langle\psi(0)|H^{d-n}|k_{\mathcal{K}}\rangle$
(64) $\displaystyle+$
$\displaystyle\frac{1}{D}\sum_{\ell=0}^{D-1}\widehat{\delta}_{2Q,\ell}\left(e^{-\frac{2\pi
i}{D}\ell p}(-i)^{d}\langle\ell_{\mathcal{K}}|\chi\rangle+e^{\frac{2\pi
i}{D}\ell p}i^{d}\langle\chi|\ell_{\mathcal{K}}\rangle\right)$
$\displaystyle+$
$\displaystyle\frac{i^{d}}{D}\sum_{k=0}^{D-1}\widehat{\delta}_{2Q,k}\left(e^{\frac{2\pi
i}{D}kp}c_{d}^{*}\langle d_{\mathcal{K}}|k_{\mathcal{B}}\rangle+e^{-\frac{2\pi
i}{D}kp}(-1)^{d}c_{d}\langle k_{\mathcal{B}}|d_{\mathcal{K}}\rangle\right).$

We can re-write this as

$W_{\mathcal{B}}^{(d)}(Q,p)=W_{\mathcal{K}}^{(d)}(Q,p)+\frac{i^{d}}{D}\sum_{k=d}^{D-1}\widehat{\delta}_{2Q,k}\left(c_{d}^{*}e^{\frac{2\pi
i}{D}kp}(\langle
d_{\mathcal{K}}|k_{\mathcal{B}}\rangle-\delta_{k,d})+c_{d}e^{-\frac{2\pi
i}{D}kp}(-1)^{d}(\langle
k_{\mathcal{B}}|d_{\mathcal{K}}\rangle-\delta_{k,d})\right).$ (65)

Note that the sum above is restricted to $k\geq d$ because for $k<d$, the
basis vectors in $\mathcal{B}$ are proportional to the Krylov basis vectors,
and thus orthogonal to $|d_{\mathcal{K}}\rangle$. Let us temporarily fix
$\phi_{d}$ by requiring that $\langle d_{\mathcal{B}}|d_{\mathcal{K}}\rangle$
be a real number in the range $[0,1]$. This just picks out a particular member
in the class of generalized Krylov bases. Our aim is to show that this basis
already has lower negativity than $\mathcal{B}$, and then minimizing over
$\phi_{d}$ will only lower the negativity. Let us denote:

$\langle
k_{\mathcal{B}}|d_{\mathcal{K}}\rangle=\delta_{d,k}+\alpha_{k},\;\;\cdots\;\;(k\geq
d)$ (66)

where $\alpha_{d}\in\mathbb{R}$ and $-1\leq\alpha_{d}\leq 0$ (be choice of
$\phi_{d}$), and

$(1+\alpha_{d})^{2}+\sum_{k=d+1}^{D-1}|\alpha_{k}|^{2}=1.$

We can thus write equation (65) as follows:

$W_{\mathcal{B}}^{(d)}(Q,p)=\begin{cases}&\displaystyle
W_{\mathcal{K}}^{(d)}(Q,p)+\frac{2(-1)^{\frac{d}{2}}}{D}|c_{d}||\alpha_{2Q}|\Theta_{2Q\geq
d}\cos\left(\frac{2\pi(2Q)p}{D}+\theta_{2Q}+\phi_{d}\right)\;\;\cdots(d\;\text{even})\\\
&\displaystyle
W_{\mathcal{K}}^{(d)}(Q,p)+\frac{2(-1)^{\frac{d-1}{2}}}{D}|c_{d}||\alpha_{2Q}|\Theta_{2Q\geq
d}\sin\left(\frac{2\pi(2Q)p}{D}+\theta_{2Q}+\phi_{d}\right)\;\;\cdots(d\;\text{odd}),\end{cases}$
(67)

where $\alpha_{k}=|\alpha_{k}|e^{i\theta_{k}}$, $c_{d}=|c_{d}|e^{-i\phi_{d}}$
and $\Theta_{2Q\geq d}$ is a step function which is one only when $2Q\geq d$
and zero otherwise. The negativity at this order is given by

$\mathcal{N}^{(d)}=\sum_{Q,p}|\overline{W}^{(d)}_{\mathcal{B}}(Q,p)|.$ (68)

Including terms up to $O(t^{d-1})$, the Wigner function has non-trivial
support till $2Q=d-1$. After including the $O(t^{d})$ terms, the Wigner
function in this region is exactly equal to that of the Wigner function in the
generalized Krylov basis, i.e., the difference in basis vectors beyond the
$(d-1)$the vector does not manifest in this region:

$\overline{W}^{(d)}_{\mathcal{B}}(Q,p)=\overline{W}^{(d)}_{\mathcal{K}}(Q,p),\;\;\;\forall\;\;2Q<d.$
(69)

So we need only look at the negativity on the region $2Q\geq d$. In this
region, there is no contribution from lower order derivatives, so we wish to
compare $\sum_{p,2Q\geq d}|W^{(d)}_{\mathcal{B}}(Q,p)|$ with $\sum_{p,2Q\geq
d}|W^{(d)}_{\mathcal{K}}(Q,p)|$.

Let us consider $d$ even; the $d$ odd case works analogously. From the
triangle inequality applied to equation (67), we get

$\left|W_{\mathcal{K}}^{(d)}(d/2,p)\right|\leq\left|W_{\mathcal{B}}^{(d)}(d/2,p)\right|+\frac{2}{D}|c_{d}||\alpha_{d}|\left|\cos\left(\frac{2\pi
dp}{D}+\theta_{d}+\phi_{d}\right)\right|.$ (70)

Thus, we find

$\displaystyle\sum_{Q\geq d/2}|W_{\mathcal{K}}^{(d)}(Q,p)|$
$\displaystyle\leq$ $\displaystyle\sum_{Q\geq
d/2}|W_{\mathcal{B}}^{(d)}(x,p)|+\frac{2}{D}|c_{d}|\Big{\\{}|\alpha_{d}||\cos\left(\frac{2\pi
dp}{D}+\theta_{d}+\phi_{d}\right)|$ (71) $\displaystyle-$
$\displaystyle\sum_{Q>d/2}|\alpha_{2Q}||\cos\left(\frac{4\pi
Qp}{D}+\theta_{2Q}+\phi_{d}\right)|\Big{\\}}.$

Now we need to sum the above inequality over the momentum direction. In the
large $D$ limit, the sums over the terms in the curly brackets can be
approximated by integrals, and this gives

$\sum_{p,Q\geq d/2}|W_{\mathcal{K}}^{(d)}(Q,p)|\leq\sum_{p,Q\geq
d/2}|W_{\mathcal{B}}^{(d)}(Q,p)|+\frac{4|c_{d}|}{\pi}\left(d|\alpha_{d}|-\sum_{k>d}k|\alpha_{k}|\right).$
(72)

Our goal now is to show that $d|\alpha_{d}|\leq\sum_{k>d}k|\alpha_{k}|$. We
know that

$(1+\alpha_{d})^{2}+\sum_{k>d}|\alpha_{k}|^{2}=1,$ (73)

so using the fact that $\alpha_{d}$ is real and negative, we get

$-2|\alpha_{d}|+|\alpha_{d}|^{2}+\sum_{k>d}|\alpha_{k}|^{2}=0.$ (74)

Since $|\alpha_{d}|$ and $|\alpha_{k}|$ are all positive numbers less than 1,
we get

$-|\alpha_{d}|+\sum_{k>d}|\alpha_{k}|\geq 0,$ (75)

which then implies $|\alpha_{d}|\leq\sum_{k>d}|\alpha_{k}|$. From here, it
immediately follows that

$d|\alpha_{d}|\leq\sum_{k>d}k|\alpha_{k}|.$ (76)

Thus, we have shown that for any given basis $\mathcal{B}$, we can choose the
phase of the $d$th generalized Krylov basis vector in a way such that we can
always make the negativity of the generalized Krylov basis smaller than the
given basis $\mathcal{B}$. Thus, if we pick the phase of the $d$th Krylov
basis vector by simply minimizing the negativity over the choice of this
phase, we are guaranteed that there exists a choice of a generalized Krylov
basis which minimizes the negativity up to this order.$\Box$

In summary, the generalized Krylov basis with appropriately chosen phases
minimizes the early growth of the Wigner negativity in the large-$D$ limit.
So, if one’s goal is to try as hard as possible to simulate quantum dynamics
on a classical computer, then it is natural to expect that the Krylov basis
(with appropriately chosen phases) would be a good bet for the computational
basis. An important caveat to bear in mind is that we have not given an
algorithm here to demonstrate efficient classical simulation of time evolution
in the Krylov basis, although we take our results as some evidence in this
direction. Still, it is natural to wonder whether this is what makes the
Krylov basis special, and in particular, a natural candidate for a “dual
gravitational basis” for strongly coupled quantum systems, with equation (56)
being the natural analog of the gravitational Hamiltonian. Remarkably, in the
DSSYK model, this is exactly what happens – the dual “gravitational”
description consists of the chord-length basis, which precisely coincides with
the Krylov basis Lin:2022rbf ; Rabinovici:2023yex . So, it seems plausible
that _the gravitational variables correspond to a choice of basis which makes
classical simulation of quantum dynamics as efficient as possible_. There are
many important details which need to be worked out to make this proposal more
concrete, some of which are listed in the discussion section.

##  4 Krylov-Wigner negativity in random matrix theory

In this section, we will study the growth of the negativity of the KW function
in random matrix theory. The DSSYK model or the SSS matrix model would perhaps
be the most natural models to consider. But as a first step, here we will
consider Gaussian random matrix theory. In the following sections, we will
consider the Krylov basis without any added phases. It would be interesting to
explore addition of phases to these results.

###  4.1 Setup

We start with some working basis, which we will call the “coordinate basis”.
We pick the initial state $\psi_{0}$ to coincide with one of the states in
this basis, say, $|0\rangle$. Next, we pick a $D\times D$ matrix from the
Gaussian Unitary ensemble (GUE) and declare that as the Hamiltonian in the
coordinate basis. The GUE corresponds to the following probability
distribution over the space of Hermitian matrices:

$p(H)=\frac{1}{Z_{\mathrm{GUE}}}e^{-\frac{D}{2}\mathrm{Tr}(H^{2})}$

where $Z_{\mathrm{GUE}}=2^{\frac{D}{2}}(\frac{\pi}{D})^{\frac{D^{2}}{2}}$ is
the normalization. Since the Hamiltonian is random, we expect that the
$|0\rangle$ state is a sufficiently generic superposition of energy
eigenstates of $H$. In the context of holography, we can think of this choice
of the initial state as a typical black hole microstate at infinite
temperature, along the lines of Kourkoulou:2017zaj . We now build the Krylov
basis with $|0\rangle$ as the first state, and then successively acting with
the Hamiltonian and orthonormalizing using the Gram-Schmidt process. Having
built a Krylov basis for the Hilbert space, we then construct the KW function
and compute its negativity as a function of time:

$\mathcal{N}(t)=\sum_{Q,p}|W_{\psi(t)}(Q,p)|,$ (77)

where $|\psi(t)\rangle=e^{-itH}|\psi_{0}\rangle$.

|
---|---
Figure 1: Negativity of the Wigner function $\mathcal{N}$ as a function of the
re-scaled time $\frac{t}{\sqrt{D}}$ for a randomly chosen Hamiltonian from the
GUE for various values of $D$.

###  4.2 Spread, ramp and plateau

In figure 1, we plot the negativity of the Wigner function as a function of
the re-scaled time $\frac{t}{\sqrt{D}}$. The numerical plots broadly show at
least three distinct phases for the negativity: (i) a gradual early-time
growth which we will call the _spreading_ phase, (ii) a sharp _ramp_ at a time
$t\sim O(\sqrt{D})$ where the negativity grows approximately linearly with a
slope (with respect to $t$) of $O(\sqrt{D})$, and (iii) saturation to a
_plateau_ at approximately $\mathcal{N}_{\text{sat.}}\sim 0.8\sqrt{D}$ at late
times ($t\gg\sqrt{D}$). Note that the saturation value is close to but less
than the upper bound of $\sqrt{D}$.

|  |
---|---|---
|  |
|  |
Figure 2: Density plots of the Krylov-Wigner function in phase space for
$D=101$ for various values of time corresponding to
$\frac{t}{\sqrt{D}}=(0.125,0.5,1,1.5,2,2.25,2.5,5,1000)$.

At $t=0$, the Wigner function is localized at $Q=0$, but spread uniformly over
the momentum direction:

$W(0,Q,p)=\frac{1}{D}\delta_{Q,0}.$ (78)

As time increases, the Wigner function spreads in phase space (see figure 2).
In the spreading phase, the Wigner negativity grows gradually with time. This
growth seems to be much slower than would be the case had we picked any other
computational basis. For instance, in figure 3 we compare the growth of the
negativity between the Krylov basis and the original coordinate basis, i.e.
the basis with respect to which the entries of the Hamiltonian are picked to
be random. Note that the negativity in the coordinate basis rises sharply
starting at $t=0$, with a slope of $O(\sqrt{D})$, and saturates within a time
of $O(1)$. We expect this to be the behaviour in any generic basis.

In contrast, the KW negativity shows a much more gradual rise, with an $O(1)$
slope at $t=0$. To compute the slop, recall that the Wigner function can be
expanded perturbatively in time as follows:

$\displaystyle W(t,Q,p)$ $\displaystyle=$
$\displaystyle\frac{1}{D}\sum_{k,\ell}\widehat{\delta}_{2Q,k+\ell}e^{\frac{2\pi
i}{D}(k-\ell)p}\left(\bra{\ell}\rho(0)\ket{k}-it\bra{\ell}[H,\rho(0)]\ket{k}+\mathcal{O}(t^{2})\right)$
$\displaystyle=$
$\displaystyle\frac{1}{D}\sum_{k,\ell}\widehat{\delta}_{2Q,k+\ell}e^{\frac{2\pi
i}{D}(k-\ell)p}\Big{(}\langle\ell|0\rangle\langle 0|k\rangle-
it(\langle\ell|H|0\rangle\langle 0|k\rangle-\langle\ell|0\rangle\langle
0|H|k\rangle)+\mathcal{O}(t^{2})\Big{)}.$

Using

$H|0\rangle=a_{0}|0\rangle+b_{1}|1\rangle,$ (80)

we get

$\displaystyle W(t,Q,p)$ $\displaystyle=$
$\displaystyle\frac{1}{D}\delta_{2Q,0}-\frac{2tb_{1}}{D}\sin\left(\frac{2\pi
p}{D}\right)\delta_{2Q,1}+O(t^{2}).$ (81)

At linear order in time, the Wigner function is only supported at $Q=0$ and
$Q=\frac{1}{2}$, so we can compute the negativity as follows:

$\mathcal{N}=1+\frac{2tb_{1}}{D}\sum_{p=0}^{D-1}\left|\sin(\frac{2\pi
p}{D})\right|+O(t^{2}).$ (82)

In the large $D$ limit, we can convert the sum into an integral, and we get

$\mathcal{N}=1+\frac{2t}{2\pi}b_{1}\int_{0}^{2\pi}d\theta|\sin(\theta)|=1+\frac{4t}{\pi}b_{1}+O(t^{2}).$
(83)

Thus, the negativity always starts rising linearly with a slope proportional
to $b_{1}$. We expect the expectation value of $\overline{b_{1}}$ in GUE to be
$O(1)$ Balasubramanian:2022tpr ; Balasubramanian:2022dnj , so that the slope
is approximately equal to $\frac{4}{\pi}$ times an $O(1)$ constant (see figure
3 for comparison with numerics). We can also give a non-linear approximation
for the early time behaviour of the KW function if we assume all the $b_{n}$
for small enough $n$ are close to one (see Appendix A).

|
---|---
Figure 3: (Left) The negativity of the Wigner function with respect to the
coordinate basis in comparison with that of the Krylov-Wigner function for
$D=101$. (Right) The KW negativity for small times and $D=151$. The blue curve
is the numerical result, while the dashed orange curve is a straight line with
slope $\frac{4}{\pi}$.

As $\frac{t}{\sqrt{D}}$ becomes $O(1)$, the Wigner function starts spreading
out over the entire phase space. Up until this point, the negativity is
smaller than the upper bound of $\sqrt{D}$. Once the spreading wave covers all
of phase space, there is a sharp, linear rise in the negativity (see figure 1)
which we call the _ramp_. During the ramp phase, the Wigner negativity rises
almost linearly with a slope of $O(\sqrt{D})$ and hits the saturation value
within an $O(1)$ amount of time.171717Presumably, the reason for this is that
the Krylov basis minimizes the growth of the Wigner negativity up to a point,
but we only have the freedom to fix $D$ basis vectors by minimization. Once
all this freedom is exhausted, the negativity shoots up. After the ramp, the
negativity saturates to a value of approximately $0.8\sqrt{D}$, close to its
maximum possible value of $\sqrt{D}$; we call this phase the _plateau_. In the
plateau phase, the Wigner function has spread over the entire phase space
uniformly, with a large number of positive and negative pockets uniformly
distributed throughout. Our numerical results are consistent with the claim
that the Krylov basis minimizes the early-time growth of Wigner negativity,
and further suggest that this basis is useful for times not scaling
exponentially in $\log\,D$. Note that the growth of the Wigner negativity for
an exponential (in $\log\,D$) amount of time followed by saturation to an
exponentially large value at exponential time follows the generally expected
behavior of measures of computational complexity Brown_2018 .

##  5 Discussion

In this work, we have introduced and studied the Krylov-Wigner function, i.e.,
the discrete Wigner function defined with respect to the Krylov basis. We
showed that with an appropriate choice of phases, the Krylov basis minimizes
the growth of the negativity of the Wigner function for large $D$ quantum
systems. We then studied the negativity of the Krylov-Wigner function in
Gaussian random matrix theory and observed broadly three phases in its time
evolution: spread, ramp and plateau. Much like circuit complexity Brown_2018 ;
Balasubramanian:2019wgd ; Balasubramanian:2021mxo , the KW negativity also
keeps growing till a time exponential in $\log\,D$, and then saturates at an
exponentially large value. We end with some remarks on the outlook and future
directions.

###  5.1 Gravity as efficient computation

The Krylov basis has appeared recently in an interesting way in the context of
the DSSYK model Berkooz:2018jqr , where it emerges naturally in the dual
gravitational description as the eigenstates of the length operator of the
bulk wormhole Lin:2022rbf ; Rabinovici:2023yex . We hope that our results shed
new light on this point. The negativity of the Wigner function is a measure of
the intrinsic non-classicality in a state. The fact that the Krylov basis
minimizes the early-time growth of the Wigner negativity in the large-$D$
limit suggests that this basis is ideally suited for a dual, semi-classical
description of chaotic quantum dynamics at large $D$. Here, the aspect of
chaos that is important is ergodicity, i.e., we want the Krylov basis to span
the entire (or at least a sufficiently large subspace of the) Hilbert space.
From this perspective, it is natural to think of the Krylov phase space and
the corresponding Hamiltonian as the dual “gravitational” theory, which seeks
to implement the boundary quantum dynamics on a classical computer to the best
extent possible, supplemented with some additional non-classical resource.
This motivates us to propose a new computational principle for the emergence
of spacetime and gravity – _the gravitational variables correspond to a choice
of basis which makes classical implementation of quantum computation as
efficient as possible_. Admittedly, at this stage, this a rather vague and
rudimentary proposal, and there are several issues which need to be addressed
in order to make it more concrete:

  1. 1.

Firstly, this “dual description” in the Krylov basis appears to come with a
discrete phase space at finite $D$. Presumably, continuum gravity emerges from
an appropriate coarse-graining in the large-$D$ limit. Roughly, one would like
to first send $D\to\infty$. If the Krylov coefficients $b_{n}$ (for $n$ not
scaling with $D$) are sufficiently smooth and slowly varying functions of $n$
(perhaps after ensemble averaging) with some effective Planck scale $\lambda$
which cuts off the large-$n$ growth of these coefficients, then one can take a
continuum limit where $\lambda\to 0,\,n\to\infty$ with $\ell=n\lambda$ held
fixed (see Lin:2022rbf ; Rabinovici:2023yex  for how JT gravity emerges from
the discrete theory in the case of the DSSYK model).

  2. 2.

The discrete phase space formalism seems to work best for prime $D$, and it
would be nice to have a more uniform treatment for any $D$. Presumably one can
get around this for non-prime dimensions by extending the Hilbert space to
include some auxiliary states so as to enlarge the Hilbert space dimension to
the nearest prime number. It would also be useful to generalize our results to
the continuum setting.

  3. 3.

The notion of “classicality” that is familiar from AdS/CFT is one of large-$N$
factorization, i.e, large $N$ suppresses quantum fluctuations in some simple
observables Jaffe . On the other hand, Wigner negativity seems like it has
more to do with efficiency of sampling the outcome probability distribution of
a quantum circuit on a classical computer. It would be good to clarify how
these two notions of classicality are related.

  4. 4.

The classical algorithms which seek to sample from the output probability
distribution of the quantum circuit Mari_2012 ; Veitch_2012 ; Pashayan_2015 ;
Wang_2019  are efficient when the Wigner quasiprobability distributions (or
more precisely, their absolute values) for each of the circuit elements can be
efficiently sampled. In a context where the circuit is made out of multiple
qubits with local gates acting upon them, this is guaranteed by locality. In
our context, one does not obviously have such a locality structure and so it
is not straightforward to establish classical efficiency. It is worth noting
that positive Wigner functions are highly constrained, and thus very special
probability distributions. It is plausible that there is a sense in which they
are easy to sample from, and that a small amount of negativity does not add
much complexity. We would like to emphasize that this is an important issue to
resolve in order to convincingly establish the connection between Wigner
negativity and classicality in our context.

  5. 5.

In our discussion, it was important that the Hamiltonian be chaotic, to the
extent that time evolution is ergodic; this makes it possible for the Krylov
basis to span the entire Hilbert space. However, it would be interesting to
understand the role of chaos v. integrability better. It would seem that in a
free theory, the original description in terms of the elementary fields of the
theory is already classical (at least for the time evolution of Gaussian
states), so from a computational standpoint, there is no reason to look for a
more efficient description.

  6. 6.

There is more to AdS/CFT than just Hamiltonian evolution – one also needs to
be able to interpret operator insertions in the boundary and the corresponding
matter fields in the bulk within this framework. When Hamiltonian evolution is
ergodic and spans the entire Hilbert space, the dual description presumably
involves pure gravity. However, when this fails, we need additional operators
(other than the Hamiltonian) in order to construct a basis for the Hilbert
space. This should clearly result in a richer structure involving bulk matter.

  7. 7.

There is a beautiful picture for the emergence of spacetime and gravity in
holography from entanglement in the dual CFT VanRaamsdonk:2010pw ;
Lashkari:2013koa ; Faulkner:2013ica ; Faulkner:2017tkh ; Lewkowycz:2018sgn
based on the Ryu-Takayanagi formula Ryu:2006bv . It would be nice to
understand how that picture fits with the complexity-based picture advocated
here.

  8. 8.

Finally, this discussion seems very tailored to $0+1$-dimensional quantum
systems, and it would be interesting to explore whether gravity in higher
dimensions could emerge from similar considerations.

###  5.2 Other future directions

|
---|---
Figure 4: Temperature dependence of the negativity of the KW function for $D=251$. |
---|---
Figure 5: Wigner negativity for the infinite temperature thermofield double
state for various values of $D$.

It would be useful to relate the Wigner negativity to the other notions of
computational complexity developed in nielsen2005geometric ; Nielsen_2006 ;
Jefferson_2017 ; Chapman_2018 ; Brown_2018 ; Balasubramanian:2019wgd ;
Balasubramanian:2021mxo . The quantity of interest in these works was the
geodesic distance between an initial state/unitary and a final state/unitary
with respect to a certain choice of complexity metric. The metric is such that
motion along local directions is taken to be “cheap”, while motion along non-
local directions is penalized with a large cost factor. The Wigner negativity
is not quite the same as circuit complexity, but one essential similarity is
the classification of Clifford operations as being cheap and non-Clifford
operations as being costly. One could imagine constructing a model of circuit
complexity where one regards Clifford directions as being cost-free, and the
remaining directions to have unit cost. It would be interesting to explore
such notions of complexity and their connection with Wigner negativity and
other magic monotones.

In this paper, we have focused on the time-dependence of the KW-negativity in
random matrix theory with the initial state being a generic pure state at
infinite temperature. We can include finite temperature by taking the initial
state to be of the form:

$|\psi_{0}\rangle=\frac{1}{\sqrt{Z_{\beta}}}e^{-\frac{\beta}{2}H}|0\rangle,$
(84)

where $\beta$ is the inverse temperature. For inverse temperature $\beta\sim
1$, the time-dependence of the negativity is not affected very much (see
figure 4 for some preliminary results). In particular, it shows the same
spread, ramp and plateau behaviour, with a very similar saturation value.
However, as the temperature is lowered, the saturation value of the negativity
decreases (although it remains proportional to $\sqrt{D}$). The distinction
between the spread and the ramp phases also seems to get washed out at smaller
temperatures, although this could be an artefact of working at low Hilbert
space dimension. It would be interesting to explore the temperature dependence
of the negativity in more detail. It would also be interesting to consider
other initial states, such as for example the thermofield double state (see
figure 5 for some preliminary results in the infinite temperature case). Here,
the saturation time seems to be $O(D)$, as opposed to the $O(\sqrt{D})$
scaling of the saturation time for pure initial states. Finally, it would be
very interesting to explore these ideas in toy models for holography such as
the SYK or the DSSYK model.

### Acknowledgements

We thank Sujay Ashok, Vijay Balasubramanian, Pawel Caputa, Jackson Fliss,
Abhijit Gadde, Prahladh Harsha, Rohit Kalloor, Gautam Mandal, Shiraz Minwalla,
Harshit Rajgadia and Sandip Trivedi for useful discussions and comments on an
earlier version of the draft. We are grateful to Pruthvi Suryadevara for his
significant help with numerical computations. We acknowledge supported from
the Department of Atomic Energy, Government of India, under project
identification number RTI 4002. We are grateful to the long term workshop
YITP-T-23-01 held at YITP, Kyoto University, where part of this work was done.

##  Appendix A Early time approximation for the Wigner function

In order to evaluate the KW function:

$W(Q,P)=\frac{1}{D}\sum_{k,\ell}\widehat{\delta}_{2Q,k+\ell}e^{\frac{2\pi
i}{D}(k-\ell)P}\langle k|e^{-itH}|0\rangle\langle 0|e^{itH}|\ell\rangle,$ (85)

we need to evaluate the matrix elements $\langle k|e^{-itH}|0\rangle$. For
$n\ll D$, the Krylov coefficients in GUE are approximately given by $b_{n}\sim
1$ and $a_{n}\sim 0$ Balasubramanian:2022tpr . With this approximation, the
Hamiltonian acting on the Krylov basis gives:

$H|n\rangle=|n+1\rangle+|n-1\rangle,\;\;\;\cdots\;\;\;(n\ll D).$ (86)

So, for $Q\ll\frac{D-1}{2}$ and small times, we can evaluate the matrix
elements in equation (85) by treating the system as a simple harmonic
oscillator. Consider the operators

$E_{+}=\sum_{n=0}^{\infty}|n+1\rangle\langle
n|,\;\;E_{-}=\sum_{n=1}^{\infty}|n-1\rangle\langle n|.$ (87)

in the simple harmonic oscillator. Then, the action of our Hamiltonian on
states with $n\ll D$ takes the form

$H\sim E_{+}+E_{-}.$ (88)

This description is, of course, approximate. Firstly, when $n$ starts scaling
with $D$, the Krylov coefficients deviate significantly from 1. Secondly, the
fact that the Hilbert space has $D$ states and is not infinite dimensional (as
in the harmonic oscillator) becomes important. Nevertheless, we expect the
approximation to be good for small times. In the harmonic oscillator, the
eigenstates of the operator $(E_{+}+E_{-})$ are known SusskindGlogower ; phase
(see also Appendix E of Rabinovici:2023yex ):

$|\theta\rangle=\sqrt{\frac{2}{\pi}}\sum_{n=0}^{\infty}\sin[(n+1)\theta]|n\rangle,$
(89) $(E_{+}+E_{-})|\theta\rangle=2\cos\,\theta|\theta\rangle.$ (90)

These states satisfy the completeness relation:

$\int_{0}^{\pi}d\theta\,|\theta\rangle\langle\theta|=1.$ (91)

Figure 6: A comparison of the Wigner negativity computed from the approximate
formula (95) and the numerical result for D = 151.

We can use these relations to approximately evaluate the matrix elements of
interest to us:

$\langle k|e^{-itH}|0\rangle\sim\int_{0}^{\pi}d\theta\,\langle
k|\theta\rangle\langle\theta|e^{-itH}|0\rangle=\frac{2}{\pi}\int_{0}^{\pi}d\theta\,e^{-2it\cos\theta}\sin[(k+1)\theta]\sin(\theta).$
(92)

Thus, the Wigner function in this approximation is given by

$W(Q,P)\sim\frac{4}{D\pi^{2}}\sum_{k,\ell}\delta_{2Q,k+\ell}e^{\frac{2\pi
i}{D}(k-\ell)P}f_{k}(-t)f_{\ell}(t)$ (93)

where

$\displaystyle f_{k}(t)$ $\displaystyle=$
$\displaystyle\int_{0}^{\pi}d\theta\,\sin\theta\sin[(k+1)\theta]e^{2it\cos\theta}$
(94) $\displaystyle=$
$\displaystyle-\frac{1}{2}\int_{0}^{\pi}d\theta\,\left(\cos[(k+2)\theta]-\cos[k\theta]\right)e^{2it\cos\theta}$
$\displaystyle=$ $\displaystyle\frac{\pi
i^{k}}{2}\left(J_{k+2}(2t)+J_{k}(2t)\right).$

Thus,

$W(Q,P)\sim\frac{i^{2Q}}{D}\sum_{k,\ell}\widehat{\delta}_{2Q,k+\ell}e^{\frac{2\pi
i}{D}(k-\ell)P}\left(J_{k+2}(-2t)+J_{k}(-2t)\right)\left(J_{\ell+2}(2t)+J_{\ell}(2t)\right).$
(95)

This is a universal formula for the early time behaviour of the KW function in
GUE. The negativity computed from this formula agrees well with the numerical
computation at early times (see figure 6).

## References

  * (1) D. Gottesman, _The heisenberg representation of quantum computers_ , 1998.
  * (2) S. Aaronson and D. Gottesman, _Improved simulation of stabilizer circuits_ , _Phys. Rev. A_ 70 (2004) 052328 [quant-ph/0406196].
  * (3) A. Mari and J. Eisert, _Positive wigner functions render classical simulation of quantum computation efficient_ , _Physical Review Letters_ 109 (2012) .
  * (4) V. Veitch, C. Ferrie, D. Gross and J. Emerson, _Negative quasi-probability as a resource for quantum computation_ , _New Journal of Physics_ 14 (2012) 113011.
  * (5) H. Pashayan, J.J. Wallman and S.D. Bartlett, _Estimating outcome probabilities of quantum circuits using quasiprobabilities_ , _Physical Review Letters_ 115 (2015) .
  * (6) X. Wang, M.M. Wilde and Y. Su, _Quantifying the magic of quantum channels_ , _New Journal of Physics_ 21 (2019) 103002.
  * (7) E. Wigner, _On the quantum correction for thermodynamic equilibrium_ , _Phys. Rev._ 40 (1932) 749.
  * (8) A. Dhar, G. Mandal and S.R. Wadia, _Classical Fermi fluid and geometric action for c=1_ , _Int. J. Mod. Phys. A_ 8 (1993) 325 [hep-th/9204028].
  * (9) A. Dhar, G. Mandal and S.R. Wadia, _Nonrelativistic fermions, coadjoint orbits of W(infinity) and string field theory at c = 1_ , _Mod. Phys. Lett. A_ 7 (1992) 3129 [hep-th/9207011].
  * (10) H. Lin, O. Lunin and J.M. Maldacena, _Bubbling AdS space and 1/2 BPS geometries_ , _JHEP_ 10 (2004) 025 [hep-th/0409174].
  * (11) V. Balasubramanian, J. de Boer, V. Jejjala and J. Simon, _The Library of Babel: On the origin of gravitational thermodynamics_ , _JHEP_ 12 (2005) 006 [hep-th/0508023].
  * (12) W.K. Wootters, _A wigner-function formulation of finite-state quantum mechanics_ , _Annals of Physics_ 176 (1987) 1.
  * (13) U. Leonhardt, _Quantum-state tomography and discrete wigner function_ , _Phys. Rev. Lett._ 74 (1995) 4101.
  * (14) S. Heiss and S. Weigert, _Discrete moyal-type representations for a spin_ , _Phys. Rev. A_ 63 (2000) 012105.
  * (15) C. Miquel, J.P. Paz and M. Saraceno, _Quantum computers in phase space_ , _Phys. Rev. A_ 65 (2002) 062309.
  * (16) K.S. Gibbons, M.J. Hoffman and W.K. Wootters, _Discrete phase space based on finite fields_ , _Physical Review A_ 70 (2004) .
  * (17) D. Gross, _Hudson’s theorem for finite-dimensional quantum systems_ , _Journal of Mathematical Physics_ 47 (2006) .
  * (18) C. Cormick, E.F. Galvão, D. Gottesman, J.P. Paz and A.O. Pittenger, _Classicality in discrete wigner functions_ , _Phys. Rev. A_ 73 (2006) 012301.
  * (19) V. Veitch, S.A.H. Mousavian, D. Gottesman and J. Emerson, _The resource theory of stabilizer quantum computation_ , _New Journal of Physics_ 16 (2014) 013009.
  * (20) C.D. White, C. Cao and B. Swingle, _Conformal field theories are magical_ , _Phys. Rev. B_ 103 (2021) 075145 [2007.01303].
  * (21) Z.-W. Liu and A. Winter, _Many-Body Quantum Magic_ , _PRX Quantum_ 3 (2022) 020333 [2010.13817].
  * (22) X. Turkeshi, A. Dymarsky and P. Sierant, _Pauli Spectrum and Magic of Typical Quantum Many-Body States_ , 2312.11631.
  * (23) P. Niroula, C.D. White, Q. Wang, S. Johri, D. Zhu, C. Monroe et al., _Phase transition in magic with random quantum circuits_ , 2304.10481.
  * (24) J.R. Fliss, _Knots, links, and long-range magic_ , _JHEP_ 04 (2021) 090 [2011.01962].
  * (25) L. Leone, S.F.E. Oliviero, Y. Zhou and A. Hamma, _Quantum Chaos is Quantum_ , _Quantum_ 5 (2021) 453.
  * (26) S.F. Oliviero, L. Leone and A. Hamma, _Transitions in entanglement complexity in random quantum circuits by measurements_ , _Physics Letters A_ 418 (2021) 127721.
  * (27) K. Goto, T. Nosaka and M. Nozaki, _Probing chaos by magic monotones_ , _Phys. Rev. D_ 106 (2022) 126009.
  * (28) V. Balasubramanian, P. Caputa, J.M. Magan and Q. Wu, _Quantum chaos and the complexity of spread of states_ , _Phys. Rev. D_ 106 (2022) 046007 [2202.06957].
  * (29) D.E. Parker, X. Cao, A. Avdoshkin, T. Scaffidi and E. Altman, _A Universal Operator Growth Hypothesis_ , _Phys. Rev. X_ 9 (2019) 041017 [1812.08657].
  * (30) S.-K. Jian, B. Swingle and Z.-Y. Xian, _Complexity growth of operators in the SYK model and in JT gravity_ , _JHEP_ 03 (2021) 014 [2008.12274].
  * (31) E. Rabinovici, A. Sánchez-Garrido, R. Shir and J. Sonner, _Operator complexity: a journey to the edge of krylov space_ , _Journal of High Energy Physics_ 2021 (2021) .
  * (32) P. Caputa, J.M. Magan and D. Patramanis, _Geometry of krylov complexity_ , 2021\.
  * (33) M.A. Nielsen, _A geometric approach to quantum circuit lower bounds_ , 2005\.
  * (34) M.A. Nielsen, M.R. Dowling, M. Gu and A.C. Doherty, _Quantum computation as geometry_ , _Science_ 311 (2006) 1133–1135.
  * (35) R.A. Jefferson and R.C. Myers, _Circuit complexity in quantum field theory_ , _Journal of High Energy Physics_ 2017 (2017) .
  * (36) S. Chapman, M.P. Heller, H. Marrochio and F. Pastawski, _Toward a definition of complexity for quantum field theory states_ , _Physical Review Letters_ 120 (2018) .
  * (37) A.R. Brown and L. Susskind, _Second law of quantum complexity_ , _Physical Review D_ 97 (2018) .
  * (38) V. Balasubramanian, M. Decross, A. Kar and O. Parrikar, _Quantum Complexity of Time Evolution with Chaotic Hamiltonians_ , _JHEP_ 01 (2020) 134 [1905.05765].
  * (39) V. Balasubramanian, M. DeCross, A. Kar, Y.C. Li and O. Parrikar, _Complexity growth in integrable and chaotic models_ , _JHEP_ 07 (2021) 011 [2101.02209].
  * (40) B. Craps, O. Evnin and G. Pascuzzi, _A relation between Krylov and Nielsen complexity_ , 2311.18401.
  * (41) J.S. Cotler, G. Gur-Ari, M. Hanada, J. Polchinski, P. Saad, S.H. Shenker et al., _Black Holes and Random Matrices_ , _JHEP_ 05 (2017) 118 [1611.04650].
  * (42) M. Berkooz, M. Isachenkov, V. Narovlansky and G. Torrents, _Towards a full solution of the large N double-scaled SYK model_ , _JHEP_ 03 (2019) 079 [1811.02584].
  * (43) H.W. Lin, _The bulk Hilbert space of double scaled SYK_ , _JHEP_ 11 (2022) 060 [2208.07032].
  * (44) J. Maldacena, D. Stanford and Z. Yang, _Conformal symmetry and its breaking in two dimensional nearly anti-de-sitter space_ , 2016.
  * (45) E. Rabinovici, A. Sánchez-Garrido, R. Shir and J. Sonner, _A bulk manifestation of Krylov complexity_ , _JHEP_ 08 (2023) 213 [2305.04355].
  * (46) L. Susskind, _Computational complexity and black hole horizons_ , 2014.
  * (47) L. Susskind, _Entanglement is not enough_ , _Fortsch. Phys._ 64 (2016) 49 [1411.0690].
  * (48) D. Stanford and L. Susskind, _Complexity and Shock Wave Geometries_ , _Phys. Rev. D_ 90 (2014) 126007 [1406.2678].
  * (49) B. Czech, _Einstein Equations from Varying Complexity_ , _Phys. Rev. Lett._ 120 (2018) 031601 [1706.00965].
  * (50) P. Caputa and J.M. Magan, _Quantum Computation as Gravity_ , _Phys. Rev. Lett._ 122 (2019) 231302 [1807.04422].
  * (51) P. Saad, S.H. Shenker and D. Stanford, _JT gravity as a matrix integral_ , 1903.11115.
  * (52) A. Kenfack and K. Życzkowski, _Negativity of the wigner function as an indicator of non-classicality_ , _Journal of Optics B: Quantum and Semiclassical Optics_ 6 (2004) 396.
  * (53) R. Hudson, _When is the wigner quasi-probability density non-negative?_ , _Reports on Mathematical Physics_ 6 (1974) 249.
  * (54) F. Soto and P. Claverie, _When is the Wigner function of multidimensional systems nonnegative?_ , _Journal of Mathematical Physics_ 24 (1983) 97 [https://pubs.aip.org/aip/jmp/article-pdf/24/1/97/11216125/97_1_online.pdf].
  * (55) F. Soto‐Eguibar and P. Claverie, _Time evolution of the Wigner function_ , _Journal of Mathematical Physics_ 24 (1983) 1104 [https://pubs.aip.org/aip/jmp/article-pdf/24/5/1104/8153406/1104_1_online.pdf].
  * (56) S.D. Bartlett, B.C. Sanders, S.L. Braunstein and K. Nemoto, _Efficient classical simulation of continuous variable quantum information processes_ , _Phys. Rev. Lett._ 88 (2002) 097904.
  * (57) S. Bravyi and A. Kitaev, _Universal quantum computation with ideal clifford gates and noisy ancillas_ , _Physical Review A_ 71 (2005) .
  * (58) I. Kourkoulou and J. Maldacena, _Pure states in the SYK model and nearly- $AdS_{2}$ gravity_, 1707.02325.
  * (59) V. Balasubramanian, J.M. Magan and Q. Wu, _Tridiagonalizing random matrices_ , _Phys. Rev. D_ 107 (2023) 126001 [2208.08452].
  * (60) L.G. Yaffe, _Large $n$ limits as classical mechanics_, _Rev. Mod. Phys._ 54 (1982) 407.
  * (61) M. Van Raamsdonk, _Building up spacetime with quantum entanglement_ , _Gen. Rel. Grav._ 42 (2010) 2323 [1005.3035].
  * (62) N. Lashkari, M.B. McDermott and M. Van Raamsdonk, _Gravitational dynamics from entanglement ’thermodynamics’_ , _JHEP_ 04 (2014) 195 [1308.3716].
  * (63) T. Faulkner, M. Guica, T. Hartman, R.C. Myers and M. Van Raamsdonk, _Gravitation from Entanglement in Holographic CFTs_ , _JHEP_ 03 (2014) 051 [1312.7856].
  * (64) T. Faulkner, F.M. Haehl, E. Hijano, O. Parrikar, C. Rabideau and M. Van Raamsdonk, _Nonlinear Gravity from Entanglement in Conformal Field Theories_ , _JHEP_ 08 (2017) 057 [1705.03026].
  * (65) A. Lewkowycz and O. Parrikar, _The holographic shape of entanglement and Einstein’s equations_ , _JHEP_ 05 (2018) 147 [1802.10103].
  * (66) S. Ryu and T. Takayanagi, _Holographic derivation of entanglement entropy from AdS/CFT_ , _Phys. Rev. Lett._ 96 (2006) 181602 [hep-th/0603001].
  * (67) L. Susskind and J. Glogower, _Quantum mechanical phase and time operator_ , _Physics Physique Fizika_ 1 (1964) 49.
  * (68) P. CARRUTHERS and M.M. NIETO, _Phase and angle variables in quantum mechanics_ , _Rev. Mod. Phys._ 40 (1968) 411.