# Limiting spectral distribution of random self-adjoint quantum channels

Cécilia Lancien and Patrick Oliveira Santos and Pierre Youssef  Cécilia
Lancien. CNRS & Institut Fourier UMR 5582, Université Grenoble Alpes,
Grenoble, France<EMAIL_ADDRESS>Patrick Oliveira
Santos. Laboratoire d’Analyse et de Mathématiques Appliquées UMR 8050,
Université Gustave Eiffel, Université Paris Est Créteil, Marne-la-Vallée,
France<EMAIL_ADDRESS>Pierre Youssef. Division of Science,
NYU Abu Dhabi, Abu Dhabi, UAE & Courant Institute of Mathematical Sciences,
New York University, New York, USA<EMAIL_ADDRESS>
###### Abstract.

We study the limiting spectral distribution of quantum channels whose Kraus
operators are sampled as $n\times n$ random Hermitian matrices satisfying
certain assumptions. We show that when the Kraus rank goes to infinity with
$n$, the limiting spectral distribution (suitably rescaled) of the
corresponding quantum channel coincides with the semi-circle distribution.
When the Kraus rank is fixed, the limiting spectral distribution is no longer
the semi-circle distribution. It corresponds to an explicit law, which can
also be described using tools from free probability.

##  1\. Introduction

In quantum physics, the state of an $n$-dimensional system is described by a
density operator on $\mathbb{C}^{n}$, i.e. $\rho\in M_{n}(\mathbb{C})$ a
Hermitian positive semidefinite matrix with trace $1$: $\rho\succeq 0$ and
$\operatorname{tr}(\rho)=1$. A transformation of such a quantum system is
described by a quantum channel on $M_{n}(\mathbb{C})$, i.e.
$\Phi:M_{n}(\mathbb{C})\to M_{n}(\mathbb{C})$ a completely positive and trace-
preserving linear map. We recall that a linear map $\Phi:M_{n}(\mathbb{C})\to
M_{n}(\mathbb{C})$ is said to be

  * •

positive if it preserves the fact of being Hermitian and positive
semidefinite: for all $X\in M_{n}(\mathbb{C})$, $X\succeq 0$ implies that
$\Phi(X)\succeq 0$;

  * •

completely positive if $\Phi\otimes\operatorname{Id}:M_{n^{2}}(\mathbb{C})\to
M_{n^{2}}(\mathbb{C})$ is positive (where
$\operatorname{Id}:M_{n}(\mathbb{C})\to M_{n}(\mathbb{C})$ denotes the
identity map);

  * •

trace-preserving if it preserves the trace: for all $X\in M_{n}(\mathbb{C})$,
$\operatorname{tr}(\Phi(X))=\operatorname{tr}(X)$.

A quantum channel $\Phi$ thus maps quantum states to quantum states (and so
does $\Phi\otimes\mathrm{id}$).

The action of a completely positive map $\Phi$ on $M_{n}(\mathbb{C})$ can
always be described in the following (non-unique) way, called a Kraus
representation of $\Phi$ (see e.g. [3, Section 2.3.2] or [25, Chapter 2]):
There exist $d\in\mathbb{N}$ and $K_{1},\ldots,K_{d}\in M_{n}(\mathbb{C})$,
called Kraus operators of $\Phi$, such that

$\Phi:X\in M_{n}(\mathbb{C})\mapsto\sum_{i\in[d]}K_{i}XK_{i}^{*}\in
M_{n}(\mathbb{C}),$ (1)

where $K_{i}^{*}$ stands for the adjoint of $K_{i}$ (and where we have used
the shorthand notation $[d]$ for the set of integers $\\{1,\ldots,d\\}$). The
fact that $\Phi$ is trace-preserving is equivalent to the following constraint
on the Kraus operators $K_{1},\ldots,K_{d}$:

$\sum_{i\in[d]}K_{i}^{*}K_{i}=\operatorname{Id}.$

The smallest $d$ such that an expression of the form of equation (1) for
$\Phi$ exists is called the Kraus rank of $\Phi$. It is always at most $n^{2}$
for a completely positive map on $M_{n}(\mathbb{C})$.

Given a completely positive map $\Phi$ on $M_{n}(\mathbb{C})$, its adjoint (or
dual) is defined by duality with respect to the Hilbert-Schmidt inner product,
i.e. as the map $\Phi^{*}$ on $M_{n}(\mathbb{C})$ (which is completely
positive as well) such that, for all $X,Y\in M_{n}(\mathbb{C})$,

$\operatorname{tr}\left(X\Phi^{*}(Y)\right)=\operatorname{tr}\left(\Phi(X)Y\right).$
(2)

$\Phi$ being trace-preserving is equivalent to $\Phi^{*}$ being unital, i.e.
such that $\Phi^{*}(\operatorname{Id})=\operatorname{Id}$.

Note that, identifying $M_{n}(\mathbb{C})$ with
$\mathbb{C}^{n}\otimes\mathbb{C}^{n}$, a linear map $\Phi:M_{n}(\mathbb{C})\to
M_{n}(\mathbb{C})$ can equivalently be seen as a linear map
$M_{\Phi}:\mathbb{C}^{n}\otimes\mathbb{C}^{n}\to\mathbb{C}^{n}\otimes\mathbb{C}^{n}$,
i.e. an element of $M_{n^{2}}(\mathbb{C})$. Concretely, a completely positive
linear map

$\Phi:X\in M_{n}(\mathbb{C})\mapsto\sum_{i\in[d]}K_{i}XK_{i}^{*}\in
M_{n}(\mathbb{C})$

can be identified with

$M_{\Phi}=\sum_{i\in[d]}K_{i}\otimes\overline{K}_{i}\in
M_{n^{2}}(\mathbb{C}),$ (3)

where $\overline{K}_{i}$ stands for the entry-wise conjugate of $K_{i}$, in
the canonical basis of $\mathbb{C}^{n}$. This identification preserves the
spectrum, i.e. $\mathrm{spec}(\Phi)=\mathrm{spec}(M_{\Phi})$. Moreover, the
matrix version of the adjoint $\Phi^{*}$ of $\Phi$ is simply the adjoint of
$M_{\Phi}$, i.e. $M_{\Phi^{*}}=M_{\Phi}^{*}$.

In this paper, we will consider quantum channels $\Phi$ whose Kraus operators
are Hermitian matrices, which ensures that $\Phi$ itself is Hermitian (in the
sense that it is equal to its adjoint $\Phi^{*}$, as defined by equation (2)),
or equivalently that its matrix version $M_{\Phi}$ is Hermitian. This may look
like a restrictive setting, but it, in fact, encompasses all Hermitian quantum
channels. Indeed, it is not hard to see that any Hermitian completely positive
map $\Phi$ on $M_{n}(\mathbb{C})$ admits a Kraus representation with Hermitian
Kraus operators, where we can additionally guarantee that the number of such
operators is at most twice the Kraus rank of $\Phi$ (i.e. in general at most
$2n^{2}$). Concretely, if $K_{1},\ldots,K_{d}\in M_{n}(\mathbb{C})$ are $d$
Kraus operators for $\Phi$, then $K_{1,R},K_{1,I},\ldots,K_{d,R},K_{d,I}\in
M_{n}(\mathbb{C})$ are $2d$ Hermitian Kraus operators for $\Phi$, where given
$K\in M_{n}(\mathbb{C})$, we set $K_{R}=(K+K^{*})/2,\,K_{I}=-i(K-K^{*})/2$, so
that $K=K_{R}+iK_{I}$ and $K_{R}^{*}=K_{R},\,K_{I}^{*}=K_{I}$. The latter
claim follows from the observation that,if $M_{\Phi}^{*}=M_{\Phi}$, then we
can re-write $M_{\Phi}=(M_{\Phi}+M_{\Phi}^{*})/2$, and thus

$M_{\Phi}=\frac{1}{2}\sum_{j\in[d]}\left(K_{j}\otimes\overline{K}_{j}+K_{j}^{*}\otimes\overline{K}_{j}^{*}\right)=\sum_{j\in[d]}\left(K_{j,R}\otimes\overline{K}_{j,R}+K_{j,I}\otimes\overline{K}_{j,I}\right).$

Now, many interesting and well-studied classes of quantum channels are
Hermitian: depolarizing channels, dephasing channels, and Pauli channels, to
name just a few. Note that such channels are, in particular, unital.

This work aims to study the spectrum of a randomly generated quantum channel
as the underlying dimension $n$ goes to infinity. More precisely, given a
quantum channel $\Phi$ on $M_{n}(\mathbb{C})$ whose Kraus operators
$K_{1},\ldots,K_{d}\in M_{n}(\mathbb{C})$ have been randomly sampled, we aim
to characterize its asymptotic spectrum (i.e. its spectrum in the limit where
$n$ goes to infinity). Our study will involve two regimes, one where the Kraus
rank $d$ is fixed and only $n$ grows, and the other one where $d=d(n)$ is a
growing function of $n$.

Previous related works were mostly concerned with identifying the spectral gap
of a random quantum channel, i.e. the difference between its largest and
second largest eigenvalues. Indeed, it is known that a quantum channel $\Phi$
on $M_{n}(\mathbb{C})$ always has its largest (in modulus) eigenvalue
$\lambda_{1}(\Phi)$ equal to $1$ (with an associated eigenvector which is a
positive semidefinite matrix), implying that $\Phi$ always has a fixed state.
It was established that random quantum channels generically have their largest
eigenvalue $1$ isolated from the rest of the spectrum. All other eigenvalues
are of order at most $1/\sqrt{d}$ (in modulus). This was shown first for
specific models, where Kraus operators were sampled either as independent Haar
unitaries [16, 22], or as blocks of a Haar isometry [15], or as independent
Ginibre matrices (i.e. matrices having i.i.d. Gaussian entries) [18], and
recently in greater generality [19].

On the other hand, much less is known concerning the asymptotic distribution
of the bulk of the spectrum. To understand this, one studies the empirical
spectral distribution, which, given a matrix $M\in M_{n}(\mathbb{C})$, is
defined as

$\displaystyle\mu_{M}:=\frac{1}{n}\sum_{k\in[n]}\delta_{\lambda_{k}(M)},$

where $\lambda_{k}(M)$ is the $k$th eigenvalue of $M$. In [8], a model where
Kraus operators are sampled as blocks of a Haar isometry was studied. It was
conjectured, from heuristic arguments and numerical simulations, that in the
regime where $d=d(n)=n^{2}$, the empirical spectral distribution (rescaled by
a factor $\sqrt{d}$) of such random quantum channel converges towards a
circular distribution as $n$ grows. It was later shown in [2] that, in the
same regime $d=d(n)$ of order $n^{2}$, the empirical singular value
distribution (again rescaled by a factor $\sqrt{d}$) of a random quantum
channel whose Kraus operators are sampled as independent Ginibre matrices
converges towards a quarter-circular distribution as $n$ grows. We believe
that for both models (and for more general ones) and in any regime
$d=d(n)\to\infty$, the limiting spectral distribution coincides with the
circular law.

In this work, we embark on the investigation of the limiting spectral
distribution of quantum channels in the Hermitian setting. This serves as a
precursor to our broader research goals involving the non-Hermitian case and
the conjecture presented above. Classical problems in Random Matrix Theory
have traditionally prioritized the exploration of Hermitian matrices before
delving into their non-Hermitian counterparts, partly due to the inherent
technical complexities introduced by non-Hermitian systems, notably the
instability of the spectrum under perturbations (see [4, Chapter 11]). For
instance, while the limiting spectral distribution of an $n\times n$ Hermitian
matrix with i.i.d. (up to symmetry) centered entries of variance $1/n$ was
shown to be Wigner’s semi-circle distribution, it is almost half a century
later that the analogous result in the non-Hermitian case, Girko’s circular
law theorem, was established in full generality [23]. Understanding the
limiting spectral distribution in the non-Hermitian case follows the
Hermitization technique invented by Girko [13, 14]. The latter requires a
quantitative control on the smallest singular value of the corresponding
random matrix model (see [6] for an introduction to the method), making the
problem significantly more involved. In the context of quantum channels, the
associated random matrix model exhibits dependencies among its entries, which
adds a layer of difficulty compared to classical random matrix models, even in
the Hermitian setting.

