# A Practical Guide to the Numerical Implementation of Tensor Networks I:
Contractions, Decompositions and Gauge Freedom

Glen Evenbly  School of Physics, Georgia Institute of Technology, Atlanta, GA
30332, USA<EMAIL_ADDRESS>
###### Abstract

We present an overview of the key ideas and skills necessary to begin
implementing tensor network methods numerically, which is intended to
facilitate the practical application of tensor network methods for researchers
that are already versed with their theoretical foundations. These skills
include an introduction to the contraction of tensor networks, to optimal
tensor decompositions, and to the manipulation of gauge degrees of freedom in
tensor networks. The topics presented are of key importance to many common
tensor network algorithms such as DMRG, TEBD, TRG, PEPS and MERA.

##  I Introduction

Tensor networks have been developed as a useful formalism for the theoretical
understanding of quantum many-body wavefunctions TN1 ; TN2 ; TN3 ; TN4 ; TN5 ;
TN6 ; TN7 ; TN8 ; TN9 ; TN10 , especially in regards to entanglement Ent1 ;
Ent2 ; Ent3 , and are also applied as powerful numeric tools and simulation
algorithms. Although developed primarily for the description of quantum many-
body systems, they have since found use in a plethora of other applications
such as quantum chemistry Chem1 ; Chem2 ; Chem3 ; Chem4 ; Chem5 , holography
Holo1 ; Holo2 ; Holo3 ; Holo4 ; Holo5 ; Holo6 , machine learning ML1 ; ML2 ;
ML3 ; ML4 ; ML5  and the simulation of quantum circuits QC1 ; QC2 ; QC3 ; QC4
; QC5 ; QC6 .

There currently exist many useful references designed to introduce newcomers
to the underlying theory of tensor networks TN1 ; TN2 ; TN3 ; TN4 ; TN5 ; TN6
; TN7 ; TN8 ; TN9 ; TN10 . Similarly, for established tensor network methods,
there often exist instructional or review articles that address the particular
method in great detail Rev1 ; Rev2 ; Rev3 ; Rev4 ; Rev5 . Nowadays, many
research groups have also made available tensor network code libraries TB1 ;
TB2 ; TB3 ; TB4 ; TB5 ; TB6 ; TB7 ; TB8 . These libraries typically allow
other researchers to make use of highly optimised tensor network routines for
practical purposes (such as for the numerical simulation of quantum many-body
systems).

Comparatively few are resources intended to help researchers that already
possess a firm theoretical grounding to begin writing their own numerical
implementations of tensor network codes. Yet such numerical skills are
essential in many areas of tensor network research: new algorithmic proposals
typically require experimentation, testing and bench-marking using numerics.
Furthermore, even researchers solely interested in the application of tensor
network methods to a problem of interest may be required to program their own
version of a method, as a pre-built package may not contain the necessary
features as to be suitable for the unique problem under consideration. The
purpose of our present work is to help fill this aforementioned gap: to aid
students and researchers, who are assumed to possess some prior theoretical
understanding of tensor networks, to learn the practical skills required to
program their own tensor networks codes and libraries. Indeed, our intent is
to arm the interested reader with the key knowledge that would allow them to
implement their own versions of algorithms such as the density matrix
renormalization group (DMRG) DMRG1 ; DMRG2 ; DMRG3 , time-evolving block
decimation (TEBD) TEBD1 ; TEBD2 , projected entangled pair states (PEPS) PEPS1
; PEPS2 ; PEPS3 , multi-scale entanglement renormalization ansatz (MERA) MERA1
, tensor renormalization group (TRG) TRG1 ; TRG2  or tensor network
renormalization (TNR) TNR1 . Furthermore, this manuscript is designed to
compliment online tensor network tutorials TenNet , which have a focus on code
implementation, with more detailed explanations on tensor network theory.

##  II Preliminaries

###  II.1 Prior knowledge

As stated above, the goal of this manuscript is to help readers that already
possess some understanding of tensor network theory to apply this knowledge
towards numeric calculations. Thus we assume that the reader has some basic
knowledge of tensor networks, specifically that they understand what a tensor
network is and have some familiarity with the standard diagrammatic notation
used to represent them. An overview of these concepts is presented in Fig. 1,
otherwise more comprehensive introductions to tensor network theory can be
found in Refs. [TN1, ; TN2, ; TN3, ; TN4, ; TN5, ; TN6, ; TN7, ; TN8, ; TN9, ;
TN10, ].

Note that we shall not assume prior knowledge of quantum many-body physics,
which is the most common application of tensor network algorithms. The skills
and ideas that we introduce in this manuscript are intended to be general for
the tensor network formalism, rather than for their use in a specific
application, thus can also carry over to other area in which tensor networks
have proven useful such as quantum chemistry Chem1 ; Chem2 ; Chem3 ; Chem4 ;
Chem5 , holography Holo1 ; Holo2 ; Holo3 ; Holo4 ; Holo5 ; Holo6 , machine
learning ML1 ; ML2 ; ML3 ; ML4 ; ML5  and the simulation of quantum circuits
QC1 ; QC2 ; QC3 ; QC4 ; QC5 ; QC6 .

Figure 1: (a-c) Diagrammatic representations of a vector $A_{i}$ (or order-$1$
tensor), a matrix $B_{ij}$ (or order-$2$ tensor) and an order-$3$ tensor
$C_{ijk}$. (d) A contraction, or summation over an index, between two tensors
is represented by a line joining two tensors.

###  II.2 Software Libraries

Currently there exists a wide variety of tensor network code libraries, which
include Refs. TB1, ; TB2, ; TB3, ; TB4, ; TB5, ; TB6, ; TB7, ; TB8, . Many of
these libraries differ greatly in not only their functioning but also their
intended applications, and may have their own specific strengths and
weaknesses (which we will not attempt to survey in the present manuscript).
Almost all of these libraries contain tools to assist in the tasks described
in this manuscript, such as contracting, decomposing and re-gauging tensor
networks. Additionally many of these libraries also contain full featured
versions of complete tensor network algorithms, such as DMRG or TEBD. For a
serious numerical calculation involving tensor networks, one where high
performance is required, most researchers would be well advised to utilise an
existing library.

However, even if ultimate intent is to use existing library, it is still
desirable that one should understand the fundamental tensor network
manipulations used in numerical calculations. Indeed, this understanding is
necessary to properly discern the limitations of various tensor network tools,
to ensure that they are applied in an appropriate way, and to customize the
existing tools if necessary. Moreover, exploratory research into new tensor
network ansatz, algorithms and applications often requires non-standard
operations and tensor manipulations which may not be present in any existing
library, thus may require the development of extensive new tensor code. In
this setting it can be advantageous to minimize or to forgo the usage of an
existing library (unless one was already intimately familiar with its inner
workings), given the inherent challenge of extending a library beyond its
intended function and the possibility of unintended behavior that this
entails.

In the remaining manuscript we aim to describe key tensor network operations
(namely contracting, decomposing and re-gauging tensor networks) with
sufficient detail that would allow the interested reader to implement tasks
numerically without the need to rely on an existing code library.

###  II.3 Programming language

Before attempting to implement tensor methods numerically one must, of course,
decide on which programming language to use. High-level languages with a focus
on scientific computation, such as MATLAB, Julia and Python (with Numpy) are
well-suited for implementing tensor network methods as they have native
support for multi-dimensional arrays (i.e. tensors), providing simple and
convenient syntax for common operations on these arrays (indexing, slicing,
scalar operations) as well as providing a plethora of useful functions for
manipulating these tensor objects. Alternatively, some tensor network
practitioners may prefer to use lower-level languages such as Fortran or C++
when implementing tensor network algorithms usually for the reason of
maximising performance. However, in many tensor network codes the bulk of the
computation time is spent performing large matrix-matrix multiplications, for
which even interpreted languages (like MATLAB) still have competitive
performance with compiled languages as they call the same underlying BLAS
routines. Nonetheless, there are some particular scenarios, such as in dealing
with tensor networks in the presence of global symmetries Sym1 ; Sym2 ; Sym3 ;
Sym4 ; Sym5 ; Sym6 ; Sym7 , where a complied language may offer a significant
performance advantage. In this circumstance Julia, which is a compiled
language, or Python, in conjunction with various frameworks which allow it to
achieve some degree of compilation, may be appealing options.

###  II.4 Terminology

Before proceeding, let us establish the terminology that we will use when
discussing tensor networks. We define the _order_ of a tensor as the number of
indices it has, i.e. such that a vector is order-$1$ and a matrix is
order-$2$. The term _rank_ (or decomposition rank) of a tensor will refer to
the number of non-zero singular values with respect to a some bi-partition of
the tensor indices. Note that many researchers also use the term _rank_ to
describe the number of tensor indices; here we use the alternative term
_order_ specifically to avoid the confusion that arises from the double usage
of _rank_. The number of values a tensor index can take will be referred to as
the dimension of an index (or bond dimension), which is most often denoted by
$\chi$ but can also be denoted by $m$, $d$ or $D$. In most cases, the use of
$d$ or $D$ to denote a bond dimension is less preferred, as this can be
confused which the spatial dimension of the problem under consideration (e.g.
when considering a model on a $1D$ or $2D$ lattice geometry).

Figure 2: (a) The internal indices $(l,m,n)$ of the network $\\{A,B,C\\}$ are
contracted to give tensor $F$. (b) The network is contracted via a sequence of
two pairwise tensor contractions, the first of which results in the
intermediate tensor $D$.

##  III Tensor Contractions

The foundation of all tensor networks routines is the contraction of a network
containing multiple tensors into a single tensor. An example of the type
problem that we consider is depicted in Fig. 2(a), where we wish to contract
the network of tensors $\\{A,B,C\\}$ to form an order-$3$ tensor $F$, which
has components defined

${F_{ijk}}=\sum\limits_{l,m,n}{{A_{ljm}}{B_{iln}}{C_{nmk}}}.$ (1)

Note that a convention for tensor index ordering is required for the figure to
be unambiguously translated to an equation; here we assumed that indices
progress clock-wise on each tensor starting from 6 o’clock. Perhaps the most
obvious way to evaluate Eq. 1 numerically would be through a direct summation
over the indices $(l,m,n)$, which could be implemented using a set of nested
‘FOR’ loops. While this approach of summing over all internal indices of a
network will produce the correct answer, there are numerous reasons why this
is not the preferred approach for evaluating tensor networks. The foremost
reason is that it is not the most computationally efficient approach
(excluding, perhaps, certain contractions involving sparse tensors, which we
will not consider here). Let us analyse the contraction cost for the example
given in Eq. 1, assuming all tensor indices are $\chi$-dimensional. A single
element of tensor $F$, which has $\chi^{3}$ elements in total, is given
through a sum over indices $(l,m,n)$, which requires $O(\chi^{3})$ operations.
Thus the total cost of evaluating tensor $F$ through a direct summation over
all internal indices is $O(\chi^{6})$.

Now, let us instead consider the evaluation of tensor $F$ broken up into two
steps, where we first compute an intermediate tensor $D$ as depicted in
Fig.2(b),

${D_{ijmn}}=\sum\limits_{l}{{A_{ljm}}{B_{iln}}},$ (2)