Now, let us shift our focus back to the specific context of this paper. Here,
we will consider the case where the Kraus operators of the quantum channel
$\Phi$ are chosen to be random Hermitian operators. This ensures that the
resulting completely positive map $\Phi$, or equivalently its matrix version
$M_{\Phi}$, is Hermitian. Conversely, as explained earlier, any Hermitian
completely positive map $\Phi$ can be written with Hermitian Kraus operators.
Moreover, we aim to keep our assumptions regarding the distribution of these
random Kraus operators as minimal as possible. Our objective is to gain a
comprehensive understanding of the spectrum of $M_{\Phi}$ as the dimension $n$
increases without imposing specific constraints on the scaling relationship
between $d$ and $n$.

To keep the introduction light, we state an informal version of our main
result and differ from the rigorous statement, requiring notions from free
probability, to Section 2 (see Theorem 2.2).

######  Theorem 1.1 (Informal Statement).

Let $W_{1},\ldots,W_{d}\in M_{n}(\mathbb{C})$ be independent Hermitian random
matrices such that, for each $i\in[d]$, $\operatorname{\mathbb{E}}(W_{i})=0$
and $\operatorname{\mathbb{E}}(W_{i}^{2})=\operatorname{Id}$. Suppose
additionally that, for each $i\in[d]$, the spectral distribution of $W_{i}$
converges to some distribution $\mu_{i}$ as $n\to\infty$ and that the joint
asymptotic distribution of the family $(W_{i})_{i\in[d]}$ is determined by the
family $(\mu_{i})_{i\in[d]}$. Next, set $K_{i}=W_{i}/\sqrt{d}$ for each
$i\in[d]$ and define $\Phi$ as the random completely positive map on
$M_{n}(\mathbb{C})$ having $K_{1},\ldots,K_{d}$ as Kraus operators, i.e.

$\Phi:X\in M_{n}(\mathbb{C})\mapsto\sum_{i\in[d]}K_{i}XK_{i}\in
M_{n}(\mathbb{C}).$

First, $\Phi$ is, on average, trace-preserving and unital. And second,

  * •

If $d$ is fixed, then the spectral distribution of
$\Phi-\operatorname{\mathbb{E}}(\Phi)$ converges to a specific distribution,
depending only on $\mu_{1},\ldots,\mu_{d}$, as $n\to\infty$.

  * •

If $d=d(n)\to\infty$ as $n\to\infty$ and $W_{i}$ are i.i.d, then the spectral
distribution of $\sqrt{d}(\Phi-\operatorname{\mathbb{E}}(\Phi))$ converges to
the semicircular distribution as $n\to\infty$.

The specific assumption underlying the family of independent matrices
$(W_{i})_{i\in[d]}$ is precisely asymptotic freeness, a concept introduced by
Voiculescu [24] within the framework of Free Probability. Asymptotic freeness
proves instrumental in a systematical examination of the asymptotic behavior
of random matrices. Importantly, this assumption is not restrictive as many
classical independent random matrices are known to be asymptotically free.
However, this property is no longer valid in the context of tensors (see
[10]), preventing a direct derivation of the asymptotic joint spectrum of
these random tensors, which is the setting of interest in this paper. To
overcome this, we analyze the moment method directly on our model in order to
derive the limiting distribution.

The required notions from Free Probability are presented in Section 2, where
the main rigorous theorem is stated. The proofs are carried out in Section 3.
Section 4 regroups a variety of examples for which our result applies, while
Section 5 offers perspectives for future work.

### Acknowledgments

Part of this work was completed during a stay of the second named author at
New York University in Abu Dhabi, partly funded by a doctoral mobility grant
delivered by Université Gustave Eiffel; he would like to thank both
institutions for their support and the excellent working assumptions. The
first named author was supported by the ANR projects ESQuisses (grant number
ANR-20-CE47-0014-01), STARS (grant number ANR-20-CE40-0008) and QTraj (grant
number ANR-20-CE40-0024-01).

##  2\. Preliminaries and rigorous statement of the main result

###  2.1. Preliminaries

We begin by recalling some notation from Free Probability; see [21]. A
noncommutative probability space is a pair $(\mathcal{A},\tau)$, where
$\mathcal{A}$ is a unital algebra equipped with a tracial state $\tau$, that
is, $\tau$ is linear, $\tau(1)=1$ and $\tau(ab)=\tau(ba)$. We say that
subalgebras $\mathcal{A}_{1},\ldots,\mathcal{A}_{d}$ are free if

$\displaystyle\tau(a_{1}\ldots a_{p})=0$ (4)

whenever $a_{i}\in\mathcal{A}_{j_{i}}$, $\tau(a_{i})=0$ for all $i\in[p]$ and
$j_{1}\neq j_{2}\neq\cdots\neq j_{p}$. We say that random variables
$a_{1},\ldots,a_{d}\in\mathcal{A}$ are free if their algebras are free. The
distribution of a self-adjoint variable $a=a^{*}$ is defined as the collection
of moments

$\displaystyle\\{\tau(a^{p}):\,p\in\mathbb{N}\\}.$

In particular, there always exists a measure $\mu$ on $\mathbb{R}$ such that

$\displaystyle\tau(a^{p})=\int x^{p}\,\text{d}\mu,$

for every $p\in\mathbb{N}$. Such a measure is also called the distribution of
$a$. A particular example of free variables is the free semicircular system
$(s_{1},\ldots,s_{d})$, whose moments satisfy

$\displaystyle\tau(s_{i_{1}}\cdots s_{i_{p}})=\sum_{\pi\in
NC_{2}(p)}\prod_{(l,k)\in\pi}\delta_{i_{l}i_{k}}.$

Here, $NC_{2}(p)$ is the set of all noncrossing pair partitions of $[p]$,
namely, each block of $\pi\in NC_{2}(p)$ has cardinality $2$ and $\pi$ does
not contain two blocks $\\{i,j\\},\\{k,l\\}$ such that $i<k<j<l$. In
particular, the moments of the semicircular variable are

$\displaystyle\tau(s^{p})=|NC_{2}(p)|=\int x^{p}f_{sc}(x)\,\text{d}x,$

where $f_{sc}$ denotes the density of the semi-circle distribution $\mu_{sc}$

$\displaystyle f_{sc}(x)=\frac{1}{2\pi}\sqrt{4-x^{2}}\mathbf{1}_{|x|\leq 2}.$
(5)

Given a matrix $M\in M_{n}(\mathbb{C})$, we define its normalized trace as

$\displaystyle\tau^{(n)}(M):=\frac{1}{n}\operatorname{tr}(M).$

In particular, for two matrices $M,N\in M_{n}(\mathbb{C})$, we have

$\displaystyle\tau^{(n^{2})}(M\otimes N)=\tau^{(n)}\otimes\tau^{(n)}(M\otimes
N)=\tau^{(n)}(M)\tau^{(n)}(N).$

We say that a random Hermitian matrix $M\in M_{n}(\mathbb{C})$ is normalized
if

$\displaystyle\operatorname{\mathbb{E}}\left(\tau^{(n)}(M)\right)=0;$
$\displaystyle\operatorname{\mathbb{E}}\left(\tau^{(n)}(M^{2})\right)=1.$

We say that a sequence $(M_{n})_{n\in\mathbb{N}}$ of random Hermitian matrices
$M_{n}\in M_{n}(\mathbb{C})$ such that $\operatorname{tr}(M_{n}^{p})$ is
integrable for all $p\geq 1$ converges weakly in probability (resp. in
expectation) to $\mu$, if $(\mu_{M_{n}})_{n\in\mathbb{N}}$ converges weakly to
$\mu$ in probability (resp. in expectation). Equivalently, for any $p\geq 1$,
we have

$\displaystyle\tau^{(n)}(M_{n}^{p})=\frac{1}{n}\text{tr}(M_{n}^{p})\underset{n\to\infty}{\to}\int
x^{p}\,\text{d}\mu\ $

in probability and

$\displaystyle\operatorname{\mathbb{E}}\left(\tau^{(n)}(M_{n}^{p})\right)=\operatorname{\mathbb{E}}\left(\frac{1}{n}\text{tr}(M_{n}^{p})\right)\underset{n\to\infty}{\to}\int
x^{p}\,\text{d}\mu$

for the convergence in expectation. Whenever $(\mu_{M_{n}})_{n\in\mathbb{N}}$
converges weakly to $\mu$, we denote it by $\mu_{M_{n}}\Rightarrow\mu$, or
simply $M_{n}\Rightarrow\mu$ as $n\to\infty$. In particular, we can find a
noncommutative random variable $a$ such that $a$ has distribution $\mu$ and
$M_{n}$ converges to $a$.

Finally, we say that $d$ random Hermitian matrices $M_{1},\ldots,M_{d}\in
M_{n}(\mathbb{C})$ are asymptotically free in probability if

$\displaystyle\lim_{n\to\infty}\tau^{(n)}\Big{(}\left(M_{i_{1}}^{p_{1}}-\tau^{(n)}(M_{i_{1}}^{p_{1}})\right)\cdots\left(M_{i_{m}}^{p_{m}}-\tau^{(n)}(M_{i_{m}}^{p_{m}})\right)\Big{)}=0$

in probability, for all $m\geq 1$, $i_{1}\neq i_{2}\neq\cdots\neq i_{m}\in[d]$
and $p_{1},\ldots,p_{m}\geq 1$. They are asymptotically free in expectation if

$\displaystyle\lim_{n\to\infty}\operatorname{\mathbb{E}}\Big{[}\tau^{(n)}\Big{(}\left(M_{i_{1}}^{p_{1}}-\operatorname{\mathbb{E}}\left(\tau^{(n)}(M_{i_{1}}^{p_{1}})\right)\right)\cdots\left(M_{i_{m}}^{p_{m}}-\operatorname{\mathbb{E}}\left(\tau^{(n)}(M_{i_{m}}^{p_{m}})\right)\right)\Big{)}\Big{]}=0,$

for all $m\geq 1$, $i_{1}\neq i_{2}\neq\cdots\neq i_{m}\in[d]$ and
$p_{1},\ldots,p_{m}\geq 1$. This is equivalent to (4) in the limit for the
algebras $\mathcal{A}(M_{i})$ generated by each $M_{i}$. As usual, we write
$M+\lambda:=M+\lambda\operatorname{Id}$, where $M\in M_{n}(\mathbb{C})$,
$\lambda\in\mathbb{C}$.

######  Remark 2.1.

We have defined weak convergence as the convergence of the moments of all
order. In particular, this requires that $\operatorname{tr}(M_{n}^{p})$ is
integrable for all $p\geq 1$. In several cases, however, such strong
integrability is not needed, and similar results can be proved for a larger
class of random matrices via truncation techniques [1, Theorem 2.1.21].

###  2.2. Rigorous results

Let $\mu_{1},\ldots,\mu_{d}$ be probability measures on $\mathbb{R}$ and
$a_{1},\ldots,a_{d}$ be free variables with distribution $\mu_{i}$,
respectively. We define the tensor measure $\mu_{i}\star\mu_{i}$ as the
measure associated with $a_{i}\otimes a_{i}$ and the tensor convolution
$\mu_{1}\circledast\cdots\circledast\mu_{d}$ as the distribution of

$\displaystyle\sum_{i\in[d]}a_{i}\otimes a_{i}=a_{1}\otimes
a_{1}+\cdots+a_{d}\otimes a_{d}.$

In other words, we have

$\displaystyle\int
x^{k}\,\text{d}\mu_{i}\star\mu_{i}=\tau\otimes\tau(a_{i}^{k}\otimes
a_{i}^{k})=\tau^{2}(a_{i}^{k}),$

and

$\displaystyle\int
x^{k}\,\text{d}\mu_{1}\circledast\cdots\circledast\mu_{d}=\tau\otimes\tau\left(\sum_{i\in[d]}a_{i}\otimes
a_{i}\right)^{k}=\sum_{i\in[d]^{k}}\tau^{2}(a_{i_{1}}\cdots a_{i_{k}}).$

Note that $\mu_{1}\circledast\cdots\circledast\mu_{d}$ is not necessarily the
free convolution of $\mu_{i}\star\mu_{i}$. Indeed, it was proved in [10] that
freeness for tensor products $a_{i}\otimes a_{i}$ does not follow from
freeness of $a_{i}$.

In the sequel, we will be considering random Hermitian matrices
$W_{1},\ldots,W_{d}\in M_{n}(\mathbb{C})$ satisfying the following
assumptions:

  1. (A.1)

For each $i\in[d]$, $W_{i}\in M_{n}(\mathbb{C})$ converges weakly in
probability and in expectation to $\mu_{i}$;

  2. (A.2)

The family $(W_{i})_{i\in[d]}$ is in probability and in expectation
asymptotically free, for each $d$ fixed.

  3. (A.3)

For each $i\in[d]$,
$\operatorname{\mathbb{E}}(W_{i}\otimes\overline{W_{i}})\in
M_{n^{2}}(\mathbb{C})$ converges weakly to $0$.

  4. (A.4)

The matrices $(W_{i})_{i\in[d]}$ are independent and identically distributed.

As we will see in Section 4, a variety of classical random matrix models
satisfy those assumptions. We are ready to state our main result.

######  Theorem 2.2.

Let $W_{1},\ldots,W_{d}\in M_{n}(\mathbb{C})$ be centered Hermitian random
matrices satisfying Assumptions (A.1),(A.2), (A.3) and $d=d(n)$. Then, the
following holds.

  1. (2.2.i)

If $d$ is fixed, then

$\displaystyle\Delta:=\frac{1}{\sqrt{d}}\sum_{i\in[d]}\left(W_{i}\otimes\overline{W}_{i}-\operatorname{\mathbb{E}}(W_{i}\otimes\overline{W}_{i})\right)$

converges weakly in probability and in expectation to the tensor convolution
of $\tilde{\mu}_{1},\ldots,\tilde{\mu}_{d}$, where $\tilde{\mu}$ denotes the
dilation of $\mu$, i.e.

$\displaystyle\int x^{p}\,\text{d}\tilde{\mu}:=\frac{1}{d^{p/2}}\int
x^{p}\,\text{d}\mu.$

This means that

$\displaystyle\mu_{\Delta}\underset{n\to\infty}{\Rightarrow}\tilde{\mu}_{1}\circledast\cdots\circledast\tilde{\mu}_{d},$

in probability and in expectation.

  2. (2.2.ii)

If $d=d(n)$ diverges, $(W_{i})_{i\in[d]}$ are also normalized and satisfy
Assumption (A.4), then $\Delta$ converges weakly in probability and in
expectation to the semicircular distribution $\mu_{sc}$. This means that

$\displaystyle\mu_{\Delta}\underset{n\to\infty}{\Rightarrow}\mu_{sc},$

in probability and in expectation.

##  3\. Proofs

We start with the following observation.

######  Lemma 3.1.

Suppose that Assumptions (A.1) and (A.3) hold. For each $i\in[d]$,
$B_{i}=W_{i}\otimes\overline{W}_{i}-\operatorname{\mathbb{E}}\left(W_{i}\otimes\overline{W}_{i}\right)$
converges weakly in probability and in expectation to $a_{i}\otimes a_{i}$,
where $a_{i}$ has distribution $\mu_{i}$.

###### Proof.

For Hermitian matrices $M_{1},M_{2}\in M_{n}(\mathbb{C})$, such that
$M_{1}\Rightarrow\mu_{1}$ and $M_{2}\Rightarrow 0$ as $n\to\infty$, we have

$\displaystyle M_{1}+M_{2}\Rightarrow\mu_{1}.$

Indeed, using Holder’s Inequality, we can bound the trace of a product of the
matrices $M_{1}$ and $M_{2}$ by

$\displaystyle\left|\tau^{(n)}\left(\prod_{j\in[m]}M_{i_{j}}\right)\right|\leq\prod_{j\in[m]}\left(\tau^{(n)}|M_{i_{j}}|^{p_{j}}\right)^{1/p_{j}},$
(6)

where $i\in\\{1,2\\}^{m}$ and $(p_{j})_{j\in[m]}$ are conjugate exponents. If
there is at least one $j\in[m]$ such that $i_{j}=2$, then the right-hand side
of (6) goes to $0$. In particular, since

$\displaystyle\tau^{(n)}\left((M_{1}+M_{2})^{p}\right)$
$\displaystyle=\sum_{i\in\\{1,2\\}^{p}}\tau^{(n)}\left(\prod_{j\in[p]}M_{i_{j}}\right),$

any term associated with $i\in\\{1,2\\}^{p}$ having at least one index
$i_{j}=2$ will asymptotically vanish. Therefore,

$\displaystyle\lim_{n\to\infty}\tau^{(n)}\left((M_{1}+M_{2})^{p}\right)=\lim_{n\to\infty}\tau^{(n)}\left(M_{1}^{p}\right).$

The result of the lemma follows by taking $M_{1}=W_{i}\otimes\overline{W_{i}}$
and $M_{2}=-\operatorname{\mathbb{E}}(W_{i}\otimes\overline{W_{i}})$, for each
$i\in[d]$. Indeed,

$\displaystyle\tau^{(n)}\otimes\tau^{(n)}\left(\left(W_{i}\otimes\overline{W}_{i}\right)^{p}\right)$
$\displaystyle=\left(\tau^{(n)}\left(W_{i}^{p}\right)\right)^{2},$

since the eigenvalues of $W_{i}$ are real. Hence, the weak convergence in
probability and in expectation of $W_{i}$ to $\mu_{i}$ implies that

$\displaystyle\tau^{(n)}\otimes\tau^{(n)}\left(\left(W_{i}\otimes\overline{W}_{i}\right)^{p}\right)\underset{n\to\infty}{\to}\left(\int
x^{p}\,\text{d}\mu_{i}\right)^{2}$

in probability and in expectation. This finishes the proof. ∎

As an immediate consequence, we have the following.

######  Corollary 3.2.

Suppose that Assumptions (A.1),(A.2), and (A.3) hold. Let
$B_{i}=W_{i}\otimes\overline{W_{i}}-\operatorname{\mathbb{E}}(W_{i}\otimes\overline{W_{i}})$.
Then, for each $d$ fixed, $(B_{i})_{i\in[d]}$ converges weakly in probability
and in expectation to $(a_{i}\otimes a_{i})_{i\in[d]}$, where $a_{i}$ are free
variables with distribution $\mu_{i}$.

###### Proof.

By Assumption (A.2), we can assume $(W_{i})_{i\in[d]}$ converges weakly in
probability and in expectation to $(a_{i})_{i\in[d]}$, where $a_{i}$ are free.
The result follows by Lemma 3.1. ∎

We are ready to prove Theorem (2.2.i).

###### Proof of Theorem (2.2.i).

By Corollary 3.2, we have

$\displaystyle\sum_{i\in[d]}B_{i}\underset{n\to\infty}{\Rightarrow}\sum_{i\in[d]}a_{i}\otimes
a_{i},$

where $a_{i}$ are free with distribution $\mu_{i}$. To conclude, it suffices
to note that $\tilde{\mu}_{i}$ is the distribution of $a_{i}\otimes
a_{i}/\sqrt{d}$. ∎

Before proving Theorem (2.2.ii), let us give a sufficient condition for
Assumption (A.3) to hold.

######  Lemma 3.3.

Let $W_{1},\ldots,W_{d}$ satisfying Assumption (A.1) and assume that

$\displaystyle\operatorname{\mathbb{E}}(\tau(W_{i}))\underset{n\to\infty}{\to}0,$
(7)

for all $i\in[d]$. Let $(W_{i}^{(l)})_{l\in[L]}$ be independent copies of
$W_{i}$. If for each $i\in[d]$ and $L\geq 1$, $(W_{i}^{(l)})_{l\in[L]}$ are
asymptotically free in probability and in expectation, then Assumption (A.3)
holds.

###### Proof.

Fix $i\in[d]$. Then, we can easily see that

$\displaystyle\frac{1}{n^{2}}\operatorname{tr}\left(\left(\operatorname{\mathbb{E}}\left(W_{i}\otimes\overline{W_{i}}\right)\right)^{p}\right)=\frac{1}{n^{2}}\operatorname{\mathbb{E}}\left(\operatorname{tr}^{2}\left(W_{i}^{(1)}\cdots
W_{i}^{(p)}\right)\right).$

Since $(W_{i}^{(l)})_{l\in[p]}$ are asymptotically free in probability and in
expectation, we get that

$\displaystyle\frac{1}{n^{2}}\operatorname{tr}\left(\left(\operatorname{\mathbb{E}}\left(W_{i}\otimes\overline{W_{i}}\right)\right)^{p}\right)\underset{n\to\infty}{\to}\tau^{2}(a_{i}^{(1)}\cdots
a_{i}^{(p)}),$

where $(a_{i}^{(l)})_{l\in[p]}$ are free i.i.d with distribution $\mu_{i}$.
The centered Assumption (7) implies that $\mu_{i}$ is centered, hence
$\tau(a_{i}^{(1)}\cdots a_{i}^{(p)})=0$ by freeness and the result follows. ∎

###  3.1. The asymptotic free Central Limit Theorem

In order to prove Theorem (2.2.ii), we begin with a simple lemma.

######  Lemma 3.4.

Suppose that Assumption (A.4) holds. For each $i\in[d]$, set
$B_{i}:=W_{i}\otimes\overline{W}_{i}-\operatorname{\mathbb{E}}(W_{i}\otimes\overline{W}_{i})$.
Let $(i_{1},\ldots,i_{p})\in[d]^{p}$ and suppose that there exists $k\in[p]$
such that, for all $l\in[p]$ with $l\neq k$, $i_{l}\neq i_{k}$. Then, we have

$\displaystyle\operatorname{\mathbb{E}}\left(\operatorname{tr}(B_{i_{1}}\cdots
B_{i_{p}})\right)=0.$

###### Proof.

We can assume without loss of generality that $k=1$, i.e.
$i_{1}\notin\\{i_{2},\ldots,i_{p}\\}$. In this case, notice that $B_{i_{1}}$
is independent of $B:=B_{i_{2}}\cdots B_{i_{p}}$ by Assumption (A.4), hence

$\displaystyle\operatorname{\mathbb{E}}\left(\operatorname{tr}(B_{i_{1}}B)\right)=\sum_{k_{1},k_{2}\in[n]^{2}}\operatorname{\mathbb{E}}\left(B_{i_{1}}(k_{1},k_{2})\right)\operatorname{\mathbb{E}}\left(B(k_{2},k_{1})\right).$

Since $B_{i_{1}}$ is centered, we have
$\operatorname{\mathbb{E}}(B_{i_{1}}(k_{1},k_{2}))=0$ and the result follows.
∎

Although Lemma 3.4 is a simple consequence of independence and centering, it
is a powerful property that will allow us to obtain a limit distribution for
$\Delta$. We are now ready to prove Theorem (2.2.ii).

###### Proof of Theorem (2.2.ii).

We begin by writing $\Delta$ as

$\displaystyle\Delta=\frac{1}{\sqrt{d}}\sum_{i\in[d]}B_{i},$

where
$B_{i}=W_{i}\otimes\overline{W}_{i}-\operatorname{\mathbb{E}}(W_{i}\otimes\overline{W}_{i})$.
We recall that, given $X\in M_{n^{2}}(\mathbb{C})$ a random matrix, we have

$\displaystyle\operatorname{\mathbb{E}}\left(\tau^{(n^{2})}(X)\right)=\operatorname{\mathbb{E}}\left(\frac{1}{n^{2}}\operatorname{tr}(X)\right).$

Then, we can compute the moments of $\Delta$ by

$\displaystyle\operatorname{\mathbb{E}}\left(\tau^{(n^{2})}(\Delta^{p})\right)=d^{-p/2}\sum_{i\in[d]^{p}}\operatorname{\mathbb{E}}\left(\tau^{(n^{2})}\left(B_{i_{1}}\cdots
B_{i_{p}}\right)\right).$