before performing a second contraction to give the final tensor $F$,

${F_{ijk}}=\sum\limits_{n,m}{{D_{ijmn}}}{C_{nm}}.$ (3)

Through similar logic as before, one finds that the cost of evaluating
intermediate tensor $D$ scales as $O(\chi^{5})$, whilst the subsequent
evaluation of $F$ in Eq. 3 is also $O(\chi^{5})$. Thus breaking the network
contraction down into a sequence of smaller contractions each only involving a
pair of tensors (which we refer to as a _pairwise tensor contraction_) is as
computationally cheap or cheaper for any non-trivial bond dimension
$(\chi>1)$. This is true in general: for any network of 3 or more (dense)
tensors it is always at least as efficient (and usually vastly more efficient)
to break network contraction into sequence of pairwise contractions, as
opposed to directly summing over all the internal indices of the network.

Two natural questions arise at this point. (i) What is optimal way to
implement a single pairwise tensor contraction? (ii) Does the chosen sequence
of pairwise contractions affect the total computational cost and, if so, how
does one decide what sequence to use? We begin by addressing the first
question.

###  III.1 Pairwise tensor contractions

Let us consider the problem of evaluating a pairwise tensor contraction,
denoted $(A\times B)$, between tensors $A$ and $B$ that are connected by one
or more common indices. A straight-forward method to evaluate such
contractions, as in the examples of Eq. 2 and Eq. 3, is by using nested ‘FOR’
loops to sum over the shared indices. The computational cost of this
evaluation, in terms of the number of scalar multiplications required, can be
expressed concisely as

$\textrm{cost}:(A\times
B)=\frac{{\left|{\dim(A)}\right|\cdot\left|{\dim(B)}\right|}}{{\left|{\dim(A\cap
B)}\right|}},$ (4)

with $\left|{\dim(A)}\right|$ as the total dimension of $A$ (i.e. the product
of its index dimensions) and $\left|{\dim(A\cap B)}\right|$ as the total
dimension of the shared indices.

Alternatively, one can recast a pairwise contraction as a matrix
multiplication as follows: first reorder the free indices and contracted
indices on each of $A$ and $B$ such that they appear sequentially (which can
be achieved in MATLAB using the ‘permute’ function) and then group the free-
indices and the contracted indices each into a single index (which can be
achieved in MATLAB using the ‘reshape’ function). After these steps the
contraction is evaluated using a single matrix-matrix multiplication, although
the final product may also need to be reshaped back into a tensor. Recasting
as a matrix multiplication does not reduce the formal computational cost from
Eq. 4. However, modern computers, leveraging highly optimized BLAS routines,
typically perform matrix multiplications significantly more efficiently than
the equivalent ‘FOR’ loop summations. Thus, especially in the limit of tensors
with large total dimension, recasting as a matrix multiplication is most often
the preferred approach to evaluate pairwise tensor contractions, even though
this requires some additional computational overhead from the necessity of
rearranging tensor elements in memory when using ‘permute’. Note that the
‘tensordot’ function in the Numpy module for Python conveniently evaluates a
pairwise tensor contraction using this matrix multiplication approach.

Figure 3: (a) A product of three tensors $\\{A,B,C\\}$ is contracted to a
tensor $F$, where all indices are assumed to be of dimension $\chi$. (b-c) The
total computational cost of contracting the network depends on the sequence of
pairwise contractions; the cost from following the sequence in (b) scales as
$(\chi^{3}+\chi^{2})$ as compared to the cost from (c) which scales as
$(2\chi^{2})$.

###  III.2 Contraction sequence

It is straight-forward to establish that, when breaking a network contraction
into a sequence of binary contractions, the choice of sequence can affect the
total computational cost. As an example, we consider the product of two
matrices $A$, $B$ with vector $C$,

${F_{i}}=\sum\limits_{j,k}{{A_{ij}}{B_{jk}}{C_{k}}},$ (5)

where all indices are assumed to be dimension $\chi$, see also Fig. 3. If we
evaluate this expression by first performing the matrix-matrix multiplication,
i.e. as $F=(A\times B)\times C$, then the leading order computational cost
scales as $O(\chi^{3})$ by Eq. 4. Alternatively, if we evaluate the expression
by first performing the matrix-vector multiplication, i.e. as
$F=A\times(B\times C)$, then the leading order computational cost scales as
$O(\chi^{2})$. Thus it is evident that the sequence of binary contractions
needs to be properly considered in order to minimize the overall computational
cost.

So how does one find the optimal contraction sequence for some tensor network
of interest? For the networks that arise in common algorithms (such as DMRG,
MERA, PEPS, TRG and TNR) it is relatively easy, with some practice, to find
the optimal sequence through manual inspection or trial-and-error. This
follows as most networks one needs to evaluate contain fewer than 10 tensors
and the tensor index dimensions take a only single or a few distinct values
within a network, which limits the number of viable contraction sequences that
need be considered. More generally, determination of optimal contraction
sequences is known to be an NP-hard problem Order1 , such that it is very
unlikely that an algorithm which scales polynomially with the number of
tensors in the network will ever be found to exist. However, numerical methods
based on exhaustive searches and/or heuristics can typically find optimal
sequences for networks with fewer than 20 tensors in a reasonable amount of
time Order1 ; Order2 ; Order4 ; Order5 , and larger networks are seldom
encountered in practice.

Note that many tensor network optimisation algorithms based on an iterative
sweep, where the same network diagrams are contracted each iteration (although
perhaps containing different tensors and with different bond dimensions). The
usual approach in this setting is to determine the optimal sequences once,
before beginning the iterative sweeps, using the initial bond dimensions and
then cache the sequences for re-use in later iterations. The contraction
sequences are then only recomputed the if the bond dimensions stray too far
from the initial values.

###  III.3 Network contraction routines

Although certainly feasible, manually writing the code for each tensor network
contraction as a sequence of pairwise contractions is not recommended. Not
only is substantial programming effort required, but this also results in code
which is error-prone and difficult to check. There is also a more fundamental
problem: contracting a network by manually writing a sequence of binary
contractions requires specifying a particular contraction sequence at the time
of coding. However, in many cases the index dimensions within networks are
variable, and the optimal sequence can change depending on the relative sizes
of dimensions. For instance, one may have a network which contains indices of
dimensions $\chi_{1}$ and $\chi_{2}$, where the optimal contraction sequence
changes dependant of whether $\chi_{1}$ is larger or smaller than $\chi_{2}$.
In order to have a program which works efficiently in both regimes, one would
have to write code separately for both contraction sequences.

Given the considerations above, the use of an automated contraction routine,
such as the ‘ncon’ (Network-CONtractor) function Order3 ; TenNet  or something
similar from an existing tensor network library TB1 ; TB2 ; TB3 ; TB4 ; TB5 ;
TB6 ; TB7 ; TB8 , is highly recommended. Automated contraction routines can
evaluate any network contraction in a single call by appropriately generating
and evaluating a sequence of binary contractions, hence greatly reducing both
the programming effort required and the risk of programming errors occurring.
Most contraction routines, such as ‘ncon’, also remove the need to fix a
contraction sequence at the time of writing the code, as the sequence can be
specified as an input variable to the routine and thus can be changed without
the need to rewrite any code. This can also allow the contraction sequence to
be adjusted dynamically at run-time to ensure that the sequence is optimal
given the specific index dimensions in use.

###  III.4 Summary: contractions

In evaluating a network of multiple tensors, it is always more efficient to
break the contraction into a sequence of pairwise tensor contractions, each of
which should (usually) be recast into a matrix-matrix multiplication in order
to achieve optimal computational performance. The total cost of evaluating a
network can depend on the particular sequence of pairwise contractions chosen.
While there is no known method for determining an optimal contraction sequence
that is efficiently scalable in the size of the network, manual inference or
brute-force numeric searches are usually viable for the relatively small
networks encountered in common tensor network algorithms. When coding a tensor
network program it is useful to utilise an automated network contraction
routine which can evaluate a tensor network in a single call by properly
chaining together a sequence of pairwise contractions. This not only reduces
the programming effort required, but also grants a program more flexibility in
allowing a contraction sequence to be easily changed.

Figure 4: Depiction of some common matrix decompositions. All indices are
assumed to be of dimension $d$ unless otherwise indicated. (a-i) The spectral
decomposition is applied to the order-4 Hermitian tensor $H$ across the
partition indicated by the dashed line, yielding a diagonal matrix of
eigenvalues $D$ and a unitary $U$. (a-ii) The unitary tensor $U$ annihilates
to identity with its conjugate, as per Eq.6. (b-i) The QR decomposition is
applied to the order-3 tensor $A$ across the partition indicated, yielding an
isometry $Q$ and an upper triangular matrix $R$. (b-ii) The isometry $Q$
annihilates to identity with its conjugate as per Eq.7, while the $R$ matrix
is upper triangular. (c-i) The singular value decomposition (SVD) is applied
to the order-3 tensor $A$ across the partition indicated, yielding an isometry
$U$, a diagonal matrix of singular values $S$, and a unitary $V$. (c-ii)
Depiction of the constraints satisfied by the isometry $U$ and unitary $V$.

##  IV Matrix Factorizations

Another key operation common in tensor network algorithms, complimentary to
the tensor contractions considered previously, are factorizations. In this
section we will discuss some of the various means by which a higher-order
tensor can be split into a product of fewer-order tensors. In particular, the
means that we consider involve applying standard matrix decompositionsMat1 ;
Mat2 , to tensor unfoldings, such that this section may serve as a review of
the linear algebra necessary before consideration of more sophisticated
network decompositions. Specifically we recount the spectral decomposition, QR
decomposition and singular value decomposition and outline their usefulness in
the context of tensor networks, particular in achieving optimal low-rank
tensor approximations. Before discussing decompositions, we define some
special types of tensor.

###  IV.1 Special tensor types

A $d$-by-$d$ matrix $U$ is said to be unitary if it has orthonormal rows and
columns, which implies that it annihilates to the identity when multiplied
with its conjugate,

${U^{\dagger}}U=U{U^{\dagger}}=I,$ (6)

where $I$ is the $d$-by-$d$ identity matrix. We define a tensor (whose order
is greater than 2) as unitary with respect to a particular bi-partition of
indices if the tensor can be _reshaped_ into a unitary matrix according to
this partition. Similarly an $d_{1}$-by-$d_{2}$ matrix $W$, with $d_{1}>d_{2}$
is said to be an isometry if

${W^{\dagger}}W=I,$ (7)

with $I$ the $d_{2}$-by-$d_{2}$ identity matrix. Likewise we say that a tensor
(order greater than 2) is isometric with respect to a particular bi-partition
of indices if the tensor can be reshaped into a isometric matrix. Notice that,
rather than equalling identity, the reverse order product does now evaluate to
a projector $P$,

$W{W^{\dagger}}=P,$ (8)

where projectors are defined as Hermitian matrices that square to themselves,

$P=P^{\dagger},\;\;P^{2}=P.$ (9)

###  IV.2 Useful matrix decompositions:

A commonly used matrix decomposition is the spectral decomposition (or eigen-
decomposition). In the context of tensor network codes it is most often used
for Hermitian, positive semi-definite matrices, such as for the density
matrices used to describe quantum states. If $H$ is a $d\times d$ Hermitian
matrix, or tensor that can be reshaped into such, then the spectral
decomposition yields