For every $i\in[d]^{p}$, we associate a partition $\pi(i)$ of $[p]$ defined by
placing any $k,\ell\in[p]$ in the same block of $\pi(i)$ whenever
$i_{k}=i_{l}$. We denote by $P(p)$ the set of all partitions of $[p]$. Notice
that for every $i\in[d]^{p}$,
$\operatorname{\mathbb{E}}\left(\tau^{(n^{2})}(B_{i_{1}}\cdots
B_{i_{p}})\right)$ only depends on the partition $\pi(i)$, since $B_{i}$ are
identically distributed by Assumption (A.4). In particular, let
$\tau^{(n^{2})}(\pi)$ be the common value of
$\operatorname{\mathbb{E}}\left(\tau^{(n^{2})}(B_{i_{1}}\cdots
B_{i_{p}})\right)$ for $\pi(i)=\pi$. Hence,

$\displaystyle\operatorname{\mathbb{E}}\left(\tau^{(n^{2})}(\Delta^{p})\right)=d^{-p/2}\sum_{\pi\in
P(p)}\tau^{(n^{2})}(\pi)\left|\\{i\in[d]^{p}:\pi(i)=\pi\\}\right|.$

The cardinality can be computed by simply choosing an index for each block of
$\pi$, namely,

$\displaystyle\left|\\{i\in[d]^{p}:\pi(i)=\pi\\}\right|=d(d-1)\cdots(d-|\pi|+1)=d^{|\pi|}\left(1+O_{p}\left(\frac{1}{d}\right)\right),$

where $f=O_{p}(1/d)$ means that $f\leq C_{p}/d$, for some constant $C_{p}$
that depends only on $p$. We then have

$\displaystyle\operatorname{\mathbb{E}}\left(\tau^{(n^{2})}(\Delta^{p})\right)=\sum_{\pi\in
P(p)}\tau^{(n^{2})}(\pi)d^{|\pi|-p/2}\left(1+O_{p}\left(\frac{1}{d}\right)\right).$

By Lemma 3.4, if $\pi\in P(p)$ has a block $V$ of size $1$, then we would have
$\tau^{(n^{2})}(\pi)=0$. Hence, we can restrict to partitions without single
blocks, yielding

$\displaystyle\operatorname{\mathbb{E}}\left(\tau^{(n^{2})}(\Delta^{p})\right)=\sum_{\begin{subarray}{c}\pi\in
P(p)\\\ \forall V\in\pi,\,|V|\geq
2\end{subarray}}\tau^{(n^{2})}(\pi)d^{|\pi|-p/2}\left(1+O_{p}\left(\frac{1}{d}\right)\right).$

Using Assumption (A.1) and Holder’s Inequality, we have that

$\displaystyle\tau^{(n^{2})}(\pi)\leq\max_{i\in[d]}\left[\operatorname{\mathbb{E}}\left(\tau^{(n^{2})}(|B_{i}|^{p})\right)\right]^{1/p}\leq
C_{p},$ (8)

where $C_{p}<\infty$ is a constant independent of $n$. Hence, whenever $\pi$
has a block of size $|V|\geq 3$, we have $|\pi|<p/2$ and

$\displaystyle\tau^{(n^{2})}(\pi)d^{|\pi|-p/2}=O_{p}\left(\frac{1}{\sqrt{d}}\right).$

In particular, the only partitions that contribute to the dominating term are
the pair partitions; that is, every block has cardinality two. Let $P_{2}(p)$
be the set of all pair partitions, then

$\displaystyle\operatorname{\mathbb{E}}\left(\tau^{(n^{2})}(\Delta^{p})\right)=\sum_{\pi\in
P_{2}(p)}\tau^{(n^{2})}(\pi)+O_{p}\left(\frac{1}{\sqrt{d}}\right).$

Note that $\tau^{(n^{2})}(\pi)$ depends only on $B_{1},\ldots,B_{p}$.
Corollary 3.2 implies that $(B_{i})_{i\in[d]}$ converges to $(a_{i}\otimes
a_{i})_{i\in[d]}$ in probability and in expectation, where $a_{i}$ are free
variables and the limit of $W_{i}$. Hence,

$\displaystyle\tau^{(n^{2})}(\pi)\underset{n\to\infty}{\to}\tau\otimes\tau(\pi):=\tau^{2}(a_{i_{1}}\cdots
a_{i_{p}}),$

In particular, freeness implies that $\tau\otimes\tau(\pi)=0$ if and only if
$\pi\in P_{2}(p)\setminus NC_{2}(p)$, where $NC_{2}(p)$ is the set of
noncrossing pair partitions. If $\pi\in NC_{2}(p)$, then
$\tau\otimes\tau(\pi)=1$ because the matrices are normalized. We deduce that

$\displaystyle\operatorname{\mathbb{E}}\left(\tau^{(n^{2})}(\Delta^{p})\right)=|NC_{2}(p)|+o_{p}(1)+O_{p}\left(\frac{1}{\sqrt{d}}\right),$

where $f=o_{p}(1)$ denotes a function that depends on $p$ and $f\to 0$ as
$n\to\infty$. This proves the convergence in expectation. To get convergence
in probability, set

$\displaystyle\tau^{(n^{2})}(j)=\tau^{(n^{2})}(B_{j_{1}}\cdots B_{j_{p}}),$

for each $j\in[d]^{p}$. Then,

$\displaystyle\operatorname{var}\left(\tau^{(n^{2})}(\Delta^{p})\right)=\frac{1}{d^{p}}\sum_{i,j\in[d]^{p}}\left\\{\operatorname{\mathbb{E}}\left(\tau^{(n^{2})}(i)\tau^{(n^{2})}(j)\right)-\operatorname{\mathbb{E}}\left(\tau^{(n^{2})}(i)\right)\operatorname{\mathbb{E}}\left(\tau^{(n^{2})}(j)\right)\right\\}.$

Now, each index $i_{l},j_{l}$ must appear at least twice in $i\cup j$, as the
matrices are centered and independent. In particular, the summation has at
most $d^{p}$ indices. We thus have

$\displaystyle\operatorname{var}\left(\tau^{(n^{2})}(\Delta^{p})\right)\leq\max_{i,j\in[d]^{p}}\left\\{\operatorname{\mathbb{E}}\left(\tau^{(n^{2})}(i)\tau^{(n^{2})}(j)\right)-\operatorname{\mathbb{E}}\left(\tau^{(n^{2})}(i)\right)\operatorname{\mathbb{E}}\left(\tau^{(n^{2})}(j)\right)\right\\}.$

Assumption (A.1) implies that $(B_{i})_{i\in[p]}$ are uniformly integrable,
then the Dominated Convergence Theorem [11, Theorem 1.5.8] and Corollary 3.2
imply that

$\displaystyle\operatorname{\mathbb{E}}\left(\tau^{(n^{2})}(i)\tau^{(n^{2})}(j)\right)\underset{n\to\infty}{\to}\tau\otimes\tau(\pi(i))\tau\otimes\tau(\pi(j)),$

The same limit holds for
$\operatorname{\mathbb{E}}\left(\tau^{(n^{2})}(i)\right)\operatorname{\mathbb{E}}\left(\tau^{(n^{2})}(j)\right)$,
and therefore

$\displaystyle\operatorname{var}\left(\tau^{(n^{2})}(\Delta^{p})\right)=o_{p}(1),$

from which convergence in probability follows. ∎

######  Remark 3.5.

As we saw in the proof, the only assumption for the existence of a limit of
$\Delta$ is the centering assumption in Lemma 3.4. In particular, the proof of
Theorem (2.2.ii) works verbatim for noncommutative variables. If
$a_{1},\ldots,a_{d}$ are centered exchangeable variables such that Lemma 3.4
holds, that is, $\tau(a_{i_{1}}\cdots a_{i_{p}})=0$ whenever there exists an
index different than the others, we have

$\displaystyle\lim_{d\to\infty}\tau\left(\frac{1}{\sqrt{d}}\sum_{i\in[d]}a_{i}\right)^{p}=\sum_{\pi\in
P_{2}(p)}\tau(\pi).$

If in addition $(a_{i})_{i\in[d]}=(a_{i}^{(n)})_{i\in[d]}$ depends on some
parameter $n$ and are asymptotically free, as $n$ grows, we get an asymptotic
free Central Limit Theorem

$\displaystyle\lim_{n,d\to\infty}\tau^{(n)}\left(\frac{1}{\sqrt{d}}\sum_{i\in[d]}a_{i}^{(n)}\right)^{p}=|NC_{2}(p)|.$

##  4\. Examples

We begin with an almost surely unital and trace-preserving Hermitian quantum
channel.

######  Example 4.1 (Rademacher diagonal matrices).

Let $R=\mathrm{diag}(\varepsilon_{k})_{k\in[n]}\in M_{n}(\mathbb{C})$ be a
Rademacher diagonal matrix, i.e. $\varepsilon_{1},\ldots,\varepsilon_{n}$ are
independent Rademacher (or symmetric Bernoulli) random variables with
parameter $1/2$. Let $U_{1},\ldots,U_{d}\in M_{n}(\mathbb{C})$ be independent
Haar unitary matrices independent of $R$ as well and set
$W_{i}=U_{i}RU_{i}^{*}$ for each $i\in[d]$. It follows from [24] that
$(W_{i})_{i\in[d]}$ are asymptotically free and Lemma 3.3 yields Assumption
(A.3). Theorem (2.2.i) thus implies that

$\displaystyle\Delta\underset{n\to\infty}{\Rightarrow}\frac{1}{\sqrt{d}}\sum_{i\in[d]}r_{i}\otimes
r_{i},$

in probability and in expectation, where $(r_{i})_{i\in[d]}$ are free with
Rademacher distribution. In particular, for a polynomial

$\displaystyle p(x):=b_{0}+\sum_{l\in[m]}b_{l}x^{l},$

we immediately compute

$\displaystyle p(r_{i})-\tau(p(r_{i}))=\left(\sum_{l\in[\lceil
m/2\rceil]}b_{2l-1}\right)r_{i};$ $\displaystyle p(r_{i}\otimes
r_{i})-\tau\otimes\tau(p(r_{i}\otimes r_{i}))=\left(\sum_{l\in[\lceil
m/2\rceil]}b_{2l-1}\right)r_{i}\otimes r_{i}.$

Since $\tau(r_{i})=0$, we can readily see that the freeness of
$(r_{i})_{i\in[d]}$ is equivalent to the freeness of $(r_{i}\otimes
r_{i})_{i\in[d]}$ in this case. Hence $(r_{i}\otimes
r_{i})_{i\in[d]}\overset{d}{=}(r_{i})_{i\in[d]}$ where the equality holds in
distribution. If we denote the limit of $\Delta$ by $z$, we then have

$\displaystyle z\overset{d}{=}\frac{1}{\sqrt{d}}\sum_{i\in[d]}r_{i}.$

For $d$ even, we can precisely compute this limit. Indeed, [21, Example
12.8.1] shows that the free convolution $(r_{1}+r_{2})$ has arcsine
distribution. Moreover, such a distribution was proved to be the law of
$u+u^{*}$ in [21, Example 1.14], where $u$ is a Haar unitary. We then deduce
that

$\displaystyle
z\overset{d}{=}\frac{1}{\sqrt{d}}\sum_{i\in[d/2]}(u_{i}+u_{i}^{*}).$

The distribution of the free sum of $u_{i}+u_{i}^{*}$ is called the Kesten-
McKay distribution with parameter $d$; see [17],[21, Exercise 12.21]. Its
density is given by

$\displaystyle
f_{KM(d)}(x)=\frac{1}{2\pi}\frac{d}{d^{2}-x^{2}}\sqrt{4(d-1)-x^{2}}\,\mathbf{1}_{|x|\leq
2\sqrt{d-1}}.$

It is also the limit spectral distribution of random $d$-regular graphs [20].
Therefore, the density of $b$ is given by

$\displaystyle\tilde{f}_{KM(d)}(x)=\frac{1}{2\pi}\frac{d}{d-x^{2}}\sqrt{4\left(1-\frac{1}{d}\right)-x^{2}}\,\mathbf{1}_{|x|\leq
2\sqrt{1-\frac{1}{d}}}.$

A direct computation shows that, for all $x\in(-2,2)$,
$\tilde{f}_{KM(d)}(x)\to f_{sc}(x)$ as $d\to\infty$. This is a local version
of Theorem (2.2.ii).

Let us interpret Example 4.1 above in terms of the corresponding random
quantum channel. Given $W=URU^{*}\in M_{n}(\mathbb{C})$ a uniformly rotated
Rademacher matrix, we have $W^{2}=R^{2}=\operatorname{Id}$. Hence, sampling
$W_{1},\ldots,W_{d}\in M_{n}(\mathbb{C})$ independently and uniformly rotated
Rademacher matrices and setting $K_{i}=W_{i}/\sqrt{d}$ for each $i\in[d]$, the
random completely positive map

$\Phi:X\in M_{n}(\mathbb{C})\mapsto\sum_{i\in[d]}K_{i}XK_{i}\in
M_{n}(\mathbb{C})$

is exactly trace-preserving and unital (not just on average), i.e. it is a
random unital quantum channel. More precisely, it is a mixture of unitary
conjugations. What is more, by what precedes, the spectral distribution of
$\sqrt{d}(\Phi-\operatorname{\mathbb{E}}(\Phi))$ converges weakly almost
surely to $\tilde{\mu}_{KM(d)}$ (for $d$ even) as $n\to\infty$. And the
spectral distribution of $\sqrt{d}(\Phi-\operatorname{\mathbb{E}}(\Phi))$
converges weakly in probability and in expectation to $\mu_{sc}$ as
$n,d\to\infty$.

######  Remark 4.2.

There is a straightforward generalization of the previous example to
deterministic matrices. Let $M_{1},\ldots,M_{d}\in M_{n}(\mathbb{C})$ be
deterministic Hermitian matrices such that, for each $i\in[d]$,
$M_{i}\Rightarrow\mu_{i}$ as $n\to\infty$ and $\mu$ is centered. A natural way
to make them asymptotically free is by conjugating them with independent Haar
unitaries. In particular, if we set, for each $i\in[d]$,
$W_{i}=U_{i}M_{i}U_{i}^{*}$, where $U_{1},\ldots,U_{d}\in M_{n}(\mathbb{C})$
are independent uniformly chosen unitaries, then $(W_{i})_{i\in[d]}$ satisfies
Assumptions (A.1),(A.2),(A.3) and (A.4). Indeed, for instance, Assumption
(A.3) follows by Lemma 3.3. So that we can apply Theorem 2.2 to this kind of
matrices.

######  Example 4.3 (Wigner matrices).

A random Hermitian matrix $W=(W_{kl})_{k,l\in[n]}\in M_{n}(\mathbb{C})$ is a
Wigner matrix if $W_{kl}$ are i.i.d. centered with variance $1/n$ for
$k,l\in[n]$ such that $k\leq l$ and $W_{kl}=\overline{W_{lk}}$ for $k,l\in[n]$
such that $l<k$. It is well-known [1, 12, 24] that independent Wigner matrices
$W_{1},\ldots,W_{d}$ satisfy Assumptions (A.1), (A.2), (A.3) (by Lemma 3.3)
and (A.4). Their joint law $(W_{i})_{i\in[d]}$ converges weakly in probability
to a free semicircular family $(s_{i})_{i\in[d]}$. Therefore, Theorem 2.2
holds and the limit in Theorem (2.2.i) is expressed as

$\displaystyle\Delta\underset{n\to\infty}{\Rightarrow}\frac{1}{\sqrt{d}}\sum_{i\in[d]}s_{i}\otimes
s_{i},$

in probability and in expectation.

Let us interpret Example 4.3 above in terms of the corresponding random
quantum channel. Given $W\in M_{n}(\mathbb{C})$ a Wigner matrix, we have
$\operatorname{\mathbb{E}}(W_{kl})=0$ and
$\operatorname{\mathbb{E}}|W_{kl}|^{2}=1/n$ for all $k,l\in[n]$ and
$(W_{kl})_{k\leq l\in[n]}$ independent. Therefore,

$\operatorname{\mathbb{E}}\left(W^{2}\right)=\sum_{k\in[n]}E_{kk}=\operatorname{Id}.$

So we can conclude that, sampling $W_{1},\ldots,W_{d}\in M_{n}(\mathbb{C})$
independent Wigner matrices and setting $K_{i}=W_{i}/\sqrt{d}$ for each
$i\in[d]$, the random completely positive map

$\Phi:X\in M_{n}(\mathbb{C})\mapsto\sum_{i\in[d]}K_{i}XK_{i}\in
M_{n}(\mathbb{C})$

is, on average, trace-preserving and unital. Moreover, the spectral
distribution of $d(\Phi-\operatorname{\mathbb{E}}(\Phi))$ converges weakly
almost surely to $\mu_{sc}^{\circledast d}$ as $n\to\infty$, and the spectral
distribution of $\sqrt{d}(\Phi-\operatorname{\mathbb{E}}(\Phi))$ converges
weakly in probability and in expectation to $\mu_{sc}$ as $n,d\to\infty$. In
this case, we can actually explicitly compute
$\operatorname{\mathbb{E}}(\Phi)$. Indeed, denoting by
$\\{e_{1},\ldots,e_{n}\\}$ the canonical basis of $\mathbb{C}^{n}$ and setting
$E_{kl}=e_{k}e_{l}^{*}$ for each $k,l\in[n]$, we have

$\operatorname{\mathbb{E}}\left(W\otimes\overline{W}\right)=\frac{1}{n}\sum_{k,l\in[n]}E_{kl}\otimes
E_{kl}+\frac{1}{n}\sum_{k\neq l\in[n]}E_{kl}\otimes
E_{lk}=\psi\psi^{*}+\frac{1}{n}\left(F-\operatorname{diag}(F)\right),$

where $\psi=\sum_{k\in[n]}(e_{k}\otimes e_{k})/\sqrt{n}$ is a maximally
entangled unit vector, $F=\sum_{k,l\in[n]}E_{kl}\otimes E_{lk}$ is the flip
operator and $\operatorname{diag}(F)=\sum_{k\in[n]}E_{kk}\otimes E_{kk}$ is
its diagonal part (with respect to the canonical product basis of
$\mathbb{C}^{n}\otimes\mathbb{C}^{n}$). $F-\operatorname{diag}(F)$ has
spectrum $\\{1,-1,0\\}$, where $1$ and $-1$ have multiplicities $n(n-1)/2$
(with associated eigenvectors $\\{(e_{k}\otimes e_{l}+e_{l}\otimes
e_{k})/\sqrt{2},\ k<l\in[n]\\}$ and $\\{(e_{k}\otimes e_{l}-e_{l}\otimes
e_{k})/\sqrt{2},\ k<l\in[n]\\}$ respectively) and $0$ has multiplicity $n$.
Hence,

$\operatorname{\mathbb{E}}(M_{\Phi})=\psi\psi^{*}+\frac{1}{n}\left(F-\operatorname{diag}(F)\right).$

This can be re-written at the level of $\Phi$ as

$\operatorname{\mathbb{E}}(\Phi):X\in
M_{n}(\mathbb{C})\mapsto\operatorname{tr}(X)\frac{\operatorname{Id}}{n}+\frac{1}{n}\left(X^{T}-\operatorname{diag}(X)\right)\in
M_{n}(\mathbb{C}),$

where $X^{T}$ denotes the transposition of $X$ and $\operatorname{diag}(X)$
its diagonal part, both with respect to the canonical basis of
$\mathbb{C}^{n}$. This means that, up to a correction that vanishes as $n$
grows, $\operatorname{\mathbb{E}}(\Phi)$ is the so-called fully randomizing
channel $\Pi:X\mapsto\operatorname{tr}(X)\operatorname{Id}/n$, which has one
eigenvalue equal to $1$, with associated eigenvector the maximally mixed state
$\operatorname{Id}/n$, and all the other eigenvalues equal to $0$. Our result
thus gives a precise understanding of how the asymptotic spectrum of $\Phi$
deviates from the flat one of $\Pi$, for $d$ either fixed or growing.

We can also consider non-homogeneous matrices as in [19, 5].

######  Example 4.4 (Non-homogeneous matrices).

Let $W_{1},\ldots,W_{d}\in M_{n}(\mathbb{C})$ be independent Gaussian matrices
such that, for each $i\in[d]$,

$\displaystyle\|\operatorname{\mathbb{E}}(W_{i})\|\underset{n\to\infty}{\to}0;$
$\displaystyle\|\operatorname{\mathbb{E}}(W_{i}^{2})-\text{Id}\|\underset{n\to\infty}{\to}0.$

Then it was proven in [5, Theorem 2.10] that $(W_{i})_{i\in[d]}$ are
asymptotically free in expectation and almost surely as long as the covariance
structure of each one goes to $0$ sufficiently fast. Concretely, we have to
impose that, for each $i\in[d]$,

$\displaystyle
v(W_{i}):=\|\operatorname{cov}(W_{i})\|=o\left(\log^{-3/2}n\right).$

The limit is a free semicircular family $(s_{i})_{i\in[d]}$. So Theorem
(2.2.i) implies that

$\displaystyle\Delta\underset{n\to\infty}{\Rightarrow}\frac{1}{\sqrt{d}}\sum_{i\in[d]}s_{i}\otimes
s_{i},$

in probability and in expectation. As a matter of fact, similar results hold
for Hermitian random matrices with bounded entries (see [7, Theorem 3.25]).

Another classical example is Wishart-type matrices.

######  Example 4.5 (Wishart matrices).

A random Hermitian matrix $W\in M_{n}(\mathbb{C})$ is a (centered) Wishart
matrix if

$\displaystyle
W=XX^{*}-\operatorname{\mathbb{E}}(XX^{*})=XX^{*}-\operatorname{Id},$

where $X\in M_{n}(\mathbb{C})$ has i.i.d. centered entries with variance
$1/n$. Independent Whishart matrices $W_{1},\ldots,W_{d}$ also satisfy
Assumptions (A.1), (A.2), (A.3) (by Lemma 3.3) and (A.4); see [9]. The limit
of a Wishart matrix $XX^{*}$ is equal to $cc^{*}$, where $c$ is a circular
element, i.e.

$\displaystyle c=\frac{s+is^{\prime}}{\sqrt{2}},$

with $s,s^{\prime}$ free semicircular variables, and it is called the quarter
circular law or Marchenko-Pastur law. The limit law in (2.2.i) is given by

$\displaystyle\Delta\underset{n\to\infty}{\Rightarrow}\frac{1}{\sqrt{d}}\sum_{j\in[d]}(c_{j}c_{j}^{*}-1)\otimes(c_{j}c_{j}^{*}-1),$

in probability and in expectation, where $(c_{j})_{j\in[d]}$ are free copies
of $c$.

##  5\. Perspectives

The results presented so far are general, and several natural classes of
random matrices satisfy Assumptions (A.1),(A.2),(A.3), and (A.4). A natural
extension would be to remove assumption (A.3). It turns out, however, that
such an assumption is necessary in order for the limit in Theorem (2.2.ii) to
be semicircular. To understand this problem, recall that the tensor
convolution $\mu_{1}\circledast\mu_{2}$ is not the free convolution of
$\mu_{1}\star\mu_{1}$ with $\mu_{2}\star\mu_{2}$. This is because if $a$ and
$b$ are free variables, then their tensor products $a\otimes a$ and $b\otimes
b$ are not necessarily free. This was proven in [10]. In particular, it is
surprising that the semicircular law can still be recovered for tensor
products under only the mean zero assumption. However, a deep understanding of
the structure of tensor products is still lacking, and while classical
freeness is lost, tensor products still enjoy “partial” freeness properties.
In a work in progress, we aim to decrypt such a structure in the framework of
Free Probability.