$H=UD{U^{\dagger}},$ (10)

where $U$ is $d\times d$ unitary matrix and $D$ is diagonal matrix of
eigenvalues, see also Fig.4(a). The numerical cost of performing the
decomposition scales as $O(d^{3})$. In the context of tensor network
algorithms the spectral decomposition is often applied to approximate a
Hermitian tensor with one of smaller rank, as will be discussed in Sect. IV.4.

Another useful decomposition is the QR decomposition. If $A$ be an arbitrary
$d_{1}\times d_{2}$ matrix with $d_{1}>d_{2}$, or tensor that can be reshaped
into such, then the QR decomposition gives

$A=QR,$ (11)

see also Fig.4(b). Here $Q$ is $d_{1}\times d_{2}$ isometry, such that
$Q^{\dagger}Q=I$, where $I$ is the $d_{2}\times d_{2}$ identity matrix, and
$R$ is $d_{2}\times d_{2}$ upper triangular matrix. Note that we are
considering the so-called economical decomposition (which is most often used
in tensor network algorithms); otherwise the full decomposition gives $Q$ as a
$d_{1}\times d_{1}$ unitary and $R$ is dimension $d_{1}\times d_{2}$. The
numerical cost of the economical QR decomposition scales as the larger
dimension times the square of the smaller dimension $O(d_{1}{d_{2}}^{2})$, as
opposed to cost $O({d_{1}}^{2}d_{2})$ for the full decomposition. The QR
decomposition is one of the most computationally efficient ways to obtain an
orthonormal basis for the column space of a matrix, thus a common application
is in orthogonalizing tensors within a network (i.e. transforming them into
isometries), which will be discussed further in Sect. V.3.1.

The final decomposition that we consider is the singular value decomposition
(SVD), which is also widely used in many areas of mathematics, statistics,
physics and engineering. The SVD allows an arbitrary $d_{1}\times d_{2}$
matrix $A$, where we assume for simplicity that $d_{1}\geq d_{2}$, to be
decomposed as

$A=USV^{\dagger}$ (12)

where $U$ is $d_{1}\times d_{2}$ isometry (or unitary if $d_{1}=d_{2}$), $S$
is diagonal $d_{2}\times d_{2}$ matrix of positive elements (called singular
values), and $V$ is $d_{2}\times d_{2}$ unitary matrix, see also Fig.4(c).
Similar to the economical QR decomposition, we have also considered the
economical form of the SVD; the full SVD would otherwise produce $U$ as a
$d_{1}\times d_{1}$ unitary and $S$ as a rectangular $d_{1}\times d_{2}$
matrix padded with zeros. The numerical cost of the economical SVD scales as
$O(d_{1}{d_{2}}^{2})$, identical to that of the economical QR decomposition.
The rank of a tensor (across a specified bi-partition) is defined as the
number of non-zero singular values that appear in the SVD. A common
application of the SVD is in finding an approximation to a tensor by another
of smaller rank, which will be discussed further in Sect. IV.4.

Notice that for any matrix $A$ the spectral decompositions of $AA^{\dagger}$
and $A^{\dagger}A$ are related to the SVD of $A$; more precisely, the
eigenvectors of $AA^{\dagger}$ and $A^{\dagger}A$ are equivalent to the
singular vectors in $U$ and $V$ respectively of the SVD in Eq. 12. Furthermore
the (non-zero) eigenvalues in $AA^{\dagger}$ or $A^{\dagger}A$ are the squares
of the singular values in $S$. It can also be seen that, for a Hermitian
positive semi-definite matrix $H$, the spectral decomposition is equivalent to
an SVD.

Figure 5: For any tensor $A$ the tensor trace Ttr of $A$ with its conjugate
$A^{\dagger}$ (drawn with opposite vertical orientation) is obtained by
contracting over all matching indices. The Frobenius norm can be defined as
the root of this tensor trace, see Eq. 14.

###  IV.3 Tensor Norms

The primary use for matrix decompositions, such as the SVD, in the context of
tensor networks is in accurately approximating a higher-order tensor as a
product lower-order tensors. However, before discussing tensor approximations,
it is necessary to define the tensor norm in use. A tensor norm that is
particularly convenient is the Frobenius norm (or Hilbert-Schmidt norm). Given
a tensor $A_{ijk\ldots}$ the Frobenius norm for $A$, denoted as
$\left\|A\right\|$, is defined as the square-root of the sum of the magnitude
of each element squared,

$\left\|A\right\|=\sqrt{\sum\limits_{ijk\ldots}{{{\left|{{A_{ijk\ldots}}}\right|}^{2}}}}.$
(13)

This can be equivalently expressed as the tensor trace of $A$ multiplied by
its conjugate,

$\left\|A\right\|=\sqrt{\textrm{Ttr}\left({{A^{\dagger}},A}\right)},$ (14)

where the tensor trace, $\textrm{Ttr}(A^{\dagger},A)$, represents the
contraction of tensor $A$ with its conjugate over all matching indices, see
Fig. 5. It also follows that Frobenius norm is related to the singular values
$s_{k}$ of $A$ across any chosen bi-partition,

$\left\|A\right\|=\sqrt{\sum\limits_{k}{{{({s_{k}})}^{2}}}}.$ (15)

Notice that Eq. 14 implies that the difference
${\varepsilon}=\left\|{A-B}\right\|$ between two tensors $A$ and $B$ of equal
dimension can equivalently be expressed as

${\left\|{A-B}\right\|^{2}}=\textrm{Ttr}({A^{\dagger}},A)-2\left|{\textrm{Ttr}({A^{\dagger}},B)}\right|+\textrm{Ttr}({B^{\dagger}},B).$
(16)

Figure 6: (a) The singular value decomposition is taken on tensor $A$ across a
bi-partition between its top two and bottom three indices, and is assumed to
produce $d$ non-zero singular values. (b) Tensor $B$ is now defined by
truncating the matrix of singular values $S\rightarrow\tilde{S}$ to retain
only $r<d$ of the largest singular values, while similarly truncating the
matrices of singular vectors, $U\rightarrow\tilde{U}$ and
$V\rightarrow\tilde{V}$, to retain only the corresponding singular vectors. By
the Eckart–Young–Mirsky theorem Low  it is known that $B$ is the optimal
rank-$r$ approximation to $A$ (across the chosen bi-partition of tensor
indices).

###  IV.4 Optimal Low-rank Approximations

Given some matrix $A$, or higher-order tensor that viewed as a matrix across a
chosen bi-partition of its indices, we now focus on the problem of finding the
tensor $B$ that best approximates $A$ according to the Frobenius norm (i.e.
that which minimizes the difference in Eq. 16), assuming $B$ has a fixed rank
$r$. Let us assume, without loss of generality, that tensor $A$ is equivalent
to a $d_{1}\times d_{2}$ matrix (with $d_{1}\geq d_{2}$) under a specified bi-
partition of its indices, and that $A$ has singular value decomposition,

${A_{ij}}=\sum\limits_{k=1}^{d_{2}}{{U_{ik}}{s_{k}}V_{kj}^{*}},$ (17)

where the singular values are assumed to be in descending order, $s_{k}\geq
s_{k+1}$. Then the optimal rank $r$ tensor $B$ that approximates $A$ is known
from the Eckart–Young–Mirsky theorem Low , which states that $B$ is given by
truncating to retain only the $r$ largest singular values and their
corresponding singular vectors,

$B_{ij}=\sum\limits_{k=1}^{r}U_{ik}s_{k}V_{kj}^{*}.$ (18)

It follows that the error of approximation $\varepsilon=\left\|A-B\right\|$,
as measured in the Frobenius norm, is related to the discarded singular values
as

$\varepsilon=\sqrt{\sum\limits_{k>r}{{{({s_{k}})}^{2}}}}.$ (19)

If the spectrum of singular values is sharply decaying then the error is well
approximated by the largest of the discarded singular values,
$\varepsilon\approx s_{(r+1)}$.

Notice that, in the case that the tensor $A$ under consideration is Hermitian
and positive definite across the chosen bi-partition, that the spectral
decomposition could instead be used in Eq. 17. The low-rank approximation
obtained by truncating the smallest eigenvalues would still be guaranteed
optimal, as the spectral decomposition is equivalent to the SVD in this case.

###  IV.5 Summary: decompositions

In this section we have described how special types of matrices, such as
unitary matrices and projections, can be generalized to the case of tensors
(which can always be viewed as a matrix across a chosen bi-partition of their
indices). Similarly we have shown how several concepts from matrices, such as
the trace and the norm, are also generalized to tensors. Finally, we have
described how the optimal low-rank approximation to a tensor can be obtained
via the SVD.

##  V Gauge Freedom

All tensor networks possess gauge degrees of freedom; exploiting the gauge
freedom allows one to change the tensors within a network whilst leaving the
final product of the network unchanged. In this section we describe methods
for manipulating the gauge degrees of freedom and discuss the utility of
fixing the gauge in a specific manner.

###  V.1 Tree tensor networks

In this manuscript we shall restrict our considerations of gauge manipulations
to focus only on tensors networks that do not possess closed loops (i.e.
networks that correspond to acyclic graphs), which are commonly referred to as
tree tensor networks (TTN) TTN1 ; TTN2 . Fig. 7(a) presents an example of a
tree tensor network. If we select a single tensor to act as the center (or
root node) then we can understand the tree tensor network as being composed of
a set of distinct branches extending from this chosen tensor. For instance,
Fig. 7(b) depicts the three branches (excluding the single trivial branch)
extending from the $4^{\textrm{th}}$-order tensor $A$. It is important to note
that connections between the different branches are not possible in networks
without closed loops; this restriction is required for the re-gauging methods
considered in this manuscript. However these methods can (mostly) be
generalized to the case of networks containing closed loops by using a more
sophisticated formalism as shown in Ref. Gauge, .

Figure 7: (a) An example of a tree tensor network (TTN), here composed of 7
tensors. (b) The three (non-trivial) branches of the tree with respect to
choosing $A$ as the root tensor.

###  V.2 Gauge transformations

Consider a tensor network of multiple tensors that, under contraction of all
internal indices, evaluates to some product tensor $H$. We now pose the
following question: is there a different choice of tensors with the same
network geometry that also evaluates to $H$? Clearly the answer to this
question is yes! As shown below in Fig. 8(a), on any of the internal indices
of the network one can introduce a resolution of the identity (i.e. a pair of
matrices $X$ and $X^{-1}$) which, by construction, does not change the final
product that the network evaluates to. However, absorbing one of these
matrices into each adjoining tensor changes the individual tensors while
leaving the product of the network unchanged. It follows that there are
infinitely many choices of tensors such that the network product evaluates to
some fixed output tensor, since the gauge change matrix $X$ can be any
invertible matrix. This ability to introduce an arbitrary resolution of the
identity on an internal index, while leaving the product of the network
unchanged, is referred to as the gauge freedom of the network.