A second extension of our results would be the non-Hermitian case, addressing,
in particular, the conjecture mentioned in the introduction. Indeed, the
methods presented here can be easily generalized to capture the convergence of
the $*$-moments of non-Hermitian matrices $W_{1},\ldots,W_{d}\in
M_{n}(\mathbb{C})$, that is, moments of the form

$\displaystyle\operatorname{\mathbb{E}}\left(\tau^{(n)}(W_{i_{1}}^{w_{1}}\cdots
W_{i_{k}}^{w_{k}})\right),$

where $i_{1},\ldots,i_{k}\in[d]$ and $w_{1},\ldots,w_{k}\in\\{1,*\\}$.
However, the spectral distribution of non-Hermitian matrices is no longer
characterized by the $*$-moments convergence or singular value convergence. A
second ingredient, which is in general more involved, is to lower bound the
smallest singular value of $W_{1},\ldots,W_{d}$ uniformly; see [6] and
references within. The first setting we might want to study is when the Kraus
operators of the random quantum channel are sampled as independent Haar
unitaries. In future work, we aim to establish the quantitative invertibility
of the model by deriving lower bounds for the corresponding smallest singular
value.

## References

  * Anderson et al. [2010] G. W. Anderson, A. Guionnet, and O. Zeitouni.  _An introduction to random matrices_.  Number 118 in Cambridge Studies in Advanced Mathematics. Cambridge University Press, 2010.
  * Aubrun and Nechita [2012] G. Aubrun and I. Nechita.  Realigning random states.  _J. Math. Phys._ , 53(10):102210, 2012.
  * Aubrun and Szarek [2017] G. Aubrun and S. Szarek.  _Alice and Bob meet Banach: The interface of asymptotic geometric analysis and quantum information theory_ , volume 223 of _Mathematical Surveys and Monographs_.  American Mathematical Society, 2017.
  * Bai and Silverstein [2010] Z. Bai and J. W. Silverstein.  _Spectral analysis of large dimensional random matrices_ , volume 20 of _Springer Series in Statistics_.  Springer, 2010.
  * Bandeira et al. [2023] A. S. Bandeira, M. T. Boedihardjo, and R. van Handel.  Matrix concentration inequalities and free probability.  _Inventiones mathematicae_ , pages 1–69, 2023.
  * Bordenave and Chafaï [2012] C. Bordenave and D. Chafaï.  Around the circular law.  2012\.
  * Brailovskaya and van Handel [2022] T. Brailovskaya and R. van Handel.  Universality and sharp matrix concentration inequalities.  _arXiv preprint arXiv:2201.05142_ , 2022.
  * Bruzda et al. [2009] W. Bruzda, V. Cappellini, H.-J. Sommers, and K. Życzkowski.  Random quantum operations.  _Phys. Lett. A_ , 373(3):320–324, 2009.
  * Capitaine and Donati-Martin [2007] M. Capitaine and C. Donati-Martin.  Strong asymptotic freeness for Wigner and Wishart matrices.  _Indiana University mathematics journal_ , 56(2):767–803, 2007.
  * Collins and Lamarre [2017] B. Collins and P. Y. G. Lamarre.  *-freeness in finite tensor products.  _Advances in Applied Mathematics_ , 83:47–80, 2017.
  * Durrett [2019] R. Durrett.  _Probability: theory and examples_ , volume 49.  Cambridge university press, 2019.
  * Dykema [1993] K. Dykema.  On certain free product factors via an extended matrix model.  _Journal of functional analysis_ , 112(1):31–60, 1993.
  * Girko [1983] V. L. Girko.  The circle law.  _Teor. Veroyatnost. i Mat. Statist._ , (28):15–21, 1983.  ISSN 0131-6982.
  * Girko [1994] V. L. Girko.  The circular law: ten years later.  _Random Oper. Stochastic Equations_ , 2(3):235–276, 1994.  ISSN 0926-6364,1569-397X.
  * González-Guillén et al. [2018] C. E. González-Guillén, M. Junge, and I. Nechita.  On the spectral gap of random quantum channels.  _arXiv:1811.08847_ , 2018.
  * Hastings [2007] M. B. Hastings.  Random unitaries give quantum expanders.  _Phys. Rev. A_ , 76:032315, 2007.  doi: 10.1103/PhysRevA.76.032315.
  * Kesten [1959] H. Kesten.  Symmetric random walks on groups.  _Transactions of the American Mathematical Society_ , 92(2):336–354, 1959.
  * Lancien and Pérez-García [2022] C. Lancien and D. Pérez-García.  Correlation length in random MPS and PEPS.  _Ann. Henri Poincaré_ , 23:141–222, 2022.  doi: 10.1007/s00023-021-01087-4.
  * Lancien and Youssef [2023] C. Lancien and P. Youssef.  A note on quantum expanders.  _arXiv preprint arXiv:2302.07772_ , 2023.
  * McKay [1981] B. D. McKay.  The expected eigenvalue distribution of a large regular graph.  _Linear Algebra and its applications_ , 40:203–216, 1981\.
  * Nica and Speicher [2006] A. Nica and R. Speicher.  _Lectures on the combinatorics of free probability_ , volume 13.  Cambridge University Press, 2006.
  * Pisier [2014] G. Pisier.  Quantum expanders and geometry of operator spaces.  _J. Eur. Math. Soc._ , 16(6):1183–1219, 2014\.  doi: 10.4171/JEMS/458.
  * Tao and Vu [2010] T. Tao and V. Vu.  Random matrices: universality of ESDs and the circular law.  _Ann. Probab._ , 38(5):2023–2065, 2010.  ISSN 0091-1798,2168-894X.  With an appendix by Manjunath Krishnapur.
  * Voiculescu [1991] D. Voiculescu.  Limit laws for random matrices and free products.  _Inventiones mathematicae_ , 104(1):201–220, 1991\.
  * Wolf [2012] M. M. Wolf.  Quantum chanels & operations: Guided tour.  2012\.  URL https://mediatum.ub.tum.de/doc/1701036/1701036.pdf.


# Velocity continuation with Fourier neural operators for accelerated
uncertainty quantification

Ali Siahkoohi, Mathias Louboutin, and Felix J. Herrmann
Georgia Institute of Technology
{alisk, mlouboutin3<EMAIL_ADDRESS>
###### Abstract

Seismic imaging is an ill-posed inverse problem that is challenged by noisy
data and modeling inaccuracies—due to errors in the background squared-
slowness model. Uncertainty quantification is essential for determining how
variability in the background models affects seismic imaging. Due to the costs
associated with the forward Born modeling operator as well as the high
dimensionality of seismic images, quantification of uncertainty is
computationally expensive. As such, the main contribution of this work is a
survey-specific Fourier neural operator surrogate to velocity continuation
that maps seismic images associated with one background model to another
virtually for free. While being trained with only 200 background and seismic
image pairs, this surrogate is able to accurately predict seismic images
associated with new background models, thus accelerating seismic imaging
uncertainty quantification. We support our method with a realistic data
example in which we quantify seismic imaging uncertainties using a Fourier
neural operator surrogate, illustrating how variations in background models
affect the position of reflectors in a seismic image.

##  1 Introduction

Seismic imaging involves estimating the short-wavelength component of the
Earth’s subsurface squared-slowness model—known as the seismic image—given
shot records and an estimation of the smooth background squared-slowness
model. This linearized imaging problem is challenged by the computationally
expensive forward operator as well as presence of measurements noise,
linearization errors, modeling errors, and the nontrivial nullspace of the
linearized forward Born modeling operator [1, 2, 3]. These challenges
highlight the importance of uncertainty quantification (UQ) in seismic
imaging, where instead of finding one seismic image estimate, a distribution
of seismic images is obtained that explains the observed data [4],
consequently reducing the risk of data overfit and enabling UQ [5, 6, 7, 8, 9,
10, 11, 12, 13, 14, 15].

The seismic imaging uncertainty can be attributed to two main sources [16, 17,
18]: (1) errors in the data, which include measurement and linearization
errors; and (2) modeling errors, which include errors in the estimation of the
background squared-slowness model. In this paper, we focus on uncertainties
with respect to the background model as it has the main contributing factor to
imaging uncertainty due to its effect on reflector positioning [19, 20, 18].
To this end, we assume there are no measurement or linearization errors and
consider the map from shot records to the seismic image to be represented by
the deterministic reverse-time migration (RTM) algorithm. In this setup,
quantifying the uncertainty in seismic imaging—due to errors in the background
model—involves computing numerous RTMs with all the background model posterior
samples [8, 9]. Due to the high-dimensionality of seismic images, the number
of seismic imaging posterior samples required to obtain accurate estimations
of the imaging posterior moments is large, making UQ with respect to the
background model computationally expensive for 2D and unfeasible for large 3D
problems. These costs can be alleviated via velocity continuation methods [21,
22, 23, 24, 25] that map seismic images associated with one background model
to another, without directly solving the imaging problem. Velocity
continuation methods are designed to be frugal compared to computing a new RTM
for each new background model, which makes them suitable for large-scale
seismic imaging uncertainty quantification [19, 20].

While several velocity continuation methods have been proposed [21, 22, 23,
24, 25], they typically involve solving partial differential equations (PDEs),
which can be still costly in the context of seismic imaging UQ. To address
this challenge, we propose a neural network surrogate for velocity
continuation that is capable of mapping seismic images associated with one
background model to another with negligible computational cost. Motivated by
the success of Fourier neural operators [FNOs, 26] in approximating the
solution operator of PDEs [27, 28, 29], we chose them as the architecture for
our neural network surrogate. Due to our main interest in accelerating
velocity continuation in the context of UQ, we train a survey-specific FNO
that acts as a surrogate for velocity continuation for the specific survey at
hand. This choice, while not offering generalizations across different survey
areas in the Earth, can speed up seismic imaging UQ for the survey at hand,
which involves computing seismic images associated with many background model
posterior samples. We show that the FNO can be trained using $200$ pairs of
background and seismic image pairs, making the survey-specific training
procedure computationally viable. To scale this method to industry size
problems, transfer learning [30] can further reduce the upfront costs of
training the FNO. After training, the FNO can be used to obtain samples from
the imaging posterior almost free of cost.

In the next section, we describe seismic imaging by introducing the forward
Born modeling operator, through which the seismic image relates to the
background model. Next, we define velocity continuation, followed by
describing our proposed FNO-based approach for accelerating seismic imaging
UQ. Finally, we evaluate the performance of the trained neural operator on a
realistic dataset, and we demonstrate how the imaging uncertainty affects the
positioning of the reflectors in the seismic image.

##  2 Theory

We introduce a deep-network surrogate for velocity continuation in order to
enable faster quantification of uncertainty—due to errors in the background
velocity model—in seismic imaging. We begin with an introduction to seismic
imaging and the linearized forward model associated with it.

###  2.1 Seismic imaging

The inverse problem that we tackle involves the process of estimating the
short-wavelength component of the Earth’s unknown subsurface squared-slowness
model given measurements recorded at the surface. This problem, also known as
seismic imaging, can be formulated as a linear inverse problem by linearizing
the nonlinear relationship between shot records and the squared-slowness
model, governed by the wave-equation. In its simplest acoustic form, the
linearization with respect to the slowness model—around a background squared
slowness model $\mathbf{m}_{0}$—leads to a linear inverse problem for
estimating the ground truth seismic image $\delta\mathbf{m}^{\ast}$ with the
following forward model,

$\delta\mathbf{d}_{i}=\mathbf{J}(\mathbf{m}_{0},\mathbf{q}_{i})\delta\mathbf{m}^{\ast}.$
(1)

In the above expression,
$\delta\mathbf{d}=\left\\{\delta\mathbf{d}_{i}\right\\}_{i=1}^{n_{s}}$ are
$n_{s}$ linearized shot records and
$\mathbf{J}(\mathbf{m}_{0},\mathbf{q}_{i})$ represents the linearized Born
scattering operator. This operator is parameterized by the source signature
$\mathbf{q}_{i}$ and the background squared-slowness model $\mathbf{m}_{0}$,
which is typically estimated in the previous inversion steps [8, 9]. In this
work, we focus on potential inaccuracies in the background model, the main
source of variability and uncertainty in seismic imaging [19, 20, 18].
Therefore, we assume there are not measurement and linearization errors, which
leads to a deterministic mapping from
$\left\\{\delta\mathbf{d}_{i}\right\\}_{i=1}^{n_{s}}$ to $\delta\mathbf{m}$
for a given background model. We use RTM for this deterministic map, which
involves applying the adjoint linearized Born scattering operator to the
linearized shot records for all source experiments,

$\delta\mathbf{m}_{\text{RTM}}=\sum_{i=1}^{n_{s}}\mathbf{J}(\mathbf{m}_{0},\mathbf{q}_{i})^{\top}\delta\mathbf{d}_{i}.$
(2)

Since this map is deterministic, the uncertainty in the background model can
be translated into uncertainty in seismic imaging by evaluating Equation 2 for
the background models sampled from the posterior distribution
$p(\mathbf{m}_{0}\mid\mathbf{d})$, where $\mathbf{d}$ represents the shot
records with no linearization (the field data) [8, 9]. Bayesian inference
involving high-dimensional posterior distributions, for example seismic
imaging and full-waveform inversion, requires many samples from the posterior
distribution for accurate Monte Carlo approximation of high-dimensional
integrals [31]. As a result, the mapping in Equation 2 must be evaluated over
numerous background model samples from $p(\mathbf{m}_{0}\mid\mathbf{d})$,
which is computationally expensive due to costs of applying
$\mathbf{J}(\mathbf{m}_{0},\mathbf{q}_{i})^{\top},\,i=1,\ldots,n_{s}$. To
address this computational challenge, our proposed method involves a FNO-based
velocity continuation approach that is capable of mapping seismic images
associated with one background model to another, effectively replacing a
costly demingration-migration. We describe this approach in the next section.

###  2.2 Velocity continuation

At its core, velocity continuation is a process which alters a seismic image
based on changes in the background model [21, 22, 24, 25]. Using this
technique, an initial seismic image associated with an initial background
model $\mathbf{m}_{\text{init}}$ is altered to approximate the seismic image
associated with a target background model $\mathbf{m}_{\text{target}}$ without
computing Equation 2. This process can be interpreted as a map [23]:

$\mathcal{T}_{\left(\mathbf{m}_{\text{init}},\mathbf{m}_{\text{target}}\right)}:\delta\mathcal{M}\to\delta\mathcal{M},$
(3)

where $\delta\mathcal{M}$ denotes the space of seismic images, and the
velocity continuation map is parameterized by $\mathbf{m}_{\text{init}}$ and
$\mathbf{m}_{\text{target}}$. We take advantage of recent advances in deep
learning to train a surrogate model for velocity continuation that
approximates
$\mathcal{T}_{\left(\mathbf{m}_{\text{init}},\mathbf{m}_{\text{target}}\right)}=\mathbf{J}(\mathbf{m}_{\text{target}},\mathbf{q}_{i})^{\top}(\mathbf{J}(\mathbf{m}_{\text{init}},\mathbf{q}_{i})^{\top})^{\dagger}$,
which can be evaluated virtually for free instead of solving four PDEs.
Through this approach, we are able to significantly accelerate velocity
continuation, trading wave-equation solves for a simple neural network
inference. This is of considerable importance in the context of seismic
imaging UQ. Before demonstrating the benefits of our approach, we introduce
FNOs in the context of velocity continuation.

###  2.3 Fourier neural operators for velocity continuation

In light of their success in learning mesh-free solution operators to PDEs
[27, 28, 29], we choose FNOs [26] as a surrogate model for velocity
continuation—a process that can be interpreted as a double linearized PDE
solve. The main components of FNOs are the Fourier layers, which involve a
Fourier transform over the spatial dimensions of their input, followed by a
learned pointwise multiplication and an inverse Fourier transform. These
layers act as long-kernel convolutional layer, akin to pseudo-spectral
methods, which explains the representation power of FNOs [29]. To train a FNO
as a surrogate model for velocity continuation, we define it as a map

$\mathcal{G}_{\mathbf{w}}:\mathcal{M}\times\delta\mathcal{M}\to\delta\mathcal{M},$
(4)

where $\mathcal{G}_{\mathbf{w}}$ denotes the FNO with weights $\mathbf{w}$,
and $\mathcal{M}$ is the space of background models. In words, we design
$\mathcal{G}_{\mathbf{w}}$ as a neural network that takes as input the target
background model and the initial seismic image, and outputs the target seismic
image. This choice is analogous to the structure of the velocity continuation
map
$\mathcal{T}_{\left(\mathbf{m}_{\text{init}},\mathbf{m}_{\text{target}}\right)}$
(cf. Equation 3), except that we fix the initial background model
$\mathbf{m}_{\text{init}}$ to an arbitrary background model posterior sample,
hence making the dependence of $\mathcal{G}_{\mathbf{w}}$ to
$\mathbf{m}_{\text{init}}$ implicit. With this choice of inputs and output for
$\mathcal{G}_{\mathbf{w}}$, the FNO’s task is to perturb the input initial
seismic image according to the provided target background model in order to
predict the target seismic image.

Due to our interest in accelerating velocity continuation in the context of
UQ, we train a survey-specific FNO that acts as a surrogate for velocity
continuation for the specific survey at hand. This choice, while not offering
generalizations across different survey areas in the Earth, can speed up
seismic imaging UQ for the survey at hand, which involves computing seismic
images associated with many posterior background model samples
$\mathbf{m}_{0}\sim p(\mathbf{m}_{0}\mid\mathbf{d})$. By training the FNO on a
small number of these background model and seismic image pairs, i.e.,
approximately $200$, we can accelerate the velocity continuation process for
the rest of the background model posterior samples while limiting the risk of
introducing generalization errors due to the strong heterogeneity of Earth. To
achieve this, we construct a set of $N$ training input-output pairs in the
form of

$\bigg{\\{}\Big{(}\big{(}\mathbf{m}_{0}^{(i)},\delta\mathbf{m}_{\text{init}}\big{)},\delta\mathbf{m}_{\text{RTM}}^{(i)}\Big{)}\
\Big{|}\ i=1,\ldots,N\bigg{\\}},$ (5)

where $(\mathbf{m}_{0}^{(i)},\delta\mathbf{m}_{\text{init}})$ is the input
target background and initial seismic image training pair, and
$\delta\mathbf{m}_{\text{RTM}}^{(i)}$ is the associated target seismic image.
Training involves minimizing the squared $\ell_{2}$-norm of the difference
between the FNO output and the target seismic image with respect to FNO
weights,

$\mathbf{w}^{\ast}=\operatorname*{arg\,min}_{\mathbf{w}}\,\frac{1}{N}\sum_{i=1}\|\mathcal{G}_{\mathbf{w}}(\mathbf{m}_{0}^{(i)},\delta\mathbf{m}_{\text{init}})-\delta\mathbf{m}_{\text{RTM}}^{(i)}\|_{2}^{2}.$
(6)

We solve optimization problem 6 with the Adam stochastic optimization
algorithm [32]. After training, the trained FNO approximates the velocity
continuation map for a fixed homogeneous initial background model, i.e.,

$\mathcal{G}_{\mathbf{w}^{\ast}}(\mathbf{m}_{\text{target}},\delta\mathbf{m}_{\text{init}})\approx\mathcal{T}_{\left(\mathbf{m}_{\text{init}},\mathbf{m}_{\text{target}}\right)}(\delta\mathbf{m}_{\text{init}})=\delta\mathbf{m}_{\text{target}},$
(7)

where $\mathbf{m}_{\text{init}}$ is fixed to an arbitrary background model
posterior sample and $\delta\mathbf{m}_{\text{target}}$ denotes the seismic
image associated with $\delta\mathbf{m}_{\text{init}}$. Given background model
samples $\mathbf{m}_{0}\sim p(\mathbf{m}_{0}\mid\mathbf{d})$ as input, the FNO
outputs samples from the imaging posterior,
$p(\mathbf{m}_{0}\mid\delta\mathbf{d})$. This process accelerates seismic
imaging uncertainty quantification as no further RTMs (Equation 2) need to be
computed. Training the FNO with as little as $200$ training pairs makes the
survey-specific training procedure computationally viable. To further
accelerate the process, training of the FNO can be started as the same time as
the background model posterior sampling phase, using the already collected
posterior samples as training data. In the next section, we show the results
of approximating the velocity continuation map with a FNO via a quasi-real
seismic experiment.

##  3 Numerical experiments

The purpose of the presented numerical experiments here is to demonstrate the
ability to approximate the velocity continuation map (Equation 3) using a FNO.
We further show how this trained FNO can be used for UQ by showing the effect
of imaging uncertainty on the positioning of the reflectors. We begin by
describing the training setup, including the seismic acquisition geometry.

###  3.1 Acquisition geometry and training configuration

Our examples involve imaging a 2D subset of the Parihaka [33, 34] prestrack
Kirchhoff migration field dataset. We use this 2D section to create linearized
data according to the linear forward model in Equation 1, where we consider no
measurement and linearization errors. The model is discretized on a
$12.5\,\mathrm{m}$ vertical and $20\,\mathrm{m}$ horizontal grid, and the data
is acquired with $102$ equally spaced sources and $204$ fixed receivers. We
use a Ricker wavelet with a central frequency of $30\,\mathrm{Hz}$ as the
source signature. For presentation purposes, we augment a $125\,\mathrm{m}$
water column on top of these models to limit the near source imaging
artifacts. We create background models by assuming access to an oracle, which
provides geologically consistent background models with the 2D section. Given
$200$ background models, we migrate the simulated seismic data via Equation 2
to obtain the corresponding seismic images. Figures 1a and 1b show the
vertical profile of five randomly selected background and seismic images at
$2.50\,\mathrm{km}$, respectively. These images, displaying strong amplitude
and phase differences, highlight the importance of quantifying the uncertainty
in imaging when dealing with errors in the background model. To reduce the
costs associated with UQ, we train the FNO surrogate using the $200$ training
pairs (cf. Equation 5) for $500$ epochs with the Adam optimizer [32].

(a)

(b)

Figure 1: Variation among five (a) background; and (b) seismic images (RTMs),
plotted as a vertical profile at $2.50\,\mathrm{km}$.

###  3.2 Results

To evaluate the accuracy of the trained FNO in predicting seismic images, we
compare its output to the target seismic image, i.e., the RTM image obtained
via Equation 2 when using the target background model (Figure 2a) to
parameterize the Born scattering forward operator. This is conducted over the
testing dataset, which is derived using the same procedure as the training
dataset using the background model creating oracle. Figures 2b and 2c show the
target and predicted seismic images, respectively, and Figures 2d includes the
difference between them. We observe that the network has accurately predicted
the target image, where errors are mostly due to amplitude differences. The
accuracy of the prediction in terms of phase and amplitude can be further
confirmed by focusing on two vertical profiles, at $2.50\,\mathrm{km}$ (Figure
2e) and $4.125\,\mathrm{km}$ (Figure 2f) horizontal locations, which include
regions with fault and torturous reflectors. We use the trained FNO for UQ by
providing testing background models and predicting the associated seismic
images. We visualize the obtained uncertainties by showing the variability in
the location of reflectors, determined via an automatic horizon tracking
software [35]. We pass the seismic images predicted by the FNO to the horizon
tracker for $25$ selected horizons. As a result, we obtain multiple instances
of each horizons, from which we compute pointwise mean and standard
deviations. Figure 3 indicates the result where the solid lines correspond to
the mean among different instances of each horizons and shaded areas indicate
the mean plus and minus the pointwise standard deviation. These shaded areas
indicate uncertainties in the location of the reflectors, which are due to the
variability in the background models. As expected, we find a general increase
of uncertainty with depth. We also observe that the areas of high uncertainty
are correlated with areas of poor illumination, faults and tortuous
reflectors.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 2: Velocity continuation with FNOs. Target (a) background; and (b)
seismic images. (c) Predicted seismic image with the FNO. (d) Difference
between target and prediction. (e) and (f) Vertical profile comparisons
between target (dashed blue) and predicted (red) seismic images. Figure 3:
Uncertainty in the tracked horizons.

In this example, we use Devito [36, 37] for the wave-equation based
simulations. We based our PyTorch FNO implementation on the original
implementation. The code to reproduce our results are made available on
GitHub.

##  4 Conclusions and discussion

Quantifying the uncertainty in seismic imaging due to errors in the background
model involves solving many seismic imaging problems that vary in the
parameterization of the background model of the forward operator. To reduce
the computational cost of this process—mainly due to the computational costs
of the forward operator—we proposed to train a survey-specific Fourier neural
operator surrogate that mimics velocity continuation. This surrogate model
maps seismic images associated with one background model to another virtually
for free, which has the benefit of accelerating uncertainty quantification. We
showed that this surrogate model can be trained with as few as $200$ training
pairs while still providing a good seismic image prediction accuracy. Further
research is required in training a reliable global surrogate, being able to
generalize across other survey areas and more realistic physics.

##  5 Acknowledgment

This research was carried out with the support of Georgia Research Alliance
and partners of the ML4Seismic Center.

## References

  * Lambaré et al. [1992] Gilles Lambaré, Jean Virieux, Raul Madariaga, and Side Jin.  Iterative asymptotic inversion in the acoustic approximation.  _Geophysics_ , 57(9):1138–1154, 1992.
  * Schuster [1993] Gerard T Schuster.  Least-squares cross-well migration.  In _SEG Technical Program Expanded Abstracts 1993_ , pages 110–113. Society of Exploration Geophysicists, 1993.
  * Nemeth et al. [1999] Tamas Nemeth, Chengjun Wu, and Gerard T Schuster.  Least-squares migration of incomplete reflection data.  _GEOPHYSICS_ , 64(1):208–221, 1999.  doi: 10.1190/1.1444517.
  * Tarantola [2005] Albert Tarantola.  _Inverse problem theory and methods for model parameter estimation_.  SIAM, 2005.  ISBN 978-0-89871-572-9.  doi: 10.1137/1.9780898717921.
  * Malinverno and Parker [2006] Alberto Malinverno and Robert L Parker.  Two ways to quantify uncertainty in geophysical inverse problems.  _GEOPHYSICS_ , 71(3):W15–W27, 2006.  doi: 10.1190/1.2194516.
  * Martin et al. [2012] James Martin, Lucas C. Wilcox, Carsten Burstedde, and OMAR Ghattas.  A Stochastic Newton MCMC Method for Large-scale Statistical Inverse Problems with Application to Seismic Inversion.  _SIAM Journal on Scientific Computing_ , 34(3):A1460–A1487, 2012.  URL http://epubs.siam.org/doi/abs/10.1137/110845598.
  * Ray et al. [2017] Anandaroop Ray, Sam Kaplan, John Washbourne, and Uwe Albertin.  Low frequency full waveform seismic inversion within a tree based Bayesian framework.  _Geophysical Journal International_ , 212(1):522–542, 10 2017.  ISSN 0956-540X.  doi: 10.1093/gji/ggx428.  URL https://doi.org/10.1093/gji/ggx428.
  * Zhu et al. [2016] Hejun Zhu, Siwei Li, Sergey Fomel, Georg Stadler, and Omar Ghattas.  A bayesian approach to estimate uncertainty for full-waveform inversion using a priori information from depth migration.  _Geophysics_ , 81(5):R307–R323, 2016.
  * Fang et al. [2018] Zhilong Fang, Curt Da Silva, Rachel Kuske, and Felix J. Herrmann.  Uncertainty quantification for inverse problems with weak partial-differential-equation constraints.  _GEOPHYSICS_ , 83(6):R629–R647, 2018.  doi: 10.1190/geo2017-0824.1.
  * Stuart et al. [2019] Georgia K Stuart, Susan E Minkoff, and Felipe Pereira.  A two-stage Markov chain Monte Carlo method for seismic inversion and uncertainty quantification.  _GEOPHYSICS_ , 84(6):R1003–R1020, 11 2019.  doi: 10.1190/geo2018-0893.1.
  * Herrmann et al. [2019] Felix J. Herrmann, Ali Siahkoohi, and Gabrio Rizzuti.  Learned imaging with constraints and uncertainty quantification.  In _Neural Information Processing Systems (NeurIPS) 2019 Deep Inverse Workshop_ , 12 2019.  URL https://arxiv.org/pdf/1909.06473.pdf.
  * Siahkoohi et al. [2020a] Ali Siahkoohi, Gabrio Rizzuti, and Felix J. Herrmann.  A deep-learning based bayesian approach to seismic imaging and uncertainty quantification.  _82nd EAGE Conference and Exhibition 2020_ , 6 2020a.  URL https://slim.gatech.edu/Publications/Public/Submitted/2020/siahkoohi2020EAGEdlb/siahkoohi2020EAGEdlb.html.
  * Siahkoohi et al. [2020b] Ali Siahkoohi, Gabrio Rizzuti, and Felix J. Herrmann.  Uncertainty quantification in imaging and automatic horizon tracking—a Bayesian deep-prior based approach.  In _SEG Technical Program Expanded Abstracts 2020_ , pages 1636–1640, 9 2020b.  doi: 10.1190/segam2020-3417560.1.
  * Siahkoohi and Herrmann [2021] Ali Siahkoohi and Felix J Herrmann.  Learning by example: fast reliability-aware seismic imaging with normalizing flows.  In _First International Meeting for Applied Geoscience & Energy_, pages 1580–1585. Society of Exploration Geophysicists, 2021.
  * Siahkoohi et al. [2021] Ali Siahkoohi, Gabrio Rizzuti, and Felix J Herrmann.  Deep bayesian inference for seismic imaging with tasks.  _arXiv preprint arXiv:2110.04825_ , 2021.
  * Thore et al. [2002] Pierre Thore, Arben Shtuka, Magali Lecour, Taoufik Ait-Ettajer, and Richard Cognot.  Structural uncertainties: Determination, management, and applications .  _Geophysics_ , 67(3):840–852, 2002.
  * Osypov et al. [2013] Konstantin Osypov, Yi Yang, Aimé Fournier, Natalia Ivanova, Ran Bachrach, Can Evren Yarman, Yu You, Dave Nichols, and Marta Woodward.  Model-uncertainty quantification in seismic tomography: method and applications.  _Geophysical Prospecting_ , 61(6-Challenges of Seismic Imaging and Inversion Devoted to Goldin):1114–1134, 2013.
  * Ely et al. [2018] Gregory Ely, Alison Malcolm, and Oleg V. Poliannikov.  Assessing uncertainties in velocity models and images with a fast nonlinear uncertainty quantification method.  _GEOPHYSICS_ , 83(2):R63–R75, 2018.  doi: 10.1190/geo2017-0321.1.  URL https://doi.org/10.1190/geo2017-0321.1.
  * Fomel and Landa [2014] Sergey Fomel and Evgeny Landa.  Structural uncertainty of time-migrated seismic images.  _Journal of Applied Geophysics_ , 101:27–30, 2014.
  * Poliannikov and Malcolm [2016] Oleg V Poliannikov and Alison E Malcolm.  The effect of velocity uncertainty on migrated reflectors: Improvements from relative-depth imaging.  _Geophysics_ , 81(1):S21–S29, 2016.
  * Fomel [2003a] Sergey Fomel.  Time-migration velocity analysis by velocity continuation.  _Geophysics_ , 68(5):1662–1672, 2003a.
  * Fomel [2003b] Sergey Fomel.  Velocity continuation and the anatomy of residual prestack time migration.  _Geophysics_ , 68(5):1650–1661, 2003b.
  * Duchkov and De Hoop [2009] Anton A Duchkov and Maarten V De Hoop.  Velocity continuation in the downward continuation approach to seismic imaging.  _Geophysical Journal International_ , 176(3):909–924, 2009.
  * van Leeuwen and Herrmann [2012] Tristan van Leeuwen and Felix J Herrmann.  Wave-equation Extended Images-Computation and Velocity Continuation.  In _74th EAGE Conference and Exhibition incorporating EUROPEC 2012_ , pages cp–293. European Association of Geoscientists & Engineers, 2012\.
  * Yang et al. [2021] Mengmeng Yang, Marie Graff, Rajiv Kumar, and Felix J Herrmann.  Low-rank representation of omnidirectional subsurface extended image volumes.  _Geophysics_ , 86(3):S165–S183, 2021.
  * Li et al. [2021] Zongyi Li, Nikola Borislavov Kovachki, Kamyar Azizzadenesheli, Burigede liu, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar.  Fourier Neural Operator for Parametric Partial Differential Equations.  In _International Conference on Learning Representations_ , 2021.  URL https://openreview.net/forum?id=c8P9NQVtmnO.
  * Toledo-Marín et al. [2021] J Quetzalcóatl Toledo-Marín, Geoffrey Fox, James P Sluka, and James A Glazier.  Deep learning approaches to surrogates for solving the diffusion equation for mechanistic real-world simulations.  _Frontiers in Physiology_ , 12, 2021.
  * Konuk and Shragge [2021] T Konuk and J Shragge.  Physics-guided deep learning using fourier neural operators for solving the acoustic vti wave equation.  In _82nd EAGE Annual Conference & Exhibition_, volume 2021, pages 1–5. European Association of Geoscientists & Engineers, 2021.
  * Kovachki et al. [2021] Nikola Kovachki, Samuel Lanthaler, and Siddhartha Mishra.  On universal approximation and error bounds for fourier neural operators.  _Journal of Machine Learning Research_ , 22:Art–No, 2021\.
  * Yosinski et al. [2014] Jason Yosinski, Jeff Clune, Yoshua Bengio, and Hod Lipson.  How Transferable Are Features in Deep Neural Networks?  In _Proceedings of the 27th International Conference on Neural Information Processing Systems_ , pages 3320–3328, 2014.  URL http://dl.acm.org/citation.cfm?id=2969033.2969197.
  * Gelman and Rubin [1992] Andrew Gelman and Donald B. Rubin.  Inference from Iterative Simulation Using Multiple Sequences.  _Statistical Science_ , 7(4):457–472, 1992\.  ISSN 08834237.  URL http://www.jstor.org/stable/2246093.
  * Kingma and Ba [2014] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.  URL https://arxiv.org/pdf/1412.6980.pdf.
  * Veritas [2005] Veritas.  Parihaka 3D Marine Seismic Survey - Acquisition and Processing Report.  Technical Report New Zealand Petroleum Report 3460, New Zealand Petroleum & Minerals, Wellington, 2005.
  * WesternGeco. [2012] WesternGeco.  Parihaka 3D PSTM Final Processing Report.  Technical Report New Zealand Petroleum Report 4582, New Zealand Petroleum & Minerals, Wellington, 2012.
  * Wu and Fomel [2018] Xinming Wu and Sergey Fomel.  Least-squares horizons with local slopes and multi-grid correlations.  _GEOPHYSICS_ , 83:IM29–IM40, 2018.  doi: 10.1190/geo2017-0830.1.  URL https://doi.org/10.1190/geo2017-0830.1.
  * Luporini et al. [2018] F. Luporini, M. Lange, M. Louboutin, N. Kukreja, J. Hückelheim, C. Yount, P. Witte, P. H. J. Kelly, F. J. Herrmann, and G. J. Gorman.  Architecture and performance of devito, a system for automated stencil computation.  _CoRR_ , abs/1807.03032, jul 2018.  URL http://arxiv.org/abs/1807.03032.
  * Louboutin et al. [2019] M. Louboutin, M. Lange, F. Luporini, N. Kukreja, P. A. Witte, F. J. Herrmann, P. Velesko, and G. J. Gorman.  Devito (v3.1.0): an embedded domain-specific language for finite differences and geophysical exploration.  _Geoscientific Model Development_ , 12(3):1165–1187, 2019.  doi: 10.5194/gmd-12-1165-2019.  URL https://www.geosci-model-dev.net/12/1165/2019/.