While in some respects the gauge freedom could be considered bothersome, as it
implies tensor decompositions are never unique, it can also be exploited to
simplify many types of operations on tensor networks. Indeed, most tensor
network algorithms require fixing the gauge in a prescribed manner in order to
function correctly. In the following sections we discuss ways to fix the gauge
degree of freedom as to create an _orthogonality center_ and the benefits of
doing so.

Figure 8: (a) Given a network of three tensors $\\{A,B,C\\}$, one can
introduce gauge change matrices $X$ and $Y$ (together with their inverses) on
the internal indices of the network while leaving the final product $D$ of the
network unchanged. (b) Definitions of the new tensors
$\\{\tilde{A},\tilde{B},\tilde{C}\\}$ after the change of gauge. Figure 9: (a)
An example of a tree tensor network. (b) A depiction of the constraints
required for the tensor $A$ to be an _orthogonality center_ : each of the
branches must annihilate to the identity when contracted with its conjugate.

###  V.3 Orthogonality centers

A given tensor $A$ within a network is said to be an _orthogonality center_ if
every branch connected to tensor $A$ annihilates to the identity when
contracted with its conjugate as shown in Fig. 9. Equivalently, each branch
must (collectively) form an isometry across the bi-partition between its open
indices and the index connected to tensor $A$. By properly manipulating the
gauge degrees of freedom, it is possible to turn any tensor with a tree tensor
network into an orthogonality centerOrth1 . We now discuss two different
methods for achieving this: a ‘pulling through’ approach, which was a key
ingredient in the original formulation of DMRGDMRG1 ; DMRG2 ; DMRG3 , and a
‘direct orthogonalization’ approach, which was an important part of the TEBD
algorithmTEBD1 ; TEBD2 .

Figure 10: A depiction of how the tensor $A$ can be made into an orthogonality
center of the network from Fig. 7 via the “pulling-though” approach. (a-b)
Tensors $D$ and $E$, which reside at the tips of a branch, are decomposed via
the QR decomposition. (c) The $R$ components of the previous QR decompositions
are absorbed into the $B$ tensor higher on the branch, which is then itself
decomposed via the QR decomposition. (d) Following this procedure, all tensors
in the network are orthogonalized (with respect to their incoming versus
outgoing indices) such that $A^{\prime}$ becomes an orthogonality center of
the network.

####  V.3.1 Creating an orthogonality center via ‘pulling through’

Here we describe a method for setting a tensor $A$ within a network as an
orthogonality center through iterative use of the QR decomposition. The idea
behind his method is very simple: if each individual tensor within a branch is
transformed into a (properly oriented) isometry, then the entire branch
collectively becomes an isometry and thus the tensor at center of the branches
becomes an orthogonality center. Let us begin by orienting each index of the
network by drawing an arrow pointing towards the desired center $A$. Then,
starting at the tip of each branch, we should perform a QR decomposition on
each tensor based on a bi-partition between its incoming and outgoing indices.
The $R$ part of the decomposition should then be absorbed into the next tensor
in the branch (i.e. closer to the root tensor $A$), and the process repeated
as depicted in Fig. 10(a-c). At the final step an $R$ part of the QR
decomposition from each branch is absorbed into the central tensor $A$ and the
process is complete, see also Fig. 10(d).

Note that the SVD could be used as an alternative to the QR decomposition:
instead of absorbing the $R$ part of the QR decomposition into the next tensor
in the branch one could absorb the product of the $S$ and $V$ part of the SVD
from Eq. 12. However, in practice, the QR decomposition is most often
preferable as it computationally cheaper than the SVD.

Figure 11: A depiction of how the tensor $A$ from the network of Fig. 7 can be
made into an orthogonality center via the “direct orthogonalization” approach.
(a) A change of gauge is made on every (non-trivial) branch connected to $A$,
such that $A$ becomes an orthogonality center. (b-d) The gauge change matrices
$\\{X_{1},X_{2},X_{3}\\}$ are obtained by contracting each branch with its
Hermitian conjugate and then taking the principle root.

####  V.3.2 Creating an orthogonality center via ‘direct orthogonalization’

Here we describe a method for setting a tensor $A$ within a network as an
orthogonality center based on use of a single decomposition for each branch,
as depicted in Fig. 11. (i) We begin by computing the positive-definite matrix
$\rho$ for each branch (with respect to the center tensor $A$) by contracting
the branch with its Hermitian conjugates. (ii) The principle square root $X$
of each of the matrices $\rho$ is then computed, i.e. such that
$\rho=XX^{\dagger}$, which can be computed using the spectral decomposition if
necessary. (iii) Finally, a change of gauge is made on each of the indices of
tensor $A$ using the appropriate $X$ matrix and its corresponding inverse
$X^{-1}$, with the $X$ part absorbed in tensor $A$ and the $X^{-1}$ absorbed
in the leading tensor of the branch such that the branch matrix transforms as
$\rho\rightarrow\tilde{\rho}=X^{-1}\rho(X^{-1})^{\dagger}$. It follows that
the tensor $A$ is now an orthogonality center as each branch matrix
$\tilde{\rho}$ of the transformed network evaluates as the identity,

$\tilde{\rho}=X^{-1}\rho(X^{-1})^{\dagger}=X^{-1}XX^{\dagger}(X^{-1})^{\dagger}=I.$
(20)

Note that, for simplicity, we have assumed that the branch matrices $\rho$ do
not have zero eigenvalues, such that their inverses exist. Otherwise, if zero
eigenvalues are present, the current method is not valid unless the index
dimensions are first reduced by truncating any zero eigenvalues.

####  V.3.3 Comparison of methods for creating orthogonality centers

Each of the two methods discussed to create an orthogonality center have their
own advantages and disadvantages, such that the preferred method may depend on
the specific application under consideration. In practice, the ‘direct
orthogonalization’ is typically computationally cheaper and easier to execute,
since this method only requires changing the gauge on the indices connected to
the center tensor. In contrast the ‘pulling through’ method involves changing
the gauge on all indices of the network. Additionally, the ‘direct
orthogonalization’ approach can easily be employed in networks of infinite
extent, such as infinite MPS TEBD1 ; TEBD2 , if the matrix $\rho$ associated
to a branch of infinite extent can be computed using by solving for a dominant
eigenvector. While ‘pulling through’ can also potentially be employed for
networks of infinite extent, i.e. through successive decompositions until
sufficient convergence is achieved, this is likely to be more computationally
expensive. However the ‘pulling through’ approach can be advantageous if the
branch matrices $\rho$ are ill-conditioned as the errors due to floating-point
arithmetic are lesser. This follows since the ‘direct orthogonalization’
requires one to square the tensors in each branch. The ‘pulling-through’
approach also results in transforming every tensor in the network (with the
exception of the center tensor) into an isometry, which may be desirable in
certain applications.

Figure 12: (a) A network of 7 tensors $\\{A,B,C,D,E,F,G\\}$ contracts to give
a tensor $H$. (b) After replacing a single tensor $A\rightarrow A^{\prime}$
the network contracts to a different tensor $H^{\prime}$. (c) The tensor
$A^{\prime}$ is decomposed into a pair of tensors $A_{L}$ and $A_{R}$, leaving
the final tensor $H^{\prime}$ unchanged.

###  V.4 Decompositions of tensors within networks

In Sect. IV it was described how the SVD could be applied to find the optimal
low-rank approximation to a tensor in terms of minimizing the Frobenius norm
between the original tensor and the approximation. In the present section we
extend this concept further and detail how, by first creating an orthogonality
center, a tensor within a network can be optimally decomposed as to minimize
the _global_ error coming from consideration of the entire network.

Let us consider a tree tensor network of tensors $\\{A,B,C,D,E,F,G\\}$ that
evaluates to a tensor $H$, as depicted in Fig. 12(a). We now replace a single
tensor $A$ from this network by a new tensor $A^{\prime}$ such that the
network now evaluates to a tensor $H^{\prime}$ as depicted in Fig. 12(b). Our
goal is to address the following question: how can we find the optimal low-
rank approximation $A^{\prime}$ to tensor $A$ such that the error from the
full network, $\|H-H^{\prime}\|$, is minimised? Notice that if we follow the
method from Sect. IV and simple truncate the smallest singular values of $A$,
then this will only ensure that the local error, $\|A-A^{\prime}\|$, is
minimised. The key to resolving this issue is through creation of an
orthogonality center, which can reduce the global norm of a network to the
norm of a single tensor. Specifically if tensor $A$ is an orthogonality center
of a network that evaluates to a final tensor $H$ then it follow from the
definition of an orthogonality center that $\|H\|=\|A\|$, as depicted in Fig.
13(a). Thus it also can be seen that under replacement of the center tensor
$A$ with a new tensor $A^{\prime}$, such that the network now evaluates to a
new tensor $H^{\prime}$, that the difference between the tensors
$\|A-A^{\prime}\|$ is precisely equal to the global difference between the
networks $\|H-H^{\prime}\|$. This follows as the overlap of $H$ and
$H^{\prime}$ equals the overlap of $A$ and $A^{\prime}$, as depicted in Fig.
13(b). In other words, by appropriately manipulating the gauge degrees of
freedom in a network, the _global_ difference resulting from changing a single
tensor in a network can become equivalent to the _local_ difference between
the single tensors. The solution to the problem of finding the optimal low-
rank approximation $A^{\prime}$ to a tensor $A$ within a network thus becomes
clear; we should first adjust the gauge such that $A$ becomes an orthogonality
center, after which we can follow the method from Sect. IV and create the
optimal global approximation (i.e. that which minimizes the global error) by
truncating the smallest singular values of $A$. The importance of this result
in the context tensor network algorithms cannot be overstated; this
understanding for how to optimally truncate a single tensor within a tensor
network, see also Ref.CanForm, , is a key aspect of the DMRG algorithm DMRG1 ;
DMRG2 ; DMRG3 , the TEBD algorithm TEBD1 ; TEBD2  and many other tensor
network algorithms.

Figure 13: (a) In the network from Fig. 12(a), if the tensor $A$ is an
orthogonality center then it follows that the norm of the network $\|H\|$ is
equal to the norm of the center tensor $\|A\|$. (b) Similarly it follows that,
in changing only the center tensor $A\rightarrow A^{\prime}$, the global
overlap between the networks $H$ and $H^{\prime}$ is equal to the local
overlap between the center tensors $A$ and $A^{\prime}$.

###  V.5 Summary: gauge freedom

In the preceding section we discussed manipulations of the gauge degrees of
freedom in a tensor network and described two methods that can be used to
create an orthogonality center. The proper use of an orthogonality center was
then demonstrated to allow one to decompose a tensor within a network in such
a way as to minimize the global error. Note that while the results in this
section were described only for tree tensor networks (i.e. networks based on
acyclic graphs), they can be generalized to arbitrary networks by using more
sophisticated methodology Gauge .

##  VI Conclusions

Network contractions and decompositions are the twin pillars of all tensor
network algorithms. In this manuscript we have recounted the key theoretical
considerations required for performing these operations efficiently and also
discussed aspects of their implementation in numeric codes. We expect that a
proper understanding of these results could facilitate an individuals effort
to implement many common tensor network algorithms, such as DMRG, TEBD, TRG,
PEPS and MERA, and also further aid researchers in the design and development
of new tensor network algorithms.

However, there are still a wide variety of additional general ideas and
methods, not covered in this manuscript, that are necessary for the
implementation of more advanced tensor network algorithms. These include (i)
strategies for performing variational optimization, (ii) methods for dealing
with decompositions in networks containing closed loops, (iii) the use of
approximations in tensor network contractions. We shall address several of
these topics in a follow-up work.

## References

  * (1) J. I. Cirac and F. Verstraete, Renormalization and tensor product states in spin chains and lattices, J. Phys. A: Math. Theor. 42, 504004 (2009).
  * (2) G. Evenbly and G. Vidal, Tensor network states and geometry, J. Stat. Phys. 145, 891-918 (2011).
  * (3) R. Orus, A practical introduction to tensor networks: Matrix product states and projected entangled pair states, Ann. Phys. 349, 117 (2014).
  * (4) J. C. Bridgeman and C. T. Chubb, Hand-waving and Interpretive Dance: An Introductory Course on Tensor Networks, J. Phys. A: Math. Theor. 50, 223001 (2017).
  * (5) S. Montangero, Introduction to Tensor Network Methods - Numerical Simulations of Low-dimensional Many-body Quantum Systems, (Springer, Berlin, 2018).
  * (6) R. Orus, Tensor networks for complex quantum systems, Nat. Rev. Phys. 1, 538-550 (2019).
  * (7) P. Silvi, F. Tschirsich, M. Gerster, J. Jünemann, D. Jaschke, M. Rizzi, and S. Montangero, The Tensor Networks Anthology: Simulation techniques for many-body quantum lattice systems, SciPost Phys. Lect. Notes 8 (2019).
  * (8) S.-J. Ran, E. Tirrito, C. Peng, X. Chen, L. Tagliacozzo, G. Su, and M. Lewenstein, Tensor Network Contractions Methods and Applications to Quantum Many-Body Systems, (Springer, 2020).
  * (9) J. I. Cirac, D. Pérez-García, N. Schuch, and F. Verstraete, Matrix product states and projected entangled pair states: Concepts, symmetries, theorems, Rev. Mod. Phys. 93, 045003 (2021).
  * (10) T.G. Kolda and B.W. Bader, Tensor decompositions and applications, SIAM review, 51(3), pp.455-500 (2009).
  * (11) G. Vidal, J. I. Latorre, E. Rico, and A. Kitaev, Entanglement in quantum critical phenomena, Phys. Rev. Lett. 90, 227902 (2003).
  * (12) M B Hastings, An area law for one-dimensional quantum systems, J. Stat. Mech. P08024 (2007)
  * (13) J. Eisert, M. Cramer, and M. Plenio, Area laws for the entanglement entropy - a review, Rev. Mod. Phys. 82, pp. 277-306 (2010).
  * (14) G. K. L. Chan and S. Sharma, The density matrix renormalization group in quantum chemistry, Annu. Rev. Phys. Chem. 62, 465 (2011).
  * (15) S. Keller, M. Dolfi, M. Troyer, and M. Reiher, An efficient matrix product operator representation of the quantum chemical Hamiltonian, J. Chem. Phys. 143, 244118 (2015).
  * (16) S. Szalay, M. Pfeffer, V. Murg, G. Barcza, F. Verstraete, R. Schneider, and Ö. Legeza, Tensor product methods and entanglement optimization for ab initio quantum chemistry, Int. J. Quantum Chem. 115, 1342 (2015).
  * (17) G. K.-L. Chan, A. Keselman, N. Nakatani, Z. Li, and S. R. White, Matrix product operators, matrix product states, and ab initio density matrix renormalization group algorithms, J. Chem. Phys. 145, 014102 (2016).
  * (18) H. Zhai and G. K.-L. Chan, Low communication high performance ab initio density matrix renormalization group algorithms, J. Chem. Phys. 154, 224116 (2021).
  * (19) B. Swingle, Entanglement renormalization and holography, Phys. Rev. D 86, 065007 (2012).
  * (20) M. Miyaji, T. Numasawa, N. Shiba, T. Takayanagi, and K. Watanabe, Continuous Multiscale Entanglement Renormalization Ansatz as Holographic Surface-State Correspondence, Phys. Rev. Lett. 115, 171602 (2015).
  * (21) F. Pastawski, B. Yoshida, D. Harlow, and J. Preskill, Holographic quantum error-correcting codes: Toy models for the bulk/boundary correspondence, JHEP 06 (2015) 149.
  * (22) P. Hayden, S. Nezami, X.-L. Qi, N. Thomas, M. Walter, and Z. Yang, Holographic duality from random tensor networks, JHEP 11 (2016) 009.
  * (23) B. Czech, L. Lamprou, S. McCandlish, and J. Sully, Tensor networks from kinematic space, JHEP 07 (2016) 100.
  * (24) G. Evenbly, Hyperinvariant Tensor Networks and Holography, Phys. Rev. Lett. 119, 141602 (2017).
  * (25) E. M. Stoudenmire and D. J. Schwab, Supervised learning with tensor networks, Adv. Neural Inf. Process. Syst. 29, pp. 4799–4807 (2016).
  * (26) J. Martyn, G. Vidal, C. Roberts, and S. Leichenauer, Entanglement and tensor networks for supervised image classification, arXiv:2007.06082 (2020).
  * (27) S. Cheng, L. Wang, and P. Zhang, Supervised learning with projected entangled pair states, Phys. Rev. B 103, 125117 (2021).
  * (28) J. Liu, S. Li, J. Zhang, and P. Zhang, Tensor networks for unsupervised machine learning, arXiv:2106.12974 (2021).
  * (29) Y. Liu, W-J Li, X. Zhang, M. Lewenstein, G. Su, and S.J Ran, Entanglement-Based Feature Extraction by Tensor Network Machine Learning, Front. Appl. Math. Stat. 7 (2021).
  * (30) E. S. Fried, N. P. D. Sawaya, Y. Cao, I. D. Kivlichan, J. Romero, and A. Aspuru-Guzik, qTorch: The quantum tensor contraction handler, PLoS ONE 13(12): e0208510. (2018).
  * (31) B. Villalonga, S. Boixo, B. Nelson, C. Henze, E. Rieffel, R. Biswas and S. Mandrà, A flexible high-performance simulator for verifying and benchmarking quantum circuits implemented on real hardware, npj Quantum Inf 5, 86 (2019).
  * (32) R. Schutski, D. Lykov, and I. Oseledets, Adaptive algorithm for quantum circuit simulation, Phys. Rev. A 101, 042335 (2020).
  * (33) F. Pan and P. Zhang, Simulation of Quantum Circuits Using the Big-Batch Tensor Network Method, Phys. Rev. Lett. 128, 030501 (2022)
  * (34) M. Levental, Tensor Networks for Simulating Quantum Circuits on FPGAs, arXiv:2108.06831 (2021).
  * (35) T. Vincent, L. J. O’Riordan, M. Andrenkov, J. Brown, N. Killoran, H. Qi, and I. Dhand, Jet: Fast quantum circuit simulations with parallel task-based tensor-network contraction, arXiv:2107.09793 (2021).
  * (36) G. Evenbly and G. Vidal, Algorithms for entanglement renormalization, Phys. Rev. B 79, 144108 (2009).
  * (37) H.-H. Zhao, Z.-Y. Xie, Q.-N. Chen, Z.-C. Wei, J. W. Cai, and T. Xiang, Renormalization of tensor-network states, Phys. Rev. B 81, 174411 (2010).
  * (38) U. Schollwoeck, The density-matrix renormalization group in the age of matrix product states, Ann. Phys. 326, 96 (2011).
  * (39) H. N. Phien, J. A. Bengua, H. D. Tuan, P. Corboz, and R. Orus, The iPEPS algorithm, improved: fast full update and gauge fixing, Phys. Rev. B 92, 035142 (2015).
  * (40) G. Evenbly, Algorithms for tensor network renormalization, Phys. Rev. B 95, 045117 (2017).
  * (41) M. Fishman, E. M. Stoudenmire, and S. R. White, ITensor Library, http://itensor.org.
  * (42) Y.-J. Kao, Y.-D. Hsieh, and P. Chen, Uni10: an open-source library for tensor network algorithms, https://uni10.gitlab.io/.
  * (43) J. Haegeman, TensorOperations, https://github.com/Jutho/TensorOperations.jl.
  * (44) J. Hauschild and F. Pollmann, Efficient numerical simulations with Tensor Networks: Tensor Network Python (TeNPy), SciPost Phys. Lect. Notes 5 (2018).
  * (45) S. Al-Assam, S. R. Clark, and D. Jaksch, The Tensor Network Theory Library, J. Stat. Mech. 093102 (2017).
  * (46) R. Olivares-Amaya, W. Hu, N. Nakatani, S. Sharma, J. Yang, and G. K.-L. Chan, The ab-initio density matrix renormalization group in practice, J. Chem. Phys. 142, 034102 (2015). Accompanying software: BLOCK: Density Matrix Renormalization Group Algorithms for Quantum Chemistry, https://sanshar.github.io/Block/
  * (47) C. Roberts, A. Milsted, M. Ganahl, A. Zalcman, B. Fontaine, Y. Zou, J. Hidary, G. Vidal, and S. Leichenauer, TensorNetwork: A Library for Physics and Machine Learning, arXiv:1905.01330 (2019).
  * (48) I. V. Oseledets, TT Toolbox, https://github.com/oseledets/TT-Toolbox
  * (49) S. R. White, Density matrix formulation for quantum renormalization groups, Phys. Rev. Lett. 69, 2863 (1992).
  * (50) S. R. White, Density-matrix algorithms for quantum renormalization groups, Phys. Rev. B 48, 10345 (1993).
  * (51) U. Schollwoeck, The density-matrix renormalization group, Rev. Mod. Phys. 77, 259 (2005).
  * (52) G. Vidal, Efficient Classical Simulation of Slightly Entangled Quantum Computations, Phys. Rev. Lett. 91, 147902 (2003).
  * (53) G. Vidal, Efficient Simulation of One-Dimensional Quantum Many-Body Systems, Phys. Rev. Lett. 93, 040502 (2004).
  * (54) F. Verstraete and J. I. Cirac, Renormalization algorithms for quantum-many-body systems in two and higher dimensions, arXiv:cond-mat/0407066.
  * (55) F. Verstraete, J.I. Cirac, and V. Murg, Matrix product states, projected entangled pair states, and variational renormalization group methods for quantum spin systems, Adv. Phys. 57, 143 (2008).
  * (56) J. Jordan, R. Orus, G. Vidal, F. Verstraete, and J. I. Cirac, Classical simulation of infinite-size quantum lattice systems in two spatial dimensions, Phys. Rev. Lett. 101, 250602 (2008).
  * (57) G. Vidal, A class of quantum many-body states that can be efficiently simulated, Phys. Rev. Lett. 101, 110501 (2008).
  * (58) M. Levin and C. P. Nave, Tensor renormalization group approach to two-dimensional classical lattice models, Phys. Rev. Lett. 99, 120601 (2007).
  * (59) Z.-Y. Xie, J. Chen, M. P. Qin, J. W. Zhu, L. P. Yang, and T. Xiang, Coarse-graining renormalization by higher-order singular value decomposition, Phys. Rev. B 86, 045139 (2012).
  * (60) G. Evenbly and G. Vidal, Tensor network renormalization, Phys. Rev. Lett. 115, 180405 (2015).
  * (61) G. Evenbly, Tensors.net website, https://www.tensors.net.
  * (62) S. Singh, R. N. C. Pfeifer, and G. Vidal, Tensor network decompositions in the presence of a global symmetry, Phys. Rev. A 82, 050301 (2010).
  * (63) S. Singh, R. N. C. Pfeifer, and G. Vidal, Tensor network states and algorithms in the presence of a global U(1) symmetry, Phys. Rev. B 83, 115125 (2011)
  * (64) A. Weichselbaum, Non-abelian symmetries in tensor networks: a quantum symmetry space approach, Ann. Phys. 327, 2972-3047 (2012).
  * (65) S. Sharma, A general non-Abelian density matrix renormalization group algorithm with application to the $C_{2}$ dimer, J. Chem. Phys. 142, 024107 (2015).
  * (66) S. Keller and M. Reiher, Spin-adapted matrix product states and operators, J. Chem. Phys. 144, 134101 (2016).
  * (67) P. Nataf and F. Mila, Density matrix renormalization group simulations of SU(N) Heisenberg chains using standard Young tableaus: Fundamental representation and comparison with a finite-size Bethe ansatz, Phys. Rev. B 97, 134420 (2018).
  * (68) P. Schmoll, S. Singh, M. Rizzi, and R. Orús, A programming guide for tensor networks with global SU(2) symmetry Ann. Phys. 419, 168232 (2020).
  * (69) R. N. C. Pfeifer, J. Haegeman, and F. Verstraete, Faster identification of optimal contraction sequences for tensor networks, Phys. Rev. E 90, 033315 (2014).
  * (70) R. N. C. Pfeifer and G. Evenbly, Improving the efficiency of variational tensor network algorithms, Phys. Rev. B 89, 245118 (2014).
  * (71) J. M. Dudek, L. Dueñas-Osorio, and M. Y. Vardi, Efficient Contraction of Large Tensor Networks for Weighted Model Counting through Graph Decompositions, arXiv:1908.04381v2 (2019).
  * (72) J. Gray and S. Kourtis, Hyper-optimized tensor network contraction, Quantum 5, 410 (2021).
  * (73) R. N. C. Pfeifer, G. Evenbly, S. Singh, and G. Vidal, NCON: A tensor network contractor for MATLAB, arXiv:1402.0939 (2014).
  * (74) R. A. Horn, C. R. Johnson, Matrix analysis, (Cambridge University Press, 1985).
  * (75) R. A. Horn, C. R. Johnson, Topics in matrix analysis, (Cambridge University Press 1991).
  * (76) C. Eckart and G. Young, The approximation of one matrix by another of lower rank, Psychometrika 1, 211–218 (1936).
  * (77) Y. Shi, L. Duan and G. Vidal, Classical simulation of quantum many-body systems with a tree tensor network, Phys. Rev. A 74, 022320 (2006).
  * (78) L. Tagliacozzo, G. Evenbly and G. Vidal, Simulation of two-dimensional quantum systems using a tree tensor network that exploits the entropic area law, Phys. Rev. B 80, 235127 (2009).
  * (79) G. Evenbly, Gauge fixing, canonical forms and optimal truncations in tensor networks with closed loops, Phys. Rev. B 98, 085155 (2018).
  * (80) S. Holtz, T. Rohwedder, and R. Schneider, The Alternating Linear Scheme for Tensor Optimization in the Tensor Train Format, SIAM J. Sci. Comput., 34(2), A683–A713.
  * (81) Y. Zhang and E. Solomonik, On Stability of Tensor Networks and Canonical Forms, arXiv:2001.01191 (2020).


Compiled

# Airy beams on incoherent background

Morteza Hajati  Department of Electrical and Computer Engineering, Dalhousie
University, Halifax, Nova Scotia, B3J 2X4, Canada  Vincent Sieben  Department
of Electrical and Computer Engineering, Dalhousie University, Halifax, Nova
Scotia, B3J 2X4, Canada  Sergey A. Ponomarenko  Department of Electrical and
Computer Engineering, Dalhousie University, Halifax, Nova Scotia, B3J 2X4,
Canada  Department of Physics and Atmospheric Science, Dalhousie University,
Halifax, Nova Scotia, B3H 4R2, Canada<EMAIL_ADDRESS>
###### Abstract

We present a class of diffraction-free partially coherent beams each member of
which is comprised of a finite-power, non-accelerating Airy bump residing on a
statistically homogeneous, Gaussian-correlated background. We examine free-
space propagation of soft apertured realizations of the proposed beams and
show that their evolution is governed by two spatial scales: the coherence
width of the background and aperture size. A relative magnitude of these
factors determines the practical range of propagation distances over which the
novel beams can withstand diffraction. The proposed beams can find
applications to imaging and optical communications through random media.

###### doi:

http://dx.doi.org/10.1364/OL.XX.XXXXXX

††journal: ol

Partially coherent light sources have attracted a great deal of interest over
the years because they generate a plethora of structured statistical light
fields with intriguing properties [1, 2] that carry promise for multiple
applications [3]. The majority of partially coherent optical fields studied to
date, though, have been statistically homogeneous [4]. Yet, statistically
nonuniform fields have also been explored, both theoretically [5, 6, 7, 8, 9]
and experimentally [11, 10, 12, 13]. In particular, propagation invariant
random fields have been shown theoretically [6, 7] and verified experimentally
[13] to necessarily possess nonuniform correlations as they manifest
themselves as bumps or dips residing atop of a statistically homogeneous
background. Such statistical diffraction free beams have been shown to possess
universal self-similar asymptotic propagation properties in random media [14].
Several generalizations of random diffraction free beams were proposed [15,
16, 17], including, notably, white-light (polychromatic) dark/antidark fields
on incoherent background [15].

At the same time, coherent diffraction free Airy waves in general [18] and
optical Airy beams in particular [19] have attracted much attention since
their experimental realization [20], not the least bit because of their
intriguing accelerating nature and a wealth of potential applications [21].
Partially coherent Airy beams have also been realized [22]. However, the
latter, albeit accelerating, are not genuine diffraction free beams, even in
the ideal, infinite-power limit since their cross-spectral density is not of
the form of a bump (or dip) on a statistically uniform background. This
observation begs a fundamental question: Can one generate bona fide
diffraction free random Airy beams, and if so, will their finite-power
realizations accelerate in free space?

In this Letter, we demonstrate how a family of genuine diffraction free,
partially coherent Airy sources can be constructed. We show that any member of
the family is represented by an Airy bump situated on a Gaussian correlated
uniform background. The coherence properties of any beam generated by such a
source are controlled by a single dimensionless parameter, a correlation
parameter given by the ratio of a coherence width of the background to a
characteristic width of the bump. We study free-space propagation of Gaussian
apertured realizations of the proposed beams and show that while their
acceleration is suppressed, the distance over which the new beams can
withstand diffraction is determined by the relative magnitude of two spatial
scales: the background coherence width and aperture size.

We start with a general representation for the cross-spectral density of any
nonnegative definite 1D partially coherent source in the form [23]

$W(x_{1},x_{2})=\int_{-\infty}^{\infty}dk\,p(k)\mathcal{A}^{\ast}(k,x_{1})\mathcal{A}(k,x_{2}),$
(1)

where we dub the nonnegative real function $p(k)\geq 0$, $p^{\ast}(k)=p(k)$, a
spectral distribution and $\\{\mathcal{A}(k,x)\\}$ generalized modes of the
source corresponding to a continuous spectrum labelled by the variable $k$.

Let us consider a Gaussian spectral distribution of the modes,

$p(k)=\frac{\sigma_{c}}{\sqrt{2\pi}}e^{-k^{2}\sigma_{c}^{2}/2},$ (2)

where $\sigma_{c}$ is the source coherence length, and the plane-wave modes
with cubic phase chirp, so that

$\mathcal{A}(k,x)=\sqrt{2I_{0}}\cos\left(\textstyle\frac{1}{6}\sigma_{I}^{3}k^{3}+\textstyle\frac{1}{2}kx\right).$
(3)

$I_{0}$ and $\sigma_{I}$, entering Eq. (3), set the scales of the source
intensity and spatial width of the bump, respectively. Numerical coefficients
are chosen for convenience. In physics terms, Eqs. (1) through (3) imply that
the source is composed of a continuum of phase chirped plane waves, with each
plane wave contribution to the overall source intensity being weighed with the
Gaussian in Eq. (2). We can also surmise that $\sigma_{c}$ controls the source
state of coherence. In particular, in the limit $\sigma_{c}\rightarrow\infty$,
we can infer at once that $p(k)\rightarrow\delta(k)$, implying the source
engenders a chirped plane wave in the fully coherent limit. As we will see in
more precise terms shortly, $\sigma_{c}$ plays the role of the background
coherence width.

It follows from the trigonometric identity,

$\cos\alpha\cos\beta=\textstyle\frac{1}{2}[\cos(\alpha-\beta)+\cos(\alpha+\beta)],$
(4)

that

$\mathcal{A}^{\ast}(k,x_{1})\mathcal{A}(k,x_{2})=I_{0}[\cos(kx_{-}/2)+\cos(\sigma_{I}^{3}k^{3}/3+kx_{+})],$
(5)

where we introduced the center-of-mass and difference coordinates as

$x_{+}=(x_{1}+x_{2})/2,\hskip 14.22636ptx_{-}=x_{1}-x_{2}.$ (6)

On substituting from Eq. (5) into Eq. (1), performing Gaussian integrations
and making use of the integral representation of the Airy function,

$\mathrm{Ai}(x/\sigma_{I})=\frac{\sigma_{I}}{2\pi}\int_{-\infty}^{\infty}ds\,\exp\left(\frac{i\sigma_{I}^{3}s^{3}}{3}+isx\right),$
(7)

we arrive at the cross-spectral density of an Airy beam on incoherent
background (ABIB) in the form

$\overline{W}_{\mathrm{A}}(X_{-},X_{+})=\overline{W}_{\mathrm{bg}}(X_{-})+\overline{W}_{\mathrm{bp}}(X_{+}),$
(8)

where

$\overline{W}_{\mathrm{bg}}(X_{-})=\exp\left(-\frac{X_{-}^{2}}{8\xi_{c}}\right),$
(9)

is the cross-spectral density of a statistically uniform, Gaussian-correlated
background and

$\overline{W}_{\mathrm{bp}}(X_{+})=P_{\mathrm{1D}}e^{\xi_{c}^{3}/12}e^{\xi_{c}X_{+}/2}\mathrm{Ai}(X_{+}-X_{0})$
(10)

is that of an Airy bump atop of the background. Here we introduced
dimensionless variables, $\overline{W}=W/I_{0}$ and
$X_{\pm}=x_{\pm}/\sigma_{I}$. Further, we introduced a correlation parameter
$\xi_{c}=\sigma_{c}^{2}/\sigma_{I}^{2}$ as well as a coordinate
$X_{0}=-\xi_{c}^{2}/4$ of the Airy function maximum and a total power
$P_{\mathrm{1D}}=\sqrt{2\pi\xi_{c}}$ of the Airy bump. The source intensity
$\overline{I}_{\mathrm{A}}(X)=\overline{W}_{\mathrm{A}}(X,X)$, can be
determined at once from Eq. (8) by setting $X_{-}=0$ and $X_{+}=X$ which
yields the expression

$\overline{I}_{\mathrm{A}}(X)=1+P_{\mathrm{1D}}e^{\xi_{c}^{3}/12}e^{\xi_{c}X/2}\mathrm{Ai}(X-X_{0}).$
(11)

At this point, let us make some instructive observations. First, Eqs. (8)
through (11) imply that the newly constructed sources generate bona fide
diffraction-free beams as the source cross-spectral density has the
propagation invariant form that was first explicitly derived in [6]. Second,
the properties of a family of ABIBs are governed—in dimensionless units—by the
single correlation parameter $\xi_{c}$. Third, ideal Airy beams on incoherent
background do not accelerate: the peak intensity of the bump is located at a
fixed position dependent on the background coherence. Fourth, each Airy bump
carries a finite power $P_{\mathrm{1D}}$ which increases with the coherence
level of the background. The last two properties of the discovered beams bring
into focus sharp contrast between ABIBs and fully coherent Airy beams.
Interestingly, the exponential cut-off of the source intensity at negative
coordinates arises naturally in the partially coherent case.

Figure 1: (a) Intensity profile of an ideal ABIB as a function of $X$ for two
values of the correlation parameter $\xi_{c}$: $\xi_{c}=0.1$ (solid curve) and
$\xi_{c}=2$ (dashed curve). (b) Contour plot of the ABIB intensity.

To visualize the spatial structure of a typical ABIB, we display the intensity
distribution of the beam for two values of the correlation parameter in Fig.
1. We can infer from the figure that while a rather correlated ABIB, such that
$\sigma_{I}\leq\sigma_{c}$ corresponding to $\xi_{c}\geq 1$, has a relatively
high-intensity bump, a nearly uncorrelated ABIB with $\sigma_{c}\ll\sigma_{I}$
features a short bump and a long oscillatory left tail. This behavior can be
explained by noticing that a weakly correlated background breaks up the beam
into a number of essentially uncorrelated beamlets, each having nearly the
same peak intensity as its neighbor. While the oscillations (beamlets) in the
left tail of the ABIB are gently cut-off at long $X$ by the exponential factor
$e^{\xi_{c}X/2}$, the ones in the right tail are sharply suppressed by the
super-exponential asymptotics of the Airy function, $\mathrm{Ai}(X)\sim
e^{-(2/3)X^{3/2}}$ as $X\rightarrow+\infty$.

Figure 2: Illustrating the hierarchy of transverse scales associated with a
Gaussian apertured Airy beam on incoherent background.

An ideal ABIB carries infinite power because it dwells on a statistically
uniform background, and hence is not realizable in the laboratory. We are then
bound to consider a “softly" apertured ABIB, which can be experimentally
realized with the aid of a Gaussian amplitude mask, for example, and examine
free-space paraxial propagation of the Gaussian apertured ABIB. We assume that
the aperture width $w_{0}$ dwarfs the transverse width of the bump,
$\sigma_{I}\ll w_{0}$. We schematically illustrate the hierarchy of transverse
scales governing free-space propagation of finite-power ABIBs in Fig. 2. The
cross-spectral density of the apertured ABIB source has the form

$W^{(0)}(x_{1},x_{2})=\exp\left(-\frac{x_{1}^{2}+x_{2}^{2}}{2w_{0}^{2}}\right)\,W_{\mathrm{A}}(x_{1},x_{2}),$
(12)

where $W_{\mathrm{A}}$ is the cross-spectral density of the ideal ABIB given
by Eqs. (8) through (10). The intensity of the apertured ABIB with the carrier
wavenumber $k_{0}$ in any cross-section $z\geq 0$ is given by a Fresnel
transform [4]

$\displaystyle I(x,z)$ $\displaystyle=\left(\frac{k_{0}}{2\pi z}\right)\int
dx_{1}\int dx_{2}\,W^{(0)}(x_{1},x_{2})e^{ik_{0}(x-x_{2})^{2}/2z}$
$\displaystyle\times e^{-ik_{0}(x-x_{1})^{2}/2z}.$ (13)

Figure 3: (a) and (c) Intensity profiles of a nearly uncorrelated Gaussian
apertured ABIB as a function of $X$ for three propagation distances:
$z=0.1z_{c}$ (solid curve), $z=z_{c}$ (dashed curve), and $z=3z_{c}$ (dash-
dotted curve) and the aperture size (a) $\epsilon=0.01$ and (c)
$\epsilon=0.1$. (b) and (d) Contour plots of the Gaussian apertured ABIB with
(b) $\epsilon=0.01$ and (d) $\epsilon=0.1$. The correlation parameter is
$\xi_{c}=0.1$.

On substituting from Eqs. (8) through (10) as well as from Eq. (12) into (Airy
beams on incoherent background) and converting to dimensionless center-of-mass
and difference coordinates, we obtain, after quite tedious but straightforward
integrations involving Gaussians and the Airy function representation, the
following expression

$\overline{I}(X,z)=\overline{I}_{\mathrm{bg}}(X,z)+\overline{I}_{\mathrm{bp}}(X,z),$
(14)

with

$\overline{I}_{\mathrm{bg}}(X,z)=\frac{1}{\sqrt{1+z^{2}/z_{c}^{2}}}\,\exp\left(-\frac{\epsilon^{2}X^{2}}{1+z^{2}/z_{c}^{2}}\right),$
(15)

and

$\displaystyle\overline{I}_{\mathrm{bp}}(X,z)$
$\displaystyle=\overline{I}_{0}(z)\exp\left[\frac{\xi_{c}X}{2}\frac{(1+z^{2}/z_{c}^{2})}{(1+z^{2}/z_{R}^{2})^{2}}\right]$
$\displaystyle\times\exp\left(-\frac{\epsilon^{2}X^{2}}{1+z^{2}/z_{R}^{2}}\right)\mathrm{Ai}\left[\frac{X-X_{0}(z)}{1+z^{2}/z_{R}^{2}}\right].$
(16)

In Eqs. (14), (15) and (Airy beams on incoherent background) we introduced a
small dimensionless parameter $\epsilon=\sigma_{I}/w_{0}\ll 1$ and

$\overline{I}_{0}(z)=\frac{P_{\mathrm{1D}}}{\sqrt{1+z^{2}/z_{R}^{2}}}\exp\left[\frac{\xi_{c}^{3}}{12}\left(\frac{1+z^{2}/z_{c}^{2}}{1+z^{2}/z_{R}^{2}}\right)^{3}\right],$
(17)

as well as

$X_{0}(z)=-\frac{\xi_{c}^{2}}{4}\frac{(1+z^{2}/z_{c}^{2})^{2}}{1+z^{2}/z_{R}^{2}}.$
(18)

The physics of apertured ABIB evolution in free space is incapsulated in two
longitudinal scales: a Rayleigh range $z_{R}=k_{0}w_{0}^{2}$ of the embedding
Gaussian beam as well as a characteristic diffraction length
$z_{c}=\sqrt{2}k_{0}\sigma_{c}w_{0}$ associated with finite background
coherence. We can infer from Eq. (15) that the Gaussian background spreads, so
that its coherence state changes appreciably at the distances of the order of
$z_{c}$. By the same token, it follows from Eqs. (Airy beams on incoherent
background), (17) and (18) that while the bump maintains its Airy shape, it
spreads slowly over a characteristic distance of the order of $z_{R}$. It can
then be viewed as practically diffraction-free over distances $z\ll z_{R}$. A
specific evolution scenario of the ABIB depends on whether the background is
nearly incoherent, $\sigma_{c}\ll\sigma_{I}$, or highly coherent,
$\sigma_{I}\ll\sigma_{c}$. In the former case, $z_{c}$ sets the scale for
diffraction, while in the latter, a characteristic distance over which the
ABIB succumbs to diffraction is determined by the Rayleigh range of the
Gaussian aperture.

To illustrate our qualitative conclusions, we exhibit the evolution of an
apertured Airy bump on a nearly incoherent background (weakly correlated ABIB)
in Fig. 3 for three propagation distances measured in the units of $z_{c}$
which sets the scale for diffraction there. We consider two cases:
$\epsilon=0.01$ and $\epsilon=0.1$ corresponding to a very large and
moderately sized apertures, respectively. In the wide-aperture limit, see Fig
3(a), the ABIB remains nearly diffraction-free. The beam maintains its ideal
structure with a long oscillatory tail, although the peak intensities of the
main and secondary lobes diminish with the propagation distance, as does the
magnitude of the background intensity. As the aperture size is reduced,
though, the ABIB structure undergoes interesting evolution. The beam tails,
which are partially cut off by the aperture at short propagation distances
from the source plane, reappear on propagation farther away, as the ABIB
gradually transfers the power from its main lobe into the tails. This scenario
is evident in Fig. 3 (c). Eventually, over long enough propagation distances,
the side lobes virtually disappear and the tightly apertured ABIB profile
becomes nearly identical to that generated by diffraction of a loosely
apertured ABIB. This conclusion can be confirmed by comparing the dash-dotted
curves corresponding to the propagation distance $z=3z_{c}$ in panels (a) and
(c) of Fig. 3. Thus, a nearly universal shape of the weakly correlated ABIB
emerges over sufficiently long propagation distances, regardless of the size
of the Gaussian aperture. The magnitude of the propagation distance over which
the reshaping toward the asymptotic profile occurs does depend on the aperture
size through the dependence of $z_{c}$ on $w_{0}$, of course. At the same
time, embedding an Airy bump on a fairly coherent background (strongly
correlated ABIB) into a Gaussian envelope results in the ABIB spreading over a
few characteristic Rayleigh distances as can be seen in Fig. 4. A strongly
correlated ABIB diffracts much like an ordinary Gaussian beam which makes it
less attractive for potential applications than its weakly correlated
counterpart. We also notice by inspecting the contour plots in Figs. 4(b) and
3(b,d) that the ABIB acceleration is strongly suppressed in the partially
coherent case, irrespective of the aperture size and coherence level of the
background.

Figure 4: (a) Intensity profiles of a strongly correlated Gaussian apertured
ABIB as a function of $X$ for three propagation distances: $z=0.1z_{R}$ (solid
curve), $z=z_{R}$ (dashed curve), and $z=3z_{R}$ (dash-dotted curve) and the
aperture size $\epsilon=0.3$. (b) Contour plot of the Gaussian apertured ABIB.
The correlation parameter is $\xi_{c}=10$.

Finally, a generalization to a 2D ABIB is straightforward. To this end, we can
represent the cross-spectral density of the 2D ABIB as

$W({\bf r}_{1},{\bf r}_{2})=\int d{\bf k}\,p({\bf k})\mathcal{A}^{\ast}({\bf
k},{\bf r}_{1})\mathcal{A}({\bf k},{\bf r}_{2}),$ (19)

where

$p({\bf
k})=\prod_{j=x,y}\frac{\sigma_{cj}}{\sqrt{2\pi}}e^{-k_{j}^{2}\sigma_{cj}^{2}/2},$
(20)

and

$\mathcal{A}({\bf k},{\bf
r})=\sqrt{2I_{0}}\cos\left(\textstyle\frac{1}{6}\sum_{j=x,y}\sigma_{Ij}^{3}k_{j}^{3}+\textstyle\frac{1}{2}{\bf
k}\cdot{\bf r}\right).$ (21)

Figure 5: Contour plots of ideal 2D ABIB corresponding to (a) $\xi_{c}=0.01$,
(b) $\xi_{c}=0.1$, (c) $\xi_{c}=1$ and (d) $\xi_{c}=3$.

Here we introduced the background coherence $\sigma_{cj}$ and bump intensity
$\sigma_{Ij}$ widths along the corresponding axes $j=x,y$ of the Cartesian
coordinate system; the source is, in general, anisotropic. Repeating exactly
the same steps as we took to derive the cross-spectral density of the 1D ABIB,
we arrive at the expression

$\overline{W}_{\mathrm{A}}({\bf R}_{-},{\bf
R}_{+})=\overline{W}_{\mathrm{bg}}({\bf
R}_{-})+\overline{W}_{\mathrm{bp}}({\bf R}_{+}),$ (22)

with

$\overline{W}_{\mathrm{bg}}({\bf
R}_{-})=\exp\left(-\frac{X_{-}^{2}}{8\xi_{cx}}-\frac{Y_{-}^{2}}{8\xi_{cy}}\right),$
(23)

and

$\displaystyle\overline{W}_{\mathrm{bp}}({\bf R}_{+})$
$\displaystyle=P_{\mathrm{2D}}\exp\left[(\xi_{cx}^{3}+\xi_{cy}^{3})/12\right]e^{\xi_{cx}X_{+}/2}e^{\xi_{cy}Y_{+}/2}$
$\displaystyle\mbox{}\times\mathrm{Ai}(X_{+}-X_{0})\mathrm{Ai}(Y_{+}-Y_{0}),$
(24)

where we introduced correlation parameters
$\xi_{cj}=\sigma_{cj}^{2}/\sigma_{Ij}^{2}$ and a total power of a 2D Airy bump
as $P_{\mathrm{2D}}=2\pi\sqrt{\xi_{cx}\xi_{cy}}$. Notice that although the
cross-spectral densities of the bump and background individually factorize in
Cartesian coordinates, the cross-spectral density of a 2D Airy beam on
incoherent background does not. We exhibit the intensity profiles of ideal 2D
ABIBs in Fig. 5. We can infer from the figure by comparing panels (a) through
(d) that as the source becomes progressively more correlated—or the background
coherence level increases—the ABIB power localizes more within the main lobe.
This picture is consistent with our findings for 1D ABIBs.

In summary, we have introduced a class of bona fide diffraction free
statistical beams, Airy beams on incoherent background. While ideal Airy beams
on incoherent background do not accelerate in free space, their soft-apertured
realizations, which carry finite power and hence are amenable to the
laboratory implementation, have their acceleration strongly suppressed. The
actual distance over which any finite-power ABIB defies diffraction is
determined by the interplay of the coherence width of the random background
and the aperture size. We anticipate the proposed beams to find applications
to imaging as well as to optical communications and information transfer
through random media.

Funding. Natural Sciences and Engineering Research Council of Canada
(RGPIN-2018-05497).

Disclosures. The authors declare no conflicts of interest.

## References

  * [1] Y. Cai, Y. Chen, J. Yu, X. Liu, and L. Liu, Prog. Opt. 62, 157 (2017).
  * [2] Y. Chen, A. Norrman, S. A. Ponomarenko, and A. T. Friberg, Prog. Opt. 65, 105 (2020).
  * [3] O. Korotkova and G. Gbur, Prog. Opt. 65, 43 (2020).
  * [4] L. Mandel and E. Wolf, _Optical Coherence and Quantum Optics_ (Cambridge University, 1995).
  * [5] S. A. Ponomarenko, J. Opt. Soc. Am. A, 18, 150 (2001).
  * [6] S. A. Ponomarenko, W. Huang, and M. Cada, Opt. Lett. 32, 2508 (2007).
  * [7] R. Borghi, F. Gori, and S. A. Ponomarenko, J. Opt. Soc. Am. A, 26, 2275 (2009).
  * [8] H. Lajunen and T. Saastamoinen, Opt. Lett. 36, 4104 (2011).
  * [9] Y. Chen , A. Norrman, S. A. Ponomarenko, and A. T. Friberg, Phys. Rev. A, 100, 053833 (2019).
  * [10] G. V. Bogatyryova, C. V. Felde, P. V. Polyanskii, S. A. Ponomarenko, M. S. Soskin, and E.Wolf, Opt. Lett. 28, 878 (2003).
  * [11] A. S. Ostrovsky, J. Garcia-Garcia, C. Rickenstorff-Parrao, and M. A. Olvera-Santamaria, Opt. Lett. 42, 5182 (2017).
  * [12] X. Chen, J. Li, S. M. H. Rafsanjani, and O. Korotkova, Opt. Lett. 43, 3590 (2018).
  * [13] X. Zhu, F. Wang, C. Zhao,Y. Cai, and S. A. Ponomarenko, Opt. Lett., 44, 2260 (2019).
  * [14] X. Li, S. A. Ponomarenko, Z. Xu, F. Wang, Y. Cai, and C. Liang, Opt. Lett., 45, 698 (2020)
  * [15] M. Yessenov, B. Bhaduri, H. E. Kondakci, M. Meem, R. Menon, and A. F. Abouraddy, Optica, 6, 598 (2019).
  * [16] M. W. Hyde and S. Avramov-Zamourovich, J. Opt. Soc. Am. A, 36 , 1058 (2019).
  * [17] M. W. Hyde, J. Phys. Commun. 4 015025 (2020).
  * [18] M. V. Berry and N. L. Balazs, Am. J. Phys. 47, 264 (1979).
  * [19] G. A. Siviloglu and D. N. Christodoulides, Opt. Lett. 32, 979 (2007).
  * [20] G. A. Siviloglu, J. Broky, A. Dogariu, and D. N. Christodoulides, Phys. Rev. Lett. 99, 213901 (2007).
  * [21] N. K. Efremidis, Z. Chen, M. Segev, and D. N. Christodoulides, Optica, 6, 686 (2019).
  * [22] Y. Lumer, Y. Liang, R. Schley, I. Kaminer, E. Greenfield, D. Song, X. Zhang, J. Xu, Z. Chen, and M. Segev, Optica, 2, 886 (2015).
  * [23] F. Gori and M. Santarsiero, Opt. Lett. 32, 3531 (2007).

## References

  * [1] Y. Cai, Y. Chen, J. Yu, X. Liu, and L. Liu, “Generation of Partially Coherent Beams", Prog. Opt. 62, 157-223 (2017).
  * [2] Y. Chen, A. Norrman, S. A. Ponomarenko, and A. T. Friberg, “Optical coherence and electromagnetic surface waves", Prog. Opt. 65, 105-166 (2020).
  * [3] O. Korotkova and G. Gbur, “Applications of optical coherence theory", Prog. Opt. 65, 43-104 (2020).
  * [4] L. Mandel and E. Wolf, _Optical Coherence and Quantum Optics_ (Cambridge University, 1995).
  * [5] S. A. Ponomarenko, ‘’A class of partially coherent beams carrying optical vortices”,J. Opt. Soc. Am. A, 18, 150-156 (2001).
  * [6] S. A. Ponomarenko, W. Huang, and M. Cada, ‘’Dark and anti-dark diffraction-free beams”, Opt. Lett., 32, 2508-2510 (2007).
  * [7] R. Borghi, F. Gori, and S. A. Ponomarenko, “On a class of electromagnetic diffraction free beams", J. Opt. Soc. Am. A, 26, 2275-2281 (2009).
  * [8] H. Lajunen and T. Saastamoinen, ‘’Propagation characteristics of partially coherent beams with spatially varying correlations”,Opt. Lett. 36, 4104-4106 (2011).
  * [9] Y. Chen , A. Norrman, S. A. Ponomarenko, and A. T. Friberg, “Partially coherent surface plasmon polariton vortex fields", Phys. Rev. A, 100, 053833 (2019).
  * [10] A. S. Ostrovsky, J. Garcia-Garcia, C. Rickenstorff-Parrao, and M. A. Olvera-Santamaria,‘’Partially coherent diffraction-free vortex beams with a Bessel-mode structure”, Opt. Lett. 42, 5182-5185 (2017).
  * [11] G. V. Bogatyryova, C. V. Felde, P. V. Polyanskii, S. A. Ponomarenko, M. S. Soskin, and E.Wolf,‘’Partially coherent vortex beams with a separable phase”, Opt. Lett., 28, 878-880 (2003).
  * [12] X. Chen, J. Li, S. M. H. Rafsanjani, and O. Korotkova, ‘’Synthesis of Im-Bessel correlated beams via coherent modes”,Opt. Lett. 43, 3590-3593 (2018).
  * [13] X. Zhu, F. Wang, C. Zhao,Y. Cai, and S. A. Ponomarenko,‘’Experimental realization of dark and antidark diffraction-free beams”, Opt. Lett., 44, 2260-2263 (2019).
  * [14] X. Li, S. A. Ponomarenko, Z. Xu, F. Wang, Y. Cai, and C. Liang, “Universal self-similar asymptotic behavior of optical bump spreading in random medium atop of incoherent background", Opt. Lett., 45, 698-701 (2020)
  * [15] M. Yessenov, B. Bhaduri, H. E. Kondakci, M. Meem, R. Menon, and A. F. Abouraddy,‘’Non-diffracting broadband incoherent space–time fields”, Optica, 6, 598-607 (2019).
  * [16] M. W. Hyde and S. Avramov-Zamourovich,‘’Generating dark and antidark beams using the genuine cross-spectral density function criterion”, J. Opt. Soc. Am. A, 36 , 1058-1063 (2019).
  * [17] M. W. Hyde, “Generating electromagnetic dark and antidark partially coherent sources", J. Phys. Commun. 4 015025 (2020).
  * [18] M. V. Berry and N. L. Balazs, “Nonspreading wave packets" Am. J. Phys. 47, 264-267 (1979).
  * [19] G. A. Siviloglu and D. N. Christodoulides,“Accelerating finite-energy Airy beams", Opt. Lett. 32, 979-981 (2007).
  * [20] G. A. Siviloglu, J. Broky, A. Dogariu, and D. N. Christodoulides, “Observation of accelerating Airy beams", Phys. Rev. Lett. 99, 213901 (2007).
  * [21] N. K. Efremidis, Z. Chen, M. Segev, and D. N. Christodoulides, “Airy beams and accelerating waves: An overview of recent advances", Optica, 6, 686-701 (2019).
  * [22] Y. Lumer, Y. Liang, R. Schley, I. Kaminer, E. Greenfield, D. Song, X. Zhang, J. Xu, Z. Chen, and M. Segev, “Incoherent self-accelerating beams", Optica, 2, 886-892 (2015).
  * [23] F. Gori and M. Santarsiero, “Devising genuine spatial correlation functions", Opt. Lett. 32, 3531-3533 (2007).