In this paper, we explore the interaction between two monoidal
structures: a multiplicative one, for the encoding of pairing, and
an additive one, for the encoding of choice. We propose a PROP to
model computation in this framework, where the choice is
parametrized by an algebraic side effect: the model can support
regular tests, probabilistic and non-deterministic branching, as
well as quantum branching, i.e. superposition.

The graphical language comes equipped with a denotational semantics
based on linear applications, and an equational theory. We prove the
language to be universal, and the equational theory to be complete
with respect to the semantics.

§ INTRODUCTION

The basic execution flow of a computation is arguably based on three
notions: sequences, tuples and branches. Sequences form the building
block of compositionality, tuples are what makes it possible to
consider multiple pieces of information together, while branches allow
the behavior to change depending on the inputs or on the state of the

Ranging from (sometimes informal) flow-chart languages <cit.> to sophisticated
structures such as interaction or proof nets <cit.>, graphical
languages are commonly used to represent the possible control flow
of a computation. On a formal level, a graphical language is a
PROP <cit.>, that is, a symmetric, strict monoidal structure
$(\mathcal{C},\top,{\boxtimes})$ whose objects are of the form
$W\boxtimes\cdots\boxtimes W$. The object $W$ is a “wire”, and any
object stands for a bunch of wires. The monoidal structure formalizes
how the bunching of wires behaves.

Two Canonical Monoidal Structures.

The monoidal structure of a PROP is very versatile. On one hand, it
can be considered in a multiplicative way, with $A\boxtimes B$
seen as the pairing of an element of type $A$ and an element of type
$B$. This approach is one followed in the design of MLL proof-nets for
instance <cit.>. On the other hand, one can consider the monoidal
structure in an additive way, with $\boxtimes$ for instance
being a co- or a bi-product. Standard examples are the category FinRel
of finite sets and relations, forming an additive PROP with $\boxtimes$
being the disjoint union, or the category of finite dimensional vector
spaces (or semi-modules) and linear maps, with $\boxtimes$ being the
cartesian product.

From a computational perspective, an additive monoidal structure can
be regarded as the possibility to choose a computational path
upon the state of the input. Depending on the underlying system, this
choice can be regarded as deterministic (if based on Set),
non-deterministic (if based on Rel), probabilistic (if based on a
suitable semi-module), etc.

To be able to handle both pairing and branching in a PROP, we cannot
uniquely identify $\boxtimes$ as being multiplicative additive. We
instead need to extend the PROP with two additional monoidal
structures, one for pairing (${\otimes}$) and one for branching

Split over coproduct

Splits over coproduct and tensor

Examples of branchings

In this paper, we focus on a framework where these two monoidal
structures are available. Graphical languages for such a setting
usually rely on a notion of sheet, or worlds, to handle
branching <cit.>. Figure <ref>
shows for instance how to represent the construction of the morphism
$f\oplus g : A\oplus A'\to B\oplus B'$ out of $f:A\to A'$ and
$g:B\to B'$. The symbol “$\oplus$” stands for the “split” of
worlds. Such a graphical language therefore comes with two distinct
“splits”: one for the monoidal structure —leaving inside one
specific world—, and one for the coproduct —splitting
worlds—. They can be intertwined, as shown in
Figure <ref>. Another approach followed by
<cit.> externalizes the two products (tensor product
and coproduct) into the structure of the diagrams themselves, at the
price of a less intuitive tensor product and a form of synchronization

However, in the state of the art this “splitting-world”
understanding has only been carried for deterministic or probabilistic
branching <cit.>. These
existing approaches do not support more exotic branchings, such as
quantum superposition.

Quantum Computation.

Conventional wisdom has it that quantum computation is about
quantum data in superposition. In the standard model, the
memory holding quantum data is encapsulated inside a coprocessor
accessed through a simple interface: The coprocessor holds
individually addressable registers holding quantum bits, on
which one can apply a fixed set of operations —gates—
specified by the interface. If some of these gates can generate
superposition of data, this is kept inside the coprocessor and opaque
to the programmer. A typical interaction with the coprocessor is a
purely classical sequence of elementary operations of the form “Apply
gate X to register $n$; apply gate Y to register $m$;
etc”. Such a sequence of instructions is usually represented
as a quantum circuit. In this model, a quantum program is then
a conventional program building a quantum circuit and sending it as a
batch-job to the coprocessor.

From a semantical perspective, the state of a quantum memory
consisting of $n$ quantum bits is a vector in a $2^n$-dimensional
Hilbert space. A (pure) quantum circuit is a linear, sequential description
of elementary operations describing a linear, unitary map on
the state space.

Coming all the way from Feyman's diagrams <cit.>,
graphical languages are commonly used for representing quantum
processes. Whether directly based on quantum
circuits <cit.>
or stemming from categorical analysis such as the
ZX-calculus <cit.>, these formal languages are
still tied to the quantum coprocessor model in the sense that the only
monoidal structure that can be applied to quantum information is the
(multiplicative) Kronecker product. The only possible branching is
based on the (probabilistic) measurement.

Quantum Control Flow.

However, one peculiar feature of quantum computation is
non-causal execution paths.  Indeed, the Janus-faced quantum
computational paradigm features two seemingly distinct notions of
control structure. On the one hand, a quantum program follows
classical control: it is hosted on the conventional computer
governing the coprocessor, and can therefore only enjoy loops, tests
and other regular causally ordered sequences of operations. On the
other hand, the lab bench turns out to be more flexible than the rigid
coprocessor model, permitting more elaborate purely quantum
computational constructs than what quantum circuits or ZX-calculus

The archetypal example of a quantum computational behavior
hardly attainable within quantum circuits or ZX-calculus is the
Quantum Switch. Consider two quantum bits $x$ and $y$ and two
unitary operations $U$ and $V$ acting on $y$. The problem consists in
generating the operation that performs $UV$ on $y$ if $x$ is in state
$\ket{0}$ and $VU$ if it is in state $\ket{1}$. As $x$ can be in
superposition, in general the operation is then sending
$(\alpha\ket{0}+\beta\ket{1})\tensor\ket{y}$ to
$\alpha\ket{0}\tensor(UV\ket{y}) + \beta\ket{1}\tensor(VU\ket{y}).$ It
is a purely quantum test: not only can we have values in
superpositions (here, $x$) but also execution orders. This is
in sharp contrast with what happens within the standard quantum
coprocessor model.

Computational models supporting superpositions of execution orders
have been studied in the literature. One strand of research consists
in proposing a suitable extension of quantum
circuits <cit.>.
These approaches typically aim at discussing the notion of quantum
channel from a quantum information theoretical standpoint.

Limitation of Current Approaches and Objective of the Paper.

Although there is a finer and finer understanding of superposition of
causal orders in the literature, none of the existing PROPs can
support both the quantum switch on complex data built from tensors
and coproducts. We claim in this paper that the same intuition
underlying probabilistic branching can be followed for quantum (and
more general) branching. In the conventional case, $1\oplus 1$ is a
regular boolean: either “left” (standing e.g. for True) or
“right” (standing e.g. for False).

In quantum computation, the sum-type $1\oplus 1$ can however be
understood as a sum of vector spaces, giving an alternative
interpretation to $1\oplus 1$: it can be regarded as the type of a
quantum bit, superposition of True and False. One should note that
this appealing standpoint should be taken cautiously: (Pure) quantum
information imposes strong constraints on the structure of the data
in superposition: orthogonality and unit-norm have to be
preserved <cit.>.

The Quantum Switch can then be naturally understood in this
framework. Consider for instance
Figure <ref>, read from left to right: as
input, a pair of an element of type $A$ and a quantum bit. Based on
the value of the qubit (True or False), the wire $A$ goes in the upper
or the lower sheet, and is fed with $U$ then $V$ or $V$ then $U$. Then
everything is merged back together.

Quantum Switch with Worlds

In this paper, we introduce a new graphical language for quantum
computation, based on compact category with biproduct
<cit.>. This language allows us to express any
process with both pairing and a general notion of algebraic branching,
encompassing deterministic, non-deterministic, probabilistic and
quantum branching. We develop a denotational semantic and an
equational theory, and prove its soundness and completeness with
respect to the semantics. As a case-study, we show how the Quantum
Switch can naturally be encoded in the language.

In the paper, the missing proofs can be found in Appendix.

§ THE MANY-WORLDS CALCULUS

Our calculus is parameterized by a commutative semiring
$(R,+,0,\times,1)$. It can be instantiated by the complex numbers
$(\mathbb{C},+,0,\times,1)$ to represent pure quantum computations,
the non-negative real numbers $(\mathbb{R}_{\geq 0},+,0,\times,1)$ for
probabilistic computations, or the booleans
for non-deterministic computations.

While the goal is to define a graphical language in which each wire
can be enabled or disabled depending on the world in which the
computation takes place, we first define the category  of
diagrams without any "worlds", and will then add the world

§.§ A First Graphical Language

We define our graphical language within the paradigm of colored PROP
<cit.>, meaning that a diagram will be composed
of nodes, or generators, linked to each other through
colored wires, that are allowed to cross each
others. Additionally, we assume that our colored PROP is compact
closed and auto-dual, i.e. we allow to bend wires to obtain a
Cup (0.8$\tikzfig{lang/cup}$) or a Cap (0.8$\tikzfig{lang/cap}$).

The generators of our language are described in <ref>
and are respectively the Identity, the Swap, the Cup, the Cap, the
Plus, the Tensor, the Unit, the $n$-ary Contraction for $n \geq 0$,
and the Scalar for $s$ ranging over the commutative semiring $R$.
Mirrored versions of those generators are defined as syntactic sugar
through the compact closure, as shown for the mirrored Plus on the
right-hand-side of <ref>.  Diagrams are read
top-to-bottom: the top-most wires are the input wires and the
bottom-most wires are the output wires. The labels $A$, $B$,
etc correspond to the colors of our PROP. There is one color for every
type generated by the syntax[We do not assume any
associativity, distributivity, etc, although types that are
semantically equivalent will be made isomorphic through the equational
theory.] $A, B ::= \one \mid A\oplus B \mid A\otimes B$. As such, the
Unit starts a wire of type $\one$, the Plus combines two wires of type
$A$ and $B$ into a wire of type $A \oplus B$, and similarly the Tensor
combines  two wires of type $A$ and $B$ into a wire of type $A \tensor

In a colored PROP, an object $\goth{A}$ is simply a list of colors
$\goth{A} = A_1 \square \dots \square A_n$ (or $\goth{A} =
\emptysquare$ for the empty collection). For example, the Plus is a
morphism from $A \square B$ to $A \oplus B$. The choice of the
notation $\square$ for wires in parallel is uncommon, we use it to put
an emphasis on the fact that contrary to languages like the
ZX-calculus, wires that are in parallel are not necessarily "in tensor
with one another". In fact, $A \boxempty B$ can be understood
semantically as "either $A \otimes B$ or $A \oplus B$".

\( D_2\circ D_1:=\tikzfig{compo-seq}\qquad D_1\boxempty D_2:=\tikzfig{compo-ten} \)

Diagrams are obtained from generators by composing them in parallel
(written $\boxempty$), or sequentially (written $\circ$). Sequential
composition requires the type (and number) of wires to match.  We
write  for the category of diagrams we defined as such.

$\tikzfig{lang/id} \quad \tikzfig{lang/swap} \quad
\tikzfig{lang/cup} \quad \tikzfig{lang/cap} \quad
\tikzfig{lang/plus} \quad \tikzfig{lang/tensor-PN} \quad
\tikzfig{lang/unit} \quad \tikzfig{lang/contraction} \quad
\tikzfig{lang/scal} \hfill \tikzfig{lang/plus-inv} :=
\tikzfig{lang/plus-cc} $
Generators of our First Graphical Language ($n \geq 0$, $s \in R$)

We consider $R = \mathbb{C}$. While  lacks the worlds
labeling[We call worlds labeling the attribution of
multiple worlds to a single wire, as defined in
<ref>.], we can already illustrate our language by
encoding some basic quantum primitive in it and show how they operate.
In <ref> we show the encoding of a quantum bit $\alpha
\ket{0} + \beta \ket{1}$ and the Hadamard unitary. In particular, the
Plus allows to "build" a new quantum bit from two scalars in parallel
or to "open" a quantum bit to recover its corresponding scalars, the
left branch corresponding to $\ket{0}$ and the right branch to
$\ket{1}$. The meaning of the Contraction is better seen when applying
Hadamard to a quantum bit as we show in <ref>: it allows
us to duplicate and sum scalars. The rewriting sequence of
<ref> is made using the equational theory defined in
<ref>, however to correctly define our equational
theory, the worlds labeling are required. So while this specific
worlds-free rewriting sequence is sound, many other similar
worlds-free rewriting sequences are unsound.

$\alpha \ket{0} + \beta \ket{1} \rightsquigarrow \tikzfig{oneqbit} \in \catD(\emptysquare,\one \oplus \one) \hfill	\begin{pmatrix}
\frac{1}{\sqrt{2}} & \frac{1}{\sqrt{2}} \\\frac{1}{\sqrt{2}} & \frac{-1}{\sqrt{2}}
\end{pmatrix}\rightsquigarrow \tikzfig{hadamard} \in \catD(\one \oplus \one,\one \oplus \one)$
A Quantum Bit and the Hadamard Unitary

$\tikzfig{h-oneqbit-1} \hfill \to \hfill  \tikzfig{h-oneqbit-2}\hfill\to\hfill\tikzfig{h-oneqbit-3}\hfill\to\hfill\tikzfig{h-oneqbit-4}$
Applying the Hadamard unitary to a quantum bit

Instead of having the Cup and the Cap as generators and
defining the mirrored version of each generator through them,
one could proceed the other way around by defining the Cap as
follows, and the Cup in a mirrored
$\phantom{.}\hfill\tikzfig{cupcap/tens} \hfill
\tikzfig{cupcap/plus} \hfill
\tikzfig{cupcap/unit}\hfill\phantom{.}$

§.§ Adding Worlds Labeling

We now label wires of our diagram with sets of worlds $w \subseteq W$ from
a given a world set $W$. For each world $a \in W$, wires labeled
by a set containing $a$ are said to be "enabled in $a$", and the
others are said "disabled in $a$". This allows us to correlate the
enabling of wires. Before making this formal, we illustrate this notion through the following example:

\[\tikzfig{explicit/CNOT-Worlds}\]

Controlled Not with world set $\{ \aa,\bb, \cc,  \sstar\}$

\[ \tikzfig{explicit/agnostic_mono} \]

Example of Worlds-Agnostic Parallel Composition

The "controlled not" on quantum bits
can be represented by the <ref>. The figure is split into
two parts: the control part on the left-hand-side, and the
computational part on the right-hand-side.  The idea is that the
control part, that uses $\oplus$, will behave as an
if-then-else and will bind the world $\aa$ to the case where
the control quantum bit is $\ket{0}$, and the worlds $\bb$ and $\cc$
to the case where the control quantum bit is  $\ket{1}$. Lastly, the
world $\sstar$ appears nowhere in the labels, and corresponds to "we
do not evaluate this circuit at all"[While not strictly
necessary, it is often practical to have a world absent from every
wire.]. On the computational part, we apply the identity within the
world $\aa$, we negate $\ket{0}$ into $\ket{1}$ within $\bb$, and we
negate $\ket{1}$ into $\ket{0}$ within $\cc$. The domain and codomain
of this labelled diagram are $\left(\one \oplus \one: \{ \aa,\bb,
\cc\}\right) \boxempty \left(\one \oplus \one : \{ \aa,\bb,
\cc\}\right)$, which we will write $(\goth{A},\ell_{\goth{A}})$ with
an object $\goth{A} = (\one \oplus \one)\boxempty (\one \oplus \one)$
and a labeling function $\ell_{\goth{A}} : 1 \mapsto  \{ \aa,\bb,
\cc\} \quad 2 \mapsto \{ \aa,\bb, \cc\} $. Similarly, the "controlled
not" above can be seen as a diagram $\D_{\text{CNOT}}$ of
$\catD(\goth{A},\goth{A})$ together with a labeling function
$\ell_{\goth{A}}$ which labels every wire with a set of worlds.

We now give a formal definition of the concept of worlds:
Given a set of worlds $W$, we define the auto-dual compact closed
colored PROP $(\ManyWorlds{W},\boxempty,\emptysquare)$ of
many-worlds calculus over $W$ as follows:

Its colors are the pairs $(A:w)$ of colors $A$ of $\catD$ and
subsets $w \subseteq W$. We write $(\goth{A},\ell_{\goth{A}})$ for
the objects, where $\goth{A}$ is an object of $\catD$ and $\ell_A$
is a labeling function from the colors composing $\goth{A}$ to the subsets
of $W$.

Its morphisms $f \in \ManyWorlds{W}((\goth{A},\ell_{\goth{A}}),(\goth{B},\ell_{\goth{B}}))$ are pairs $(\D_f,\ell_f)$ of a morphism
$\D_f \in \catD(\goth{A},{\goth{B}})$ and a labeling function
$\ell_f$ from the wires of $\D_f$ to the subsets of $W$, satisfying
the following constraints: The label on an input or output wire of color $(A:w)$ must be equal
to $w$, and

\[\tikzfig{explicit/id-w}\hfill\tikzfig{explicit/swap-w}\hfill\tikzfig{explicit/cap-w}\hfill\tikzfig{explicit/cup-w}
\hfill\tikzfig{explicit/scalar-w} \hfill
\tikzfig{explicit/plus-w}\hfill\tikzfig{explicit/tensor-PN-w}\hfill\tikzfig{explicit/tensor-unit-w}\hfill\forall
n \geq 0,~\tikzfig{explicit/contraction-w}\]
where $\sqcup$
denotes disjoint set-theoretic union. The constraints for the
mirrored versions are similar.  The sequential composition $\circ$
and the parallel composition $\boxempty$ preserve the labels.

We write $[f]_W : \goth{A} \to \goth{B}$ for $f = (\D_f,\ell_f) \in
\ManyWorlds{W}((\goth{A},\ell_{\goth{A}}),(\goth{B},\ell_{\goth{B}}))$
where $\ell_{\goth{A}}$ and $\ell_{\goth{B}}$ can be deduced from
$\ell_f$ using the first restriction on labels. We note that when
considering $[f]_W : \goth{A} \to \goth{B}$ and $[g]_W : \goth{B} \to
\goth{C}$, there is no reason for the labels deduced on $\goth{B}$ to
be the same, so there is no reason for $[g \circ f]_W : \goth{A} \to
\goth{B}$ to be defined. More generally, when building two diagrams
$[f]_W$ and $[g]_V$, there is no reason for the worlds sets $W$ and
$V$ to be equal, but we might still want to be able to compose them
with one another.

We define the auto-dual compact closed colored PROP
$(\ManyWorlds{\forall},\boxempty,\emptysquare)$ of many-worlds
calculus as follows:

Its colors and objects are the same as the ones of .

Its morphisms from $\goth{A}$ to $\goth{B}$ are simply
morphisms $[f]_W : \goth{A} \to \goth{B}$ of $\ManyWorlds{W}$ for
any finite set $W$. See <ref> for the canonical
labelings on the identity, swap, cup and cap. Morphisms are
considered up to renaming of the worlds[More precisely,
$[(\D_f,\ell_f)]_W = [(\D_f,\sigma\circ \ell_f)]_{\sigma(W)}$ for
any $\sigma$ describing a bijection between $W$ and $\sigma(W)$.
Without this quotient, we would not have $[f]_W \circ \id =

The parallel composition is given by $[f]_W \boxempty [g]_V
:= [f^{- \times V} \boxempty g^{W \times -}]_{W \times V}$ where
$f^{\sigma(-)}$ has the same diagram $\D_f$ and has for labels
$\ell_{f^{\sigma(-)}}(x) = \sigma(\ell_f(x))$.

The sequential composition is given by $[g]_V \boxempty
[g]_W := [g^{(W \times -)\cap Z} \circ f^{(- \times V) \cap
Z}]_{Z}$ where $Z \subseteq W \times V$ is the greatest subset
such that this composition is well-defined, as explained in the
following example.

\[ \begin{matrix}
W=\{a,\star\} && W = \{a,b,c,\star\} && W = \{a,\star\} && W = \{a,\star\} \\
\tikzfig{explicit/id-can} && \tikzfig{explicit/swap-can} && \tikzfig{explicit/cup-can} &&  \tikzfig{explicit/cap-can}
\end{matrix} \]

Canonical labelings in $\ManyWorlds{\forall}$

[World Agnostic Composition]
In <ref>, we show the result of the parallel
composition of $[f]_{\{a,\star\}} = \id_{A:\{a\}}$ and
$[g]_{\{b,\star\}} = \id_{A:\{b\}}$ for a color $A$: the world
$(\aa,\sstar)$ corresponds to the left wire being enabled and the
right one disabled, the world $(\sstar,\bb)$ is the opposite, the
world $(\aa,\bb)$ is both enabled and the world $(\sstar,\sstar)$
is both disabled.

Then, in <ref> we continue by composing
with the Cup over $A:\{\cc,\sstar\}$. To compute the composition,
we proceed in two steps: first we handle the situation as if it
were a parallel composition, leading to a diagram labelled over $W
\times V \times U$, but with multiple contradictory labels on the
wire. Then, we eliminate as few worlds as possible to make those
labels compatible:

* We eliminate $(\aa,\bb,\sstar)$ and
$(\aa,\sstar,\sstar)$ which are on the left label but
not on the bottom one.

* We eliminate $(\sstar,\bb,\cc)$ and
$(\sstar,\sstar,\cc)$ which are on the bottom label but
not on the left one.

* Looking at the right, we eliminate
$(\sstar,\bb,\sstar)$ and $(\aa,\sstar,\cc)$ for similar

We eliminated six worlds, with the only remaining ones being $Z =
\{(\aa,\bb,\cc),(\sstar,\sstar,\sstar)\}$.

\[\tikzfig{explicit/agnostic_compo}  \]

Example of Worlds-Agnostic Sequential Composition

§.§ Representing Computation

As a motivational example, we may see how this Many-Worlds diagrams can
be used to represent computations expressed in a language that
explicitely uses the two compositions $\otimes$ (through pairs), and
$\oplus$ (through pattern-matching), as in <cit.>.
As this translation is not the focus of this paper, a lot of
technicalities are glossed over.

§.§.§ Syntax of the Language

The language we present here is adapted from
<cit.>, where we consider the values and types
together with the enrichment that is linear combinations of terms, but
without abstraction nor recursion. The syntax that is used in the
language is given as follows, with scalars $s$ ranging over the
commutative semiring $R$:

(Base types)       A, B          ::=         A ⊕B
A ⊗B

(Isos, first-order)     α         ::=      AB

(Values)     v          ::=       ⟨⟩x v v
v_1…, v_n

(Expressions)     e          ::=       v x_1…,x_nω x_1…,x_ne e+e

(Isos)              ::=

(Terms)     t          ::=       ⟨⟩x t t

αt t + t  t x_1…,x_ntt

The language in particular features branching through the $\inl{}$ and
$\inr{}$ constructors, linear combinations of terms and expressions,
and crucially isos that have type $A\leftrightarrow B$: they
turn a term of type $A$ to a term of type $B$ using pattern-matching.
The language comes with a predicate (not presented here, although it
could be given a diagrammatic meaning), used in the typing rule of
isos, to ensure exhaustivity and the non-overlapping character of the
left-hand and right expressions of the clauses, allowing in particular
to define unitaries (in the complex setting). Constraints on the
linear combinations may also be used to enforce probabilistic
constraints (i.e. that states are normalised in the quantum setting).

There are two different typing judgements, one for terms $\vdash_t$
and one specific for isos $\vdash_\omega$.

For a more complete description of the language, we refer the reader
to <cit.>. In the following, we will use the
shorthands $\textbf{ff}:=\inl\langle\rangle$ and

In the case where $R = \mathbb C$, one can encode the Hadamard gate by:

\scalebox{0.8}{$\left\{
\begin{array}{l@{~}c@{~}l}
\textbf{ff} & {\isot} & \frac{1}{\sqrt{2}} (\textbf{ff} + \textbf{tt}) \\
\textbf{tt} & {\isot} & \frac{1}{\sqrt{2}} (\textbf{ff} - \textbf{tt}) \\
\end{array}
\right\} : \one\oplus\one\isot\one\oplus\one$}

In this setting, any quantum circuit can be encoded by an iso.

§.§.§ Encoding into the Many-Worlds

One can encode any term of the language into a Many-Worlds diagram.
Given some typing derivation $\xi$ of a term $x_1 : A_1, \dots, x_n :
A_n \vdash_t t : B$ we write  $\pparenthesis{\xi}$ for the function
that maps $\xi$ to a diagram in the Many-Worlds Calculus with $n$ input
wires of type $A_1, \dots, A_n$ and one output wire of type $B$.
For the typing derivation $\xi$ of an iso $\vdash_\omega \omega :A
\isot B$, $\pparenthesis{\xi}$ gives a diagram with one input wire of
type $A$ and one output wire of type $B$.

$\pparenthesis{\xi}$ is defined inductively over $\vdash_t$ and
$\vdash_\omega$ as shown in appendix in <ref>,
where  the worlds sets are handled by compositions of world-agnostic

Representation of the term $t$ from

We can represent the term

\[t:=\left\{\begin{matrix} \langle\textbf{tt},x\rangle
\leftrightarrow & \textbf{let }y=\textup{H }x\textup{ in
}\frac{1}{\sqrt{2}}\langle \textbf{tt},y\rangle +
\frac{1}{\sqrt{2}} \langle \textbf{ff},y\rangle\\
\langle\textbf{ff},x\rangle \leftrightarrow & \textbf{let
}y=\textup{Id }x\textup{ in }\frac{1}{\sqrt{2}}\langle
\textbf{tt},y\rangle - \frac{1}{\sqrt{2}} \langle
\textbf{ff},y\rangle \end{matrix}
\right\}:(\one\oplus\one)^{\otimes2}\leftrightarrow(\one\oplus\one)^{\otimes2}\]

(already given in <cit.>) as shown in
<ref>. The yellow box stands for the Hadamard gate
(which one can build following <ref>).
Each line of this isomorphism corresponds to a column of the
figure. Each columns start by matching the input as
$\langle\textbf{tt},x\rangle$ or $ \langle\textbf{ff},x\rangle$,
then compute $y$ from $x$, and then build the output by following
the syntax. The world set is computed by composing each of the
blocks of this term, using the world-agnostic composition of
$\ManyWorlds{\forall}$. It can be seen as a subset of $\{a,b\}
\times \{c,c'\} \times \{0,1,2,3\}$ where $\{a,b\}$ corresponds to
being on the first or second line of the matching, $\{c,c'\}$
being on the left or right of the sum, and $\{0,1,2,3\}$ being the
world set of the Hadamard gate.

§.§ Comparison with Other Graphical Languages

The distinctive feature of the Many-Worlds Calculus is that it
graphically puts the tensor and the biproduct on an equal footing. By
comparison, other graphical language for quantum computing are
inherently centered around either one of them. The
ZX-calculus <cit.> and cousin languages ZW- and
ZH-Calculi <cit.>, as
well as Duncan's Tensor-Sum Logic <cit.>, use
the tensor product as the default monoid, while more recent language
– particularly for linear optics <cit.> – use the
biproduct. We have a closer look at each of them in the following, and
show how – at least part of – each language can be encoded naturally
in the Many-Worlds Calculus. Most of them comes equipped with an
equational theory. By completeness of our language
(<ref> to follow), all the equations expressible in
the fragments we consider can be derived in our framework.

§.§.§ ZX-Calculus

\(\tikzfig{trad-zx-gn} \)

The first difference is the restrictions of the ZX-calculus to
computations between qubits, in other words linear map from
$\mathbb{C}^{2^n} \mapsto \mathbb{C}^{2^m}$, while our language can
encode any linear map from $\mathbb{C}^{n} \mapsto \mathbb{C}^{m}$.
The Tensor generator allowing the decomposition of $\mathbb{C}^{2^n}$ into
instances of $\mathbb{C}^2$ was already present in the scalable
extension of the ZX-calculus <cit.>, but the main
difference comes from the Plus (and the Contraction).

Additionally, every ZX-diagram can be encoded in our graphical
language. The identity, swap, cup and cap of the ZX calculus are
encoded by the similar generators over the type $\one \oplus \one$,
the Hadamard gate is encoded as in <ref>, and the green
spider is encoded as shown on the right. An encoding for the red
spider can then be deduced from those. Diagrams are composed them
together with the world-agnostic composition of

§.§.§ Tensor-Sum Logic

\(\tikzfig{box-contraction} \)

The core difference between their work and ours is the presence of the
contraction in our graphical language. They instead rely on an
enrichment of their category by a sum, which they represent
graphically with boxes. We show on the right how the morphism $f+g$
would be encoded in both their and our language. More generally, their
boxes correspond to uses of our contraction generator in a
"well-bracketed" way. Another point of difference is their approach
to quantum computation, as we do not assign the same semantics to
those superpositions of morphisms. In their approach, the superposition
is a classical construction and corresponds to the measurement and the
classical control flow, while in our approach the superposition is a
quantum construction and corresponds to the quantum control.

§.§ PBS-Calculus

The PBS Calculus <cit.> allow one to represent coherent quantum
control by the use of polarising beam splitters (pbs): whenever
a qbit enter a pbs node, depending on the polarity of the qbit it will
either go through or be reflected. By making implicit the target
system, controlled by the optical system represented by the diagram,
the PBS-Calculus allows one to encode the Quantum Switch (depicted on
the right). The pbs generator is related to the $\oplus$ of the

The first main difference with
our language is that, since the generators of the PBS-Calculus represent
physical components, any PBS-diagram is by construction physical, while our
approch is more atomic and decomposes physical components into abstract
smaller ones. The second main difference lie in the trace: while they can
allow a particle to
pass through a wire at most twice, in our system, each wire can be
used at most once: more formally, their trace is based on the
coproduct while ours is on the tensor product. If we are assured that
each wire can only be used once during the computation, any PBS-diagram can
be translated to the Many-Worlds calculus direclty, with the transformation on the right,
where we distinguish the control system (the part of the diagram
connected to [style=plus] ; s) from the target system
(connected to [style=contraction] ; s) which is
implicit in the PBS-Calculus.

§.§ LOv-Calculus

In the PBS-Calculus, the qubit in the control system (the one explicitely
represented) cannot be put in arbitrary superpositions of $\ket0$ and $\ket1$
during the computation. To allow this feature, we may add some linear
optical components to the language's generators, and end up with the
LOv-Calculus <cit.>. In this language, there is no trace and there is a
unique photon traveling the circuit, which relieves us of
the previous constraint. There is also no need for an implicit target system
anymore. All wires at the interface between the generators are of type
$\one\oplus\one$, and parallel wires have disjoint sets of worlds. Each
generator can then be interpreted as follows:

\[
\tikzfig{lov-to-mwc-wave-plates}
\qquad\qquad
\tikzfig{lov-to-mwc-BS}\]

\[\tikzfig{lov-to-mwc-PBS}
\qquad\qquad
\tikzfig{lov-to-mwc-phase-shifter}
\qquad\qquad
\tikzfig{lov-to-mwc-vacuum-source}\]

§.§ Path-Calculus

The Path-Calculus is another recent graphical language for linear optical
circuits <cit.>. Its generators correspond directly to a subset of the
Many-Worlds' with path-to-mwc-unary-contraction,
path-to-mwc-binary-contraction and path-to-mwc-scalar;
where each wire has type $\one$ and where parallel wires are on disjoint sets
of worlds. This language is then used as the core for a more
expressive language called QPath, which this time cannot be direclty encoded
in our language, except when restricting the set of generators (specifically
to $n=1$), in which case path-to-mwc-ket1.

§ SEMANTICS OF THE MANY-WORLDS CALCULUS

Our calculus represent linear operators between finite dimension
$R$-modules $(\FdM)$. In particular, in the case $R = \mathbb{C}$ we use linear
operator between finite dimension Hilbert spaces, which correspond to
pure quantum computations. More precisely, we will define two
semantics, a world-dependent semantics $\interp{-}_a$ for every world
$a \in W$, which will be a monoidal functor from $\ManyWorlds{W}$ to
$\FdM$, and a worlds-agnostic semantics $\interp{-}$ from
$\ManyWorlds{\forall}$ to $\FdM$.

We start by defining those semantics on the objects. For every object
$\goth{A}$ of , we define its enablings $\enab{\goth{A}}$
as "replacing any number of wire type by $\bullet$". For example
$\enab{(A \boxempty B)} = \{ \bullet \boxempty \bullet, \bullet
\boxempty B, A \boxempty \bullet, A \boxempty B \}$.
Then, for any object $(\goth{A},\ell_{\goth{A}})$ of $\ManyWorlds{W}$
and any world $a \in W$, we define
$\enabW{(\goth{A},\ell_{\goth{A}})}{a}$ to be the enabling of
$\goth{A}$ in which every $(A:w)$ with $a \in w$ is preserved and
every $(A:w)$ with $a \notin w$ is replaced by $\bullet$. For example
$\enabW{(A:\{a\}\boxempty B:\{b\})}{a} = A \boxempty \bullet$.
To each enabling
$\goth{E} \in \enab{\goth{A}}$ we associate a $R$-module
$\mathcal{M}_{\goth{E}}$ as follows:

\begin{matrix}
\mathcal{M}_{A_1\boxempty \dots \boxempty A_n} := \mathcal{M}_{A_1}
\otimes \dots \otimes \mathcal{M}_{A_n}
\mathcal{M}_{A \otimes B} := \mathcal{M}_{{A}} \otimes
\mathcal{M}_{{B}}
\\
\mathcal{M}_{\emptysquare} := R \qquad  \mathcal{M}_{\bullet}
:= R \qquad \mathcal{M}_{\one} := R
\mathcal{M}_{A \oplus B} := \mathcal{M}_{{A}} \oplus \mathcal{M}_{{B}}
\end{matrix}$

We can
then define the semantics $\interp{-}_a:\ManyWorlds{W} \to
\FdM$
and $\interp{-}:\ManyWorlds{\forall} \to \FdM$ on objects
as $\interp{(\goth{A},\ell_{\goth{A}}))}_a :=
\mathcal{M}_{\enabW{(\goth{A},\ell_{\goth{A}})}{a}}$
$\interp{\goth{A}} := \bigoplus_{\goth{E} \in \enab{\goth{A}}} \mathcal{M}_{\goth{E}}$.

Then, for the morphisms, we proceed by compositionality for
$\interp{-}_a$, meaning that we define $\interp{-}_a$ on every
generator and compute the semantics of a diagram by decomposing it
$\interp{g \circ f}_a := \interp{g}_a \circ \interp{f}_a$ and
$\interp{f \boxempty g}_a := \interp{f}_a \otimes \interp{g}_a$.

The semantics of most generators is given in <ref>; for
all the other generators their semantics is simply the identity operator.

\[\interp{\tikzfig{explicit/swap-w}}_a = \left\lbrace
\begin{matrix}
\textup{Id}
&\in \FdM(\mathcal{M}_{A},\mathcal{M}_A) &\text{ if }a \in w\backslash v \\
\textup{Id}
&\in \FdM(\mathcal{M}_{B},\mathcal{M}_B) &\text{ if }a \in v\backslash w \\
h \otimes h' \mapsto h' \otimes h &\in \FdM(\mathcal{M}_{A\otimes B},\mathcal{M}_{B \otimes A}) &\text{ if }a \in w \cap v\\

(1) &\in \FdM(R,R) &\text{ otherwise} \end{matrix}\right. \]

\[\interp{\tikzfig{explicit/cup-w}}_a = \left\lbrace
\begin{matrix}
h \otimes h' \mapsto \langle h | h' \rangle &\in \FdM(\mathcal{M}_{A\otimes A},R) &\text{ if }a \in w \\
(1) &\in \FdM(R,R) &\text{ otherwise} \end{matrix}\right. \qquad  \text{\parbox{3cm}{and the (conjugate) transposed operator for the Cap}}
\]

\[\interp{\tikzfig{explicit/scalar-w}}_a = \left\lbrace
\begin{matrix}
s \cdot \textup{Id} &\in \FdM(\mathcal{M}_A,\mathcal{M}_A) &\text{ if }a \in w \\
(1) &\in \FdM(R,R) &\text{ otherwise} \end{matrix}\right.
\]

\[\interp{\tikzfig{explicit/plus-w}}_a = \left\lbrace
\begin{matrix} \begin{pmatrix}
\textup{Id} \\ 0
\end{pmatrix} &\in \FdM(\mathcal{M}_{A\oplus B},\mathcal{M}_A) &\text{ if }a \in w \\
\begin{pmatrix}
0 \\ \textup{Id}
\end{pmatrix} &\in \FdM(\mathcal{M}_{A\oplus B},\mathcal{M}_B) &\text{ if }a \in v\\ (1) &\in \FdM(R,R) &\text{ otherwise} \end{matrix}\right. \qquad \qquad \text{\parbox{3cm}{and the (conjugate) transposed operator for the mirrored Plus}} \]

Semantics of the Generators of $\ManyWorlds{W}$ in a World $a \in W$.

The worlds-agnostic semantics is defined from the world-dependent
semantics, as follows. Consider
$f\in\ManyWorlds{W}((\goth{A},\ell_{\goth{A}}),(\goth{B},\ell_{\goth{B}}))$, we define:

\[\interp{[f]_W} := \left\{
\text{\raisebox{15pt}{$\sum_{\substack{a\in W
\\\enabW{(\goth{A},\ell_{\goth{A}})}{a}=\goth{A}'\\\enabW{(\goth{B},\ell_{\goth{B}})}{a}=\goth{B}'}}
\interp{f}_a$}}\right\}_{\goth{A}' \in \goth{A}^{\bullet},
\goth{B}' \in \goth{B}^{\bullet}} \in \FdM\left(
\bigoplus_{\goth{E} \in \goth{A}^{\bullet}}
\mathcal{M}_{\goth{E}}, \bigoplus_{\goth{F} \in
\goth{B}^{\bullet}} \mathcal{M}_{\goth{B}} \right)\]

For example, the worlds-agnostic semantics of the Tensor and the Plus
are simply the collection of all their world-dependent semantics
assembled into a single linear operator:

\[\interp{\tikzfig{explicit/tensor-PN}} = \semTensorBin\]

\[\interp{\tikzfig{explicit/plus}} = \semPlus\]

The worlds-agnostic semantics is functorial (see <ref>) and is universal in the following sense:

For every $\goth{A},\goth{B}$ objects of $\catD$, and for every
linear operator in $U \in \FdM\left( \bigoplus_{\goth{E} \in
\goth{A}^{\bullet}} \mathcal{H}_{\goth{E}}, \bigoplus_{\goth{F}
\in \goth{B}^{\bullet}} \mathcal{H}_{\goth{B}} \right)$, there
exists a world set $W$, and a morphism $[f]_W \in
\ManyWorlds{\forall}(\goth{A},\goth{B})$ such that $\interp{f} =

§ THE EQUATIONAL THEORY

Equations with a Fixed World Set $W$

\[\tikzfig{eq/all_complete}\]

Equations with Side-Effects on World Sets

Similarly to how our semantics is defined in two steps, the
equational theory is also defined in two steps:

* A set of equations within $\ManyWorlds{W}$ for a fixed set of worlds
$W$, which will not be complete, but will be sound for $\interp{-}_a$
for every $a \in W$, hence sound for $\interp{-}$ too. We write
$\equiv_{W}$ for the induced congruence[In other words the
smallest equivalence relation satisfying those equations and such that
$f \equiv_{W} f' \implies \forall g,h,l, g \circ (f \boxempty h) \circ
k \equiv_{W} g \circ (f' \boxempty h) \circ k$.] over
$\ManyWorlds{W}$. We list those equations in <ref>.
Quite notably, the last two rows describe the fact that the
contraction is a natural transformation.

* Five additional equations with side effects on the set of worlds,
which will be sound and complete for $\interp{-}$ (but not for $\interp{-}_a$). We write $\equiv$ for the induced equivalence
relation, which is a congruence over $\ManyWorlds{\forall}$.  We have: One equation
that allows us to rename the worlds: for every morphism
$(\D_f,\ell_f)$ of $\ManyWorlds{W}$, and for every bijection $i : W
\to V$, we have $[(\D_f,\ell_f)]_W \equiv [(\D_f,i \circ \ell_f)]_V$; Two equations allowing the annihilation (or creation,
when looking at them from right to left) of worlds due to coproducts
or scalars (first row of <ref>); Two equations
allowing the splitting (or merging, when looking at them from right to
left) of worlds due to coproducts or scalars (second row of

For $f$ a morphism of $\ManyWorlds{W}$ and $g$ of
$\ManyWorlds{V}$, whenever $[f]_W \equiv [g]_V$ we have $\interp{[f]_W} =
\interp{[g]_V}$. Additionally if $W = V$, whenever $f \equiv_W g$
we have $\forall a \in W, \interp{f}_a = \interp{g}_a$.

For every $[f]_W : {\goth{A}} \to {\goth{B}}$ and $[g]_V: {\goth{A}} \to {\goth{B}}$, $\interp{[f]_W} = \interp{[g]_V} \iff [f]_W \equiv

The soundness can be proved by a case-by-case analysis on every
equation.  The completeness theorem follows from the existence of a
normal form for $\equiv$, as described in <ref>.

[The Quantum Switch] Similarly to the
"controlled not" in <ref> where a quantum bit controls
whether the identity or a negation is applied to some data, one
can consider the case where a quantum bit control whether $U \circ
V$ or $V \circ U$ is applied to some data (for $U$ and $V$ two
quantum operators on the same type $A$).

On the rightmost part of <ref>, one can see the
most direct implementation of this "quantum switch", but this
implementation uses two copies of $U$ and $V$. On the leftmost
part of the figure, we show another implementation which only
relies on one copy of each, and both can be rewritten into another
using the equational theory.

In those diagrams, the world set is $W = w \sqcup v$ and we rely on
violet, blue and red wires to indicate respectively worlds labels
$w \sqcup v$, $w$ and $v$. Each figure has a control side which
operates on a quantum bit (type $\one \oplus \one$) and binds the
world $w$ to $\ket{0}$ and the world $v$ to $\ket{1}$; and a
computational side which operates on some data of an arbitrary type
$A$, on which could be applied $U$ and/or $V$ which stand for two
morphisms of $\ManyWorlds{W}(A:W,A:W)$.

The first rewriting step relies on the two lemmas on the right,
both of which being deducible from the equational theory (see
<ref>). The second rewriting step is simply using
the properties of a compact closed category.
$\tikzfig{explicit/lemma_nat} \qquad \tikzfig{explicit/lemma_c}$

\[ \tikzfig{explicit/legend}   \hfill \tikzfig{explicit/QS_unique}
\quad \equiv_{w \sqcup v} \quad \tikzfig{explicit/QS_middle} \quad
\equiv_{w \sqcup v} \quad \tikzfig{explicit/QS_double} \]

Rewriting the Quantum Switch

§ CONCLUSION

We introduced a new sound and complete graphical language based on
compact categories with biproducts, along with an equational theory and
a worlds system, helping us to build a denotational semantics of our

This language is a first step towards the unification of languages
based on the tensor $\otimes$ and those based on the biproduct
$\oplus$. This allows us to reason about both systems in parallel, and
superposition of executions, as shown by the encoding of the Quantum
Switch in Example <ref> and the translation from term of the
language in <ref>.

Following this translation, a natural development of the Many-Worlds
calculus consists in accomodating function abstraction and recursion
in the language. The question of a complete equational theory for the
language on mixed states (e.g. via the discard construction
<cit.>) is also left open.

Finally, while our language allows for quantum control, it does not
entirely capture another language that aims at formalizing quantum
control, namely the PBS-Calculus <cit.>. How and in which context
could we capture the PBS-Calculus is left for future work.

§ ENCODING THE ISO LANGUAGE WITH THE MANY WORLDS CALCULUS

For simplicity of the encoding, we give the following syntactic

$\pparenthesis{\infer{x : A \vdash x : A}{}} = \tikzfig{terms/var}$$\pparenthesis{\infer{\vdash \langle\rangle : \one}{}} = \tikzfig{terms/unit}$$\pparenthesis{\infer{\Delta\vdash\alpha t : A}{\infer{\Delta\vdash t : A}{\xi}}} = \tikzfig{terms/alpha}$

$\pparenthesis{\infer{\Delta \vdash \inr{t} : A\oplus
B}{\infer{\Delta \vdash t : A}{\xi}}} = \tikzfig{terms/right}$$\pparenthesis{\infer{\Delta \vdash \inl{t} : A\oplus
B}{\infer{\Delta \vdash t : B}{\xi}}} = \tikzfig{terms/left}$

$\pparenthesis{\infer{\Delta_1,\Delta_2\vdash \pv{t_1}{t_2} : A\otimes
B}{\infer{\Delta_1\vdash t_1 : A}{\xi_1}\qquad \infer{\Delta_2\vdash
t_2 : B}{\xi_2}}}= \tikzfig{terms/tensor}$$\pparenthesis{\infer{\Delta\vdash t_1 + t_2 : A}{\infer{\Delta\vdash t_1 :
A}{\xi_1} \qquad \infer{\Delta\vdash t_2 : A}{\xi_2}}} =
\tikzfig{terms/plus}$

$\pparenthesis{\infer{\Delta_1, \Delta_2 \vdash \letpv{x_1, \dots, x_n}{t_1}{t_2} :
B}{\infer{\Delta_1 \vdash t_1 : A_1\otimes \dots \otimes A_n}{\xi_1}
\qquad\infer{x_1 : A_1, \dots, x_n : A_n, \Delta_2 \vdash t_2 :
B}{\xi_2}}} = \tikzfig{terms/let}$

$\pparenthesis{\infer{\vdash \isobasique : A \isot B}{\infer{\Delta_i \vdash v_i :
A}{\xi_i}\qquad \infer{\Delta_i \vdash e_i : B}{\xi'_i}}} =
\tikzfig{terms/iso}$

$\pparenthesis{\infer{\Delta \vdash \omega~t : B}{\infer{\vdash_\omega \omega :
A\isot B}{\xi_1} \qquad \infer{\Delta\vdash t : A}{\xi_2}}} =
\tikzfig{terms/app}$
Inductive translation of terms from <ref> into Many-Worlds diagrams.

§ NORMAL FORM, UNIVERSALITY AND COMPLETENESS

We can prove that the equational theory of <ref> is
complete, by defining a normal form on the morphisms, and by showing
that all diagrams can be put in this normal form, which is unique.

§.§ Normal Form

It will be practical for our normal forms to define the following
syntactic sugar, which we call the unitor, its unit and
its generalized form:
\hfill\tikzfig{unit-of-unitor}\hfill
\tikzfig{generalized-unitor}~:=~\tikzfig{generalized-unitor-def} $

\(
\tikzfig{NF-matrix-block-gen} \)

We define the short-hand on the right with the assumption that the
world set is in bijection with the set of scalars, in other words
the scalars $\lambda_{i\,j}$ have for worlds label
$\{a_{i\,j}\}$. In
particular, all the input (resp. output) wires live in mutually
exclusive worlds.

A first important observation is that we can put this diagram in the following alternative form, with exactly the same assumption on the world set:

\begin{equation}
\label{eq:NF-matrix-block-alt}
\tikzfig{NF-matrix-block-gen-alt}
\end{equation}

We will freely change between the two representations in the following.
Another important observation is that any permutation
of wires (all mutually exclusive, and of type $\one$) can easily be
put in this form using the following equations:

$\phantom{.}\hfill\tikzfig{NF-dangling-scalar-0} \hfill
\tikzfig{NF-0-branch} \hfill
\tikzfig{1-1-contraction-is-identity}\hfill\phantom{.}$
The normal
form of a morphism $f:\goth{A}\to \goth{B}$ is defined as the form of
the diagram on the left of <ref>, where the
morphism $\iso_{\goth{A}}$ is defined inductively right.

$f = \tikzfig{NF-gen} \quad $ 0.9$\tikzfig{iso-all}$
The Normal Form (left) and the Inductive Definition of
$\iso_{\goth{A}}$ (right)

The output wires of $\iso_{\goth{A}}$ for any $\goth{A}$ live in
mutually exclusive worlds, but once again, we don't overload the
diagrams with unitors or world names encoding this information,
although it will be used in the following.

Notice that graphically, there is no difference between $\goth
A\boxempty(\goth B\boxempty \goth C)$ and $(\goth A\boxempty \goth
B)\boxempty \goth C$, in other words $\boxempty$ is strictly
associative. However $\iso_{\goth A\boxempty(\goth B\boxempty \goth
C)}$ and $\iso_{(\goth A\boxempty \goth B)\boxempty \goth C}$ are
different, but they are equivalent up to a rearranging of the output

There exists a wire permutation $\sigma$ such that $\iso_{\goth
A\boxempty(\goth B\boxempty \goth C)}=\sigma\circ\iso_{(\goth
A\boxempty \goth B)\boxempty\goth C}$.

notice that in both $\iso_{\goth A\boxempty(\goth B\boxempty \goth
C)}$ and $\iso_{(\goth A\boxempty \goth B)\boxempty \goth C}$, we can
use the bialgebra between contractions and unitors, followed by their
respective fusions in the following way:

\begin{align*}
\tikzfigc{00}
\equiv\tikzfigc{01}
\equiv\tikzfigc{02}\\
\equiv\tikzfigc{03}
\end{align*}

and similarly for $\iso_{(A\boxempty B)\boxempty C}$. It then suffices
to check which contractions the bottom unitors are linked to. Naming
the $i$th contraction existing $\iso_A$ as $a_i$, and similarly for
$B$ and $C$, we can see that for each triple $(a_i,b_j,c_k)$ there is
exactly one unitor connected to precisely contraction $a_i$, $b_j$ and
$c_k$, in both diagrams. The same is true for every pair $(a_i,b_j)$,
$(a_i,c_k)$ and $(b_j,c_k)$, as well as for every 1-tuple $(a_i,)$,
$(b_j,)$ and $(c_k,)$. This shows that both diagrams are equal up to
rearranging of the outputs.

We hence have a choice to make here for canonicity, and choose
$\iso_{A_0\boxempty A_1\boxempty A_2\boxempty
...}:=\iso_{(...((A_0\boxempty A_1)\boxempty A_2)\boxempty ...)}$.

We define $\iso_{\goth{A}}^{\inv}$ inductively in the same way, but upside-down.

We note that $\iso_{\goth{A}}^{-1}\circ\iso_{\goth{A}}$ is the normal form of $\textup{id}_{\goth{A}}$.

Notice that iso-plus-wrt-boxempty and iso-tensor-wrt-boxempty.

We have the following identities:

\[\tikzfig{isos-are-isos}~~\text{ and }~~\tikzfig{ok-isos-are-actually-pseudo-isos} \text{where $\ell(s_i)\cap\ell(s_j)=\varnothing$ when $i\neq j$.}\]

We will use the following identities:

\begin{align}
\tikzfig{isos-are-isos-boxempty} \begin{array}{l} \text{ when all the wires }a_{ij}, b_i,c_j\\\text{ are mutually exclusive}\end{array}
\end{align}

We can show this result by induction on $n$ and $m$. Case $(0,m)$ is obvious. Case $(1,1)$ can be proven easily using world sets:

\begin{align*}
\tikzfigc{00}
\end{align*}

For any $n$ and $m$, we can then prove the case $(n+1,m)$ using the cases $(n,m)$ and $(1,m)$ (the case $(n,m+1)$ is completely symmetric):

\begin{align*}
\tikzfigc{00}
\end{align*}

In a similar way, it is possible to show the following three identities:

\[\tikzfig{isos-are-isos-boxempty-2}\begin{array}{l} \text{ when all the wires }a_{ij}, b_i,c_j\\\text{ are mutually exclusive}\end{array}\]

\[\tikzfig{isos-are-isos-tensor}\begin{array}{l} \text{ when all the wires }a_{ij}\\\text{ are mutually exclusive}\end{array}\]

\[\tikzfig{isos-are-isos-tensor-2}\begin{array}{l} \text{ when all the wires }a_{ij}\\\text{ are mutually exclusive}\end{array}\]

*  [$\iso_{\goth A}^{\inv}\circ\iso_{\goth A}$]: The result is obvious in cases $\one$ and $\emptysquare$. For $A\oplus B$:

\begin{align*}
\tikzfigc{00}
\equiv\tikzfigc{01}
\equiv\tikzfigc{02}
\equiv\tikzfigc{03}
\end{align*}

The proof is similar for $\otimes$ and $\boxempty$ using the previous identities.

*  [$\iso_{\goth A}\circ\iso_{\goth A}^{\inv}$]: The result is again obvious for $\one$ and $\emptysquare$. The general result is easy to prove by induction using the above identities.

§.§ Universality

Now that we have defined a normal form, we may use it to show the universality (<ref>) of the language, i.e. that for every $\goth{A},\goth{B}$ objects of $\catD$, and for every linear operator

\[ \Lambda = \begin{pmatrix}
\lambda_{1\,1} & \cdots & \lambda_{n\,1} & \lambda_{10} \\
\vdots & & \vdots & \vdots \\
\lambda_{1\,m} & \cdots & \lambda_{n\,m} & \lambda_{m0} \\
\lambda_{01} & \cdots & \lambda_{0n} & \lambda_{00}
\end{pmatrix} \in \FdM\left(
\bigoplus_{\goth{E} \in \goth{A}^{\bullet}}
\mathcal{H}_{\goth{E}},
\bigoplus_{\goth{F} \in \goth{B}^{\bullet}}
\mathcal{H}_{\goth{B}}
\right)\]

then the morphism $f = \iso_{\goth{B}}^{-1} \circ \lambda \circ \iso_{\goth{A}}$ defined in <ref> satisfies $\interp{f} = \Lambda$.

As stated in the corresponding section, there is one world for each
scalar in $\lambda$, so we write $a_{i\,j}$ the world associated to
$\lambda_{i\,j}$  and $W$ the set of all those worlds. Let us write
$\one^k_0=\bullet^{k}$ for $\bullet \boxempty \dots \boxempty \bullet$
(with $k$ elements) and $\one^k_i$ for $\bullet^k$ where the $i$-th
$\bullet$ has been replaced by $\one$ if $1\leq i\leq k$. Beware that
we consider $\one^k_0$ to be the last element of the canonical
basis. We then have $\interp{\lambda}_{a_{i\,j}} :\one^n_i \to
\one^m_j : x \mapsto \lambda_{i\,j}  \cdot x$ where:

\[\tikzfig{NF-matrix-block-gen}\]
Additionally, one can show by
induction that $\interp{\iso_{\goth{A}}}_{a_{i\,j}}$ is simply the
projection on the $i$-th element of the canonical basis, and
$\interp{\iso^{-1}_{\goth{B}}}_{a_{i\,j}}$ is simply the injection on
the $j$-th element of the canonical basis. Since we have $\interp{f}_a
= \interp{\iso^{-1}_{\goth{B}}}_{a} \circ \interp{\lambda}_a \circ
\interp{\iso_{\goth{A}}}_{a}$ for every $a \in W$, and $\interp{f}$
being the collection of all the $\interp{f}$, we obtain that
$\interp{f} = \Lambda$.

§.§ Uniqueness of the Normal Form

An crucial feature of the normal form for the completeness is its uniqueness:

The normal form is unique.

Let $f$ and $g$ be two diagrams in normal form (with respectively
$\lambda$ and $\mu$ as inner block), such that $\interp{f}=\interp{g}$
(the naming of the worlds is taken to be the same in both diagrams,
and is the same as in the previous proof). By the definition of
$\interp{.}$, we have $\interp{f}_a=\interp{g}_a$ for every $a\in W$.
We hence have $\interp{\iso^{-1}_{\goth{B}}}_{a} \circ
\interp{\lambda}_a \circ \interp{\iso_{\goth{A}}}_{a} =
\interp{\iso^{-1}_{\goth{B}}}_{a} \circ \interp{\mu}_a \circ
\interp{\iso_{\goth{A}}}_{a}$.

Denoting $e_i^{\goth A}$ (resp. $e_i^{\goth B}$) the $i$-th element of
the basis of $\goth A$ (resp. $\goth B$), we have:

\[\lambda_{ij} = e_j^{\goth
B\dagger}\interp{\iso^{-1}_{\goth{B}}}_{a_{ij}} \circ
\interp{\lambda}_{a_{ij}} \circ
\interp{\iso_{\goth{A}}}_{a_{ij}}e_i^{\goth A} = e_j^{\goth
B\dagger}\interp{\iso^{-1}_{\goth{B}}}_{a_{ij}} \circ
\interp{\mu}_{a_{ij}} \circ
\interp{\iso_{\goth{A}}}_{a_{ij}}e_i^{\goth A} = \mu_{ij}\]
Hence, all
coefficients in the scalars of $f$ and $g$ are the same. Since the
structure is otherwise the same for $f$ and $g$, they are the same

§.§ Completeness

We can now use this normal form to show that our equational theory is
complete for arbitrary morphisms. To do so, we need to show that all
the generators can be put in normal form, and then that any
composition of morphisms in normal form can be put in normal form.
To do so we will first derive a few lemmas:

\equiv\tikzfigc{05}$ when $s_1\cap s_4 = s_2\cap s_3 = \varnothing$

Single-colored isos distribute over the contraction:

\[\tikzfig{isos-distribute-over-c}\]

Scalars distribute over single-colored isos:

\[\tikzfig{scalars-distribute-over-isos}\]

The result is obvious for $\one$ and $\emptysquare$. For $A\oplus B$:

\begin{align*}
\tikzfigc{00}
\equiv\tikzfigc{02}
\equiv\tikzfigc{03}
\end{align*}

For $A\otimes B$:

\begin{align*}
\tikzfigc{00}
\equiv\tikzfigc{02}
\equiv\tikzfigc{03}\\
\equiv\tikzfigc{05}
\end{align*}

\begin{align*}
\tikzfigc{00}
\equiv\tikzfigc{03}
\end{align*}

with $s_0\cap s_1=\varnothing$.

\begin{align*}
\equiv\tikzfigc{05}
\end{align*}

with $\nu_{ij} = \sum_k\lambda_{ik}\mu_{kj}$.

Suppose the worlds of $\lambda$ are the $\{a_{ij}\}_{ij}$ and that of
$\mu$ are the $\{b_{k\ell}\}_{k\ell}$. We count inputs/outputs
starting at $1$, hence $p\geq1$ in the following.

The $p$-th output of $\lambda$ has world set $\{a_{ip}\}_{i}$ and the
$p$-th input of $\mu$ has world set $\{b_{p\ell}\}_{\ell}$. When
composing the two in sequence, in the first step, each singleton
world $\{a_{ij}\}$ (resp. $\{b_{k\ell}\}$) is mapped to
(resp. $\{(a_{ij},b_{k\ell})\}_{ij}$), and unions of singleton worlds
to unions of the worlds each is mapped to.

In particular, the $p$-th output of $\lambda$ now has world set
$\{(a_{ip},b_{k\ell})\}_{ik\ell}$, and the $p$-th input of $\mu$ now
has world set $\{(a_{ij},b_{p\ell})\}_{ij\ell}$. After composition,
the two sets have to match. They become
$\{(a_{ip},b_{p\ell})\}_{i\ell}$, and all pairs of world that were
removed from these two sets are removed globally.

We do so for all $p\geq 1$, which means all worlds in:

\[\{(a_{ij},b_{k\ell})\mid j\geq1, j\neq
k\}\cup\{(a_{ij},b_{k\ell})\mid k\geq1, j\neq
k\}=\{(a_{ij},b_{k\ell})\mid j\neq k\}\]
are removed. Crucially, this
means that the world sets $\{a_{i0}\}_{i}$ and $\{b_{0\ell}\}_{\ell}$
are mapped to the same world set:

\begin{align*}
\{a_{i0}\}_{i}
\mapsto &~\{(a_{i0},b_{k\ell})\}_{ik\ell}\setminus\{(a_{ij},b_{k\ell})\mid j\neq k\}\\
=\{(a_{ij},b_{0\ell})\}_{ij\ell}\setminus\{(a_{ij},b_{k\ell})\mid j\neq k\}
\mapsfrom\{b_{0\ell}\}_{\ell}
\end{align*}

Hence, the sequencial composition of $\lambda$ and $\mu$ becomes:

\begin{align*}
\tikzfigc{00}
\equiv\tikzfigc{01}
\equiv\tikzfigc{02}\\
\equiv\tikzfigc{03}
\equiv\tikzfigc{04}
\end{align*}

with $\nu_{ij} = \sum_k\lambda_{ik}\mu_{kj}$.

For any "matrix block" $\lambda$, there exists a "matrix block" $\nu$ such that:

\begin{align*}
\tikzfigc{00}
\equiv\tikzfigc{05}
\end{align*}

We show the result through several steps. First we show the following:

\begin{align}
\label{eq:dangling-scalar-through-unitor}
\tikzfigc{00}
\equiv\tikzfigc{01}
\equiv\tikzfigc{02}
\equiv\tikzfigc{03}
\end{align}

Then, we deal with the top dangling scalars and prove that:

\begin{align}
\label{eq:dangling-scalar-through-iso-boxempty}
\tikzfigc{00}
\equiv\tikzfigc{04}
\end{align}

\begin{align*}
\tikzfigc{00}
\underset{(\ref{eq:dangling-scalar-through-unitor})}\equiv\tikzfigc{02}\\
\equiv\tikzfigc{04}
\end{align*}

This case generalizes to any bottom wire the dangling scalar is
applied on. When we have several of them, we can make them go through
the top part in turns, then aggregate them under a single "matrix
block" using Lemma

We then derive the following equation:

\begin{align}
\label{eq:NF-tensor-lemma}
\tikzfigc{01}
\equiv\tikzfigc{02}
\equiv\tikzfigc{03}
\equiv\tikzfigc{04}
\end{align}

We can finally prove the lemma:

\begin{align*}
\equiv\tikzfigc{01}
\underset{(\ref{eq:dangling-scalar-through-iso-boxempty})}\equiv\tikzfigc{02}\\
\underset{(\ref{eq:NF-tensor-lemma})}\equiv\tikzfigc{04}\\
\end{align*}

We can now move on to show that generators can be put in normal form:

The generators can be put in normal form.

\[\tikzfig{NF-tensor-simp-gen}\]

\[\tikzfig{NF-plus-simp-gen}\]

\[\tikzfig{NF-swap-boxempty}\]

with $\sigma$ a simple permutation of wires.

\[\tikzfig{NF-contraction-gen}\]

\[\tikzfig{NF-scalar}\]

\[\tikzfig{NF-unit}\]

\[\tikzfig{NF-cap}\]
The upside-down versions of the generators
are provided in exactly the same way (but upside-down).

Compositions of diagrams in normal form can be put in normal form.

It is then possible to show that compositions of diagrams in normal form can be put in normal form:

In the case of sequential composition:

\begin{align*}
\tikzfigc{00}
\underset{\ref{lem:isos-are-isos}}\equiv\tikzfigc{01}
\underset{\ref{lem:sequential-composition-of-matrix-blocks}}\equiv\tikzfigc{02}
\end{align*}

In the case of parallel composition:

\begin{align*}
\tikzfigc{00}
\underset{\ref{lem:isos-are-isos}}\equiv\tikzfigc{01}
\underset{\ref{lem:isos-are-isos}}\equiv\tikzfigc{02}
\underset{\ref{lem:matrix-block-through-iso-boxempty}}\equiv\tikzfigc{03}\\
\underset{\ref{lem:matrix-block-through-iso-boxempty}}\equiv\tikzfigc{04}
\underset{\substack{\ref{lem:isos-are-isos}\\\ref{lem:sequential-composition-of-matrix-blocks}}}\equiv\tikzfigc{05}
\equiv\tikzfigc{06}
\underset{\ref{lem:sequential-composition-of-matrix-blocks}}\equiv\tikzfigc{07}
\end{align*}

where $\overline{\goth A\boxempty \goth C}$ represents the canonical choice of composition with $\boxempty$ (and similarly for $\overline{\goth B\boxempty \goth D}$).

This finally allows us to prove the completeness theorem claimed above:

The right-to-left direction of the equivalence can be directly
checked by verifying that all the axioms preserve the

Let $f_1$ and $f_2$ be two morphisms such that
$\interp{f_1}=\interp{f_2}$. Both morphisms can be put in
normal form, resp. $f_1^{NF}$ and $f_2^{NF}$, with $f_i\equiv
f_i^{NF}$ and thus $\interp{f_i^{NF}}=\interp{f_i}$. By
uniqueness of the normal form, and since
$\interp{f_1^{NF}}=\interp{f_2^{NF}}$, we get $f_1^{NF}\equiv
f_2^{NF}$, which ends the proof that $f_1\equiv f_2$.


# A Parallel Integrated Computational-Statistical Platform
for Turbulent Transport Phenomena

Ali Akhavan-Safaei  Mohsen Zayernouri<EMAIL_ADDRESS>Department of Mechanical
Engineering, Michigan State University, East Lansing, MI 48824, USA
Department of Computational Mathematics, Science and Engineering, Michigan
State University, East Lansing, MI 48824, USA  Department of Statistics and
Probability, Michigan State University, East Lansing, MI 48824, USA

###### Abstract

In this paper, we present an open-source, automated, and multi-faceted
computational-statistical platform to obtain synthetic homogeneous isotropic
turbulent flow and passive scalar transport. A parallel implementation of the
well-known pseudo-spectral method in addition to the comprehensive record of
the statistical and small-scale quantities of the turbulent transport are
offered for executing on distributed memory CPU-based supercomputers. The
user-friendly workflow and easy-to-run design of the developed package is
disclosed through an extensive and step-by-step example. The resulting low-
and high-order statistical records vividly verify well-established and fully-
developed turbulent state as well as the seamless statistical balance of
conservation laws. Post-processing tools provided in this platform would let
the user to readily construct multiple important transport quantities from the
primitive turbulent fields.

###### keywords:

Turbulent Transport , Passive Scalars , Statistical Analysis , Pseudo-spectral
Method , High-Performance Computing

Code Metadata Version | v1.0
---|---
Link to code/repository | github.com/FMATH-Group/PSc_HIT3D
Legal Code Licence | General Public License (GPLv3)
Code versioning system used | git
Software code languages | Python, MPI
Requirements & dependencies | NumPy, SciPy, MPI4Py
Support email for questions |<EMAIL_ADDRESS>

##  1 Motivation and Significance

Developing an open-source, sustainable, portably parallel, and integrated
computational-statistical platform with high-order spatial and temporal
accuracy provides a useful academic ground for better understanding of complex
standard to anomalous turbulent transport across a multitude of scales [1].
Moreover, from the educational point of view, developing such a user-friendly
scientific software will essentially fill the existing training gap in the
subjective trinity, i.e., fluid mechanics, computational fluid dynamics, and
turbulent transport; hence, leading to a more cohesive ramp up in training the
future generation of researchers in a variety of academic-to-industrial
disciplines.

Direct numerical simulation (DNS) of turbulent transport as a rigorous
scientific tool is supposed to fully resolve the smallest scales of the motion
resulted from the fluctuating fields in spatial domain while maintaining a
high-order temporal accuracy as turbulence evolves in time [2]. Among the
current open-source computational platforms, Nektar++ [3, 4] (the
spectral/$hp$ element method flow software), HERCULIS [5] and Xcompact3D [6]
(the high-order finite difference flow solvers), GRINS [7] (the adaptive mesh
refinement finite element method software), spectralDNS [8] (the spectral
method computational package for DNS), and OpenFOAM [9] are the notable
contributions to the DNS of turbulent transport. On the other hand, the random
nature of turbulence requires a thorough statistical analysis on the
fluctuating fields and their gradients so that one can identify when the
realistic and fully-developed turbulent state is obtained in DNS during an
ongoing simulation. This necessitates development of a comprehensive
computational platform that includes computing and recording of such
statistical quantities of turbulent transport as time-series format.

(a)

(b)

Figure 1: Schematic of the architecture of the software. Fig. (a) Illustrates
the pseudo-spectral NS solver to archive fully-developed turbulent state,
statistical records, velocity output and restarting the simulation from output
file. Fig. (b) Shows the pseudo-spectral NS and AD solvers to reach fully-
developed turbulent scalar state.

In the current work, our goal is to offer an extensible open-source
computational platform that carries out the high-fidelity simulations of
homogeneous isotropic turbulence (HIT) for an incompressible flow and also
obtains the transport of a passive scalar (temperature or concentration of
species) in such flow [10, 11] while it keeps track of statistical quantities
of turbulent transport. Here, we numerically solve the incompressible Navier-
Stokes (NS) equations in addition to the advection-diffusion (AD) to
sufficiently resolve the fluid velocity and passive scalar fields,
respectively. Spatial homogeneity of fluctuating fields makes this problem
well-suited for pseudo-spectral implementation of the NS and AD equations
based on the Fourier collocation discretization as employed in our work. This
computational platform is based upon programming in PYTHON and leveraging MPI
library for parallel implementation.

The rest of this paper is organized as follows: in section 2, we describe the
details and capabilities of the developed platform as a scientific software
and we point out the theoretical backgrounds briefly. Furthermore, in section
3, we go over a comprehensive example illustrating the results of a fully-
developed turbulent flow and passive scalar field with proper statistical
testing and verification. In section 4, we outline broader impacts of the
current work onto the research in turbulent transport. Finally, in section 5,
we conclude this paper with a summary and conclusions.

##  2 Software Description

###  2.1 Governing Equations

The incompressible HIT considered in the present software is governed by the
NS equations

$\displaystyle\frac{\partial\bm{U}}{\partial t}+\bm{U}\cdot\nabla\bm{U}$
$\displaystyle=-\nabla p+\nu\,\Delta\bm{U}+\mathcal{A}\,\bm{U},$ (1)

subject to the continuity $\nabla\cdot\bm{U}=0$. In (1),
$\bm{U}=(U_{1},U_{2},U_{3})$ and $p$ are the instantaneous velocity and
modified pressure (pressure divided by the constant density of fluid) fields
in the Cartesian coordinate system $\bm{x}=(x_{1},x_{2},x_{3})$, respectively.
Moreover, $\nu$ denotes the dynamic viscosity of the Newtonian fluid, and
$\mathcal{A}$ is a dynamically evaluated coefficient corresponding to the
artificial forcing scheme we employ in order to obtain statistically
stationary and fully-turbulent state. From the Reynolds decomposition of
instantaneous velocity field,
$\bm{U}(\bm{x},t)=\left\langle\bm{U}(\bm{x},t)\right\rangle+\bm{u}(\bm{x},t)$,
where $\langle\cdot\rangle$ represents the ensemble-averaging operator, and
$\bm{u}(\bm{x},t)$ denotes the fluctuating part of the velocity field. In HIT,
$\langle\bm{U}(\bm{x},t)\rangle=0$; therefore, the instantaneous velocity
field equals the fluctuating part that is governed by (1). Introducing a
passive scalar $\Phi(\bm{x},t)$ transported in the considered fully-developed
HIT flow, the AD equation governing the passive scalar concentration may be
formulated as

$\displaystyle\frac{\partial\Phi}{\partial t}+\bm{u}\cdot\nabla\Phi$
$\displaystyle=\mathcal{D}\,\Delta\Phi,$ (2)

where $\mathcal{D}$ denotes the diffusivity of passive scalar. Applying the
Reynolds decomposition on the total passive scalar,
$\Phi=\langle\Phi\rangle+\phi$, and $\phi$ is the fluctuating part of the
scalar concentration. Considering a uniform mean gradient for the passive
scalar as $\nabla\langle\Phi\rangle=(0,\beta,0)$, where $\beta$ is a constant,
the AD equation in (2) is rewritten as

$\displaystyle\frac{\partial\phi}{\partial t}+\bm{u}\cdot\nabla\phi$
$\displaystyle=-\beta\,u_{2}+\mathcal{D}\,\Delta\phi.$ (3)

###  2.2 Fourier Pseudo-Spectral Method

Here, we consider spatial homogeneity for the fluctuating velocity and scalar
concentration, which allows periodic boundary conditions for these fluctuating
fields as

$\displaystyle\bm{u}(\bm{x}+\mathcal{L}\,\mathrm{e}_{i},t)=\bm{u}(\bm{x},t),\quad\phi(\bm{x}+\mathcal{L}\,\mathrm{e}_{i},t)=\phi(\bm{x},t),$
(4)

where $\mathrm{e}_{i},\>_{i=1,2,3}$, is the unit vector for the $i$-th
direction of the Cartesian coordinate, and $\mathcal{L}$ is the periodicity
length that defines the spatial domain as $\bm{\Omega}=[0,\mathcal{L}]^{3}$.
Discretizing $\bm{\Omega}$ using a uniform three-dimensional grid returns
$N^{3}$ grid points with grid spacing along each direction as $\Delta
x=\mathcal{L}/N$. Transforming this discretization into spectral domain let us
have a standard pseudo-spectral representation of the governing equations (1)
and (3). Subsequently, $\bm{k}=(k_{1},k_{2},k_{3})$ represents the coordinate
system in the spectral space and using Fourier collocation method the
discretized representation of $\bm{k}$ would be $k_{i}=(-N/2+1,\dots,N/2)$,
$i=1,2,3$. Accordingly, the discrete Fourier transform of any field variable
such as $\phi(\bm{x},t)$ is written as

$\displaystyle\phi(\bm{x},t)=\frac{1}{N^{3}}\sum_{\bm{k}}\hat{\phi}_{\bm{k}}(t)\,e^{\mathfrak{i}\bm{k}\cdot\bm{x}},$
(5)

where $\mathfrak{i}=\sqrt{-1}$, and $e^{\mathfrak{i}\bm{k}\cdot\bm{x}}$ are
the Fourier basis functions. Subsequently, the Fourier coefficients associated
with $\bm{k}$ are represented as
$\hat{\phi}_{\bm{k}}(t)=\sum_{\bm{x}}\phi(\bm{x},t)\,e^{-\mathfrak{i}\bm{k}\cdot\bm{x}}$.
Standard pseudo-spectral formulation of the NS equations based upon the
Fourier collocation method is obtained after taking the Fourier transform of
(1),

$\displaystyle\frac{d\hat{\bm{u}}_{\bm{k}}}{dt}+(\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("\bm{u}\cdot\nabla\bm{u}")]{\kern
0.1pt\mathchar 866\relax\kern
0.1pt}{\rule{0.0pt}{505.89pt}}}{}}{2.4ex}}\stackon[-6.9pt]{\bm{u}\cdot\nabla\bm{u}}{\tmpbox})$
$\displaystyle=-\mathfrak{i}\bm{k}\>\hat{p}_{\bm{k}}-\nu|\bm{k}|^{2}\hat{\bm{u}}_{\bm{k}}+\mathcal{A}\hat{\bm{u}}_{\bm{k}};$
(6) $\displaystyle\mathfrak{i}\bm{k}\cdot\hat{\bm{u}}_{\bm{k}}=0,$

By taking the divergence of momentum equation in (6) and applying the
continuity, modified pressure is explicitly represented in terms of the
velocity field. Considering that ${\bm{k}}\cdot{\bm{k}}=|\bm{k}|^{2}$, one can
derive
$\hat{p}_{\bm{k}}=\mathfrak{i}\bm{k}\cdot(\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("\bm{u}\cdot\nabla\bm{u}")]{\kern
0.1pt\mathchar 866\relax\kern
0.1pt}{\rule{0.0pt}{505.89pt}}}{}}{2.4ex}}\stackon[-6.9pt]{\bm{u}\cdot\nabla\bm{u}}{\tmpbox})_{\bm{k}}/|\bm{k}|^{2}$;
hence, (6) may be reformulated as

Similarly, the pseudo-spectral representation of the AD equation for passive
scalar (3), is written as

$\displaystyle\frac{d\hat{\phi}_{\bm{k}}}{dt}+(\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("\bm{u}\cdot\nabla\phi")]{\kern
0.1pt\mathchar 866\relax\kern
0.1pt}{\rule{0.0pt}{505.89pt}}}{}}{2.4ex}}\stackon[-6.9pt]{\bm{u}\cdot\nabla\phi}{\tmpbox})$
$\displaystyle=-\beta\,\hat{u}_{k_{2}}-\mathcal{D}|\bm{k}|^{2}\hat{\phi}_{\bm{k}}.$
(7)

Employing the fourth-order Runge-Kutta (RK4) scheme, the time-stepping for
both NS and AD equations is explicitly done since the nonlinear (advetctive)
terms are evaluated in the physical space and then transformed into the
spectral space.

Table 1: Statistical characteristics of turbulent flow to be recorded from the
NS solver as time-series within user-defied time-intervals.

TKE | Dissipation | Kolmogorov length-scale | Taylor-scale Reynolds | Large-eddy turnover time
---|---|---|---|---
$K=\frac{1}{2}\langle\bm{u}\cdot\bm{u}\rangle$ | $\varepsilon=2\,\nu\,\langle\bm{S}:\bm{S}\rangle$ | $\eta=\left(\nu^{3}/\varepsilon\right)^{1/4}$ | $Re_{\lambda}=\lambda\,u_{rms}/\nu$ | $T_{e}=l_{o}/u_{rms}$

  * $\ast$

$\bm{S}=\frac{1}{2}\left(\nabla\bm{u}+\nabla\bm{u}^{\mathrm{T}}\right),\qquad
u_{rms}=\sqrt{2K/3},\qquad\lambda=u_{rms}\sqrt{15\nu/\varepsilon},\qquad
l_{o}=u_{rms}^{3}/\varepsilon$.

###  2.3 Software Architecture

(a)

(b)

(c)

Figure 2: (a) Snapshot of fully-developed turbulent velocity field, $u_{1}$
component. (b) Time-averaged TKE spectrum. (c) Time-series of $Re_{\lambda}$
(red dashed line), and VGT skewness factors, $\mathcal{S}_{u_{1,1}}$,
$\mathcal{S}_{u_{2,2}}$, $\mathcal{S}_{u_{3,3}}$ .

The structure of the prepared software is schematically illustrated in Figure
1. According to Figure 1, a user starts from a pre-processing step, where the
isotropic and random velocity initial condition (IC) is constructed based on a
prescribed spectrum for turbulent kinetic energy (TKE). The procedure is the
straightforward implementation of the well-known work by Rogallo to generate
divergence-free isotropic velocity state [12]. According to Lamorgese et al.
[13], the initial TKE spectrum is chosen to be

$\displaystyle
E(\kappa,0)=\frac{u_{rms}^{2}}{k_{F}}\times\begin{cases}(\kappa/k_{F})^{2},&\text{if}\
\ \ \kappa\leq k_{F},\\\ (\kappa/k_{F})^{-5/3},&\text{if}\ \ \
\kappa>k_{F}.\end{cases}$ (10)

where $\kappa$ represents the wavenumber associated with spherical shells,
$k_{F}$ denotes the maximum wavenumber of TKE shell we apply artificial
forcing to, and $u_{rms}$ specifies the initial root-mean-square (rms)
intensity of velocity fluctuations. In construction of velocity IC, $u_{rms}$
is set to be unity while $k_{F}$ and the number of Fourier collocation points,
$N$ are taken as input parameters. In a UNIX/LINUX environment, these inputs
are taken as arguments in the execution commandline that are imported through
sys library in PYTHON. Once the velocity IC is obtained, it is partitioned
into $N_{p}$ slabs according to the slab decomposition method. We completely
adopted the implementation of Mortensen and Langtangen [8] for domain
decomposition in addition to the parallel implementation of forward and
inverse three-dimensional fast Fourier transform (FFT) in PYTHON programming
language. Here, the MPI communications depend on the mpi4py library [14, 15,
16].

Having the partitioned velocity IC prepared in the pre-processing step, it is
fed to the main body of the software where the initial velocity field might be
magnified by a user-defined input argument so that a target TKE is considered
for the simulation. The viscosity of fluid, $\nu$, is also is taken as another
user-defined input argument. Next, the magnified velocity field is passed into
the solver where $\hat{\bm{u}}_{\bm{k}}$ and
$(\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("\nabla\bm{u}")]{\kern
0.1pt\mathchar 866\relax\kern
0.1pt}{\rule{0.0pt}{505.89pt}}}{}}{2.4ex}}\stackon[-6.9pt]{\nabla\bm{u}}{\tmpbox})_{\bm{k}}=\mathfrak{i}\bm{k}\hat{\bm{u}}_{\bm{k}}$
are separately transformed back into the physical space so that
$\bm{u}\cdot\nabla\bm{u}$ is simply computed and then transformed into Fourier
space. The aliasing error that appears due to this procedure is removed by
phase-shifting and truncation according to $2\sqrt{2}N/3$ as the maximum
wavenumber [17]. Afterwards, all of the terms in equation (2.2) are directly
evaluated in every stage of RK4 time-integration; however, the last term in
the right-hand side of (2.2) is only evaluated after the last stage during the
artificial forcing by keeping the energy of the low wavenumbers constant,
which is associated with the sphere of $0<|\bm{k}|\leq k_{F}$. In this
procedure, $\mathcal{A}$ is computed in a way that the the dissipated energy
of turbulent motion is injected to the large-scales. This scheme prevents the
flow to undergo decay process before the realistic and fully-developed
turbulent state is achieved, nevertheless, the artificial forcing scheme could
be turned off through a user-defined input argument if one seek to obtain
decaying HIT data. The forcing coefficient could be determined either
deterministically [18] or stochastically [19, 20] and both of these methods
are supported in the software as would be specified as an input option.
Moreover, regarding the stable time-integration, the Courant-Friedrichs-Lewy
(CFL) number is dynamically checked through a user-defined time-frequency.
According to Eswaran and Pope’s work [19], CFL for this problem is
demonstrated as

$\displaystyle\mathrm{CFL}:=\frac{\Delta t}{\Delta
x}\max\left(|u_{1}|+|u_{2}|+|u_{3}|\right).$ (11)

In (11), $\Delta x$ is the uniform grid spacing in each direction and $\Delta
t$ is the user-defined constant time-interval used in the RK4 time-stepping.
In practice, CFL is required to be less than unity to ensure a stable time-
integration.

(a)

(b)

Figure 3: (a) Time-records of production over dissipation of scalar variance
(blue dashed line), and the flatness factor for the scalar gradient vector
component along the direction of mean scalar gradient (green solid line). (b)
Snapshot of fully-developed turbulent passive scalar field.

Since the fully-developed turbulent state is characterized by a meticulous
tracking of statistical quantities of the flow, the present software provides
a comprehensive framework for computing and recording the statistical
quantities of turbulent flow. Given the homogeneity of the fluctuating fields,
spatial averaging is employed for computing these records at user-defined time
intervals. These statistical quantities are categorized into turbulent
characteristics of small-scale motion reported in Table 1, and high-order
central moments of diagonal components of velocity gradient tensor (VGT),
$\nabla\bm{u}$. For instance,

$\displaystyle\mathcal{S}_{u_{1,1}}=$
$\displaystyle\left.\left\langle\Big{(}\frac{\partial u_{1}}{\partial
x_{1}}\Big{)}^{3}\right\rangle\middle/\left\langle\Big{(}\frac{\partial
u_{1}}{\partial x_{1}}\Big{)}^{2}\right\rangle^{3/2}\right.,$ (12)
$\displaystyle\mathcal{K}_{u_{1,1}}=$
$\displaystyle\left.\left\langle\Big{(}\frac{\partial u_{1}}{\partial
x_{1}}\Big{)}^{4}\right\rangle\middle/\left\langle\Big{(}\frac{\partial
u_{1}}{\partial x_{1}}\Big{)}^{2}\right\rangle^{2}\right.,$ (13)

where $\mathcal{S}_{u_{1,1}}$ and $\mathcal{K}_{u_{1,1}}$ indicate the
skewness factor and flatness factor (or kurtosis) associated with the first
diagonal component of VGT, $u_{1,1}=\partial u_{1}/\partial x_{1}$,
respectively. Fully-turbulent flow state would be identified when the time-
series of these records reach to a statistically stationary state after long
enough time-integration, i.e. approximately 10 to 15 large-eddy turnover times
(see Table 1). The parallel implementation for computing and collecting these
statistical quantities and later recording them in time as different time
series which were performed by point-to-point and collective MPI directives.

Furthermore, the velocity and pressure fields might be written as output files
stored in directories named Out$\\_\ast$ based on a user-defined time-interval
that might be useful for any post-processing after the flow reaches to fully-
developed turbulent state. A “restart from file” capability is also designated
so that once the statistical record is written out on file, the latest state
of velocity field and its related time-integration information are also output
on files, which are stored in a directory named Restart. Starting a simulation
from either a prescribed IC or restarting it to continue an ongoing simulation
that was stopped is specified by a user-defined input argument. All the
parallel I/O to store the velocity and pressure fields is done by employing
scipy.io library and using loadmat() and savemat() routines for the partition
of data that is resolved inside each MPI process in the compressed format and
with the machine precision accuracy. This is a fast and efficient I/O while it
maintains the simplicity for user for any post-processing step such as time-
averaging on the statistically stationary turbulent data.

According to Figure 1, once the fully-turbulent velocity state is achieved,
the user might be able to use a restart or output instance as the velocity IC
to introduce a passive scalar transport with a directional constant mean-
gradient as described in (3) while the fluctuating concentration is assumed to
be zero, $\phi_{0}(\bm{x})=0$. Here, the goal is to resolve the fluctuating
scalar concentration field, equation (7), transported on the fully-turbulent
incompressible flow for long enough time-span so that the fully-developed and
realistic turbulent state for the passive scalar is obtained. Subsequently,
similar procedure as what described based on the schematic Figure 1 is
followed while the pseudo-spectral AD solver is fed by the resolved velocity
field from the NS solver. The diffusivity of passive scalar, $\mathcal{D}$, is
specified by a user-defined input argument for Schmidt number,
$Sc=\nu/\mathcal{D}$. Accordingly, similar to the NS solver, the advective
scalar flux,
$(\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("\bm{u}\cdot\nabla\phi")]{\kern
0.1pt\mathchar 866\relax\kern
0.1pt}{\rule{0.0pt}{505.89pt}}}{}}{2.4ex}}\stackon[-6.9pt]{\bm{u}\cdot\nabla\phi}{\tmpbox})_{\bm{k}}$,
is computed in the physical space by the inverse FFT of
$\hat{\bm{u}}_{\bm{k}}$ and
$(\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("\nabla\phi")]{\kern
0.1pt\mathchar 866\relax\kern
0.1pt}{\rule{0.0pt}{505.89pt}}}{}}{2.4ex}}\stackon[-6.9pt]{\nabla\phi}{\tmpbox})_{\bm{k}}$
and forward FFT computation of $\bm{u}\cdot\nabla\phi$. Similar dealiasing
procedure as described for NS solver is employed in pseudo-spectral AD solver
and the RK4 time-integration scheme is utilized to numerically perform
explicit time-stepping.

In homogeneous scalar turbulence, time evolution of the scalar variance
$\langle\phi^{2}\rangle$ is governed by

$\displaystyle\frac{d}{dt}\langle\phi^{2}\rangle=$
$\displaystyle-2\,\langle\bm{q}\rangle\cdot\nabla\langle\Phi\rangle-\langle\chi\rangle,$
(14)

where $\bm{q}=\phi\,\bm{u}$ denotes the scalar flux vector, and the scalar
dissipation is defined as $\chi=2\,\mathcal{D}\,\nabla\phi\cdot\nabla\phi$
[21]. According to (3), the first term in the right-hand side of (14) is
simplified to $-2\,\beta\,\langle\phi\,u_{2}\rangle$ that denotes the scalar
variance production (by uniform mean scalar gradient, $\beta$). The present
software is capable of computing and recording of rate of scalar variance in
addition to the production and dissipation terms. This is useful in terms of
the checking if the balance of both sides of equation (14) holds throughout a
simulation so that one ensures that the implementation of the solver works
seamlessly. On the other hand, as a measure to evaluate that the statistically
stationary state for the passive scalar is achieved is to check if
$-2\,\beta\,\langle\phi\,u_{2}\rangle/\langle\chi\rangle\sim 1$ throughout the
simulation. Moreover, recording the skewness and flatness factors for the
components of fluctuating scalar gradient vector (e.g.,
$\mathcal{S}_{\phi_{,2}}$ and $\mathcal{K}_{\phi_{,2}}$ for
$\phi_{,2}=\partial\phi/\partial x_{2}$ similar to 12 and 13 for VGT) is
another statistical indicator measure for fully-developed turbulent passive
scalar state. Therefore, in the current computational platform, the user would
be able to recognize the statistically stationary state through monitoring the
explained time-series data that is written out according to the user-defined
time-interval as a software input.

The field data output and restart capability for the AD solver is designated
similar to the described strategy for the NS solver so that the user would be
able to resume an interrupted/stopped simulation and use the output field data
for desired applications or post-processing.

In the following section, we present a comprehensive example that step-by-step
walks through using the present software.

##  3 Illustrative Example

This comprehensive example is mainly consisted of construction of isotropic
velocity IC, obtaining well-resolved fully-turbulent velocity field, and
simulating well-resolved passive scalar turbulence with imposed mean scalar
gradient.

###  3.1 IC construction and DNS of the HIT

Table 2: Input arguments for the PScHIT.py and the specified values for the example case. The order of the arguments in the execution command-line are as listed here. Input Argument |  | Value
---|---|---
t_end |  | 40
Output_frequency |  | 1000
Stats_frequency |  | 100
TKE_magnification |  | 6.0
$\nu$ |  | 0.0008
$k_{F}$ |  | 2
forcing_type |  | deterministic
$N$ |  | 520
$\Delta t$ |  | 0.0005
If_Restart |  | 0 or 1

According to the descriptions in section 2.3, the isotropic and divergence-
free velocity IC is constructed based upon a prescribed energy spectrum given
in (10). Considering the periodicity length $\mathcal{L}=2\pi$, this pre-
processing step is done through serial execution of Gen_IC.py script that
takes the following input arguments, respectively: $N$ (spatial resolution
along each direction), $k_{F}$ (forcing wavenumber), and $N_{p}$ (number of
slab partitions). We need to emphasize that $N_{p}$ must be chosen in a way
that $N$ be a multiple of $N_{p}$. The resulting velocity field is located in
a directory named IC, where $N_{p}$ number of .mat velocity files are stored.
In this example, we take $N=520$, $k_{F}=2$, and $N_{p}=40$. All the
components of velocity IC in addition to the VGT components have Gaussian
distribution. This velocity IC is being passed into the NS solver written in
PScHIT.py script that takes the following input arguments given in Table 2.
Here, Output_frequency and Stats_frequency are multiplied by the specified
$\Delta t$. Moreover, If_Restart argument could be either 0 or 1, where 0
indicates it is a simulation starting from the constructed IC while 1
specifies resuming a simulation from restart files. In this example, we
perform the simulation for $t/T_{e}\sim 15$ to ensure the fully-turbulent flow
state is achieved and Figure 2 portrays the first component of the velocity
field. Figure 2 shows the radial TKE spectrum averaged over 5 large-eddy
turnover times. Moreover, Figure 2 includes the time-records of the Taylor-
scale Reynolds number and VGT skewness factor for diagonal components computed
and recorded over 40 large-eddy turnover times. This shows that the
statistically stationary state is achieved through the long-time DNS where
$Re_{\lambda}\sim 240$,
$\mathcal{S}_{u_{1,1}}=\mathcal{S}_{u_{2,2}}=\mathcal{S}_{u_{3,3}}\sim-0.55$,
and $\mathcal{K}_{u_{1,1}}=\mathcal{K}_{u_{2,2}}=\mathcal{K}_{u_{3,3}}\sim
6.8$ at fully-turbulent state. The statistical records of VGT clearly show
that the resolved velocity field is isotropic. Finally, $k_{max}\eta>1.45$
ensures that the small scale turbulent motions are well-resolved
($k_{max}=\sqrt{2}N/3$) [22].

###  3.2 DNS of passive scalar transport

Similar to starting the NS solver from a prescribed velocity IC, we take a
fully-turbulent velocity output (velocity state at $t/T_{e}=15$ in section
3.1) and continue the simulation under the artificial forcing while we
introduce a passive scalar field where its fluctuating part is initialized at
zero. The Schmidt number, $Sc$, is specified by user through an input
argument. According to the problem setting for the mean scalar gradient, we
let $\beta=1$ (mean scalar gradient along $x_{2}$ direction). Therefore, for a
passive scalar with $Sc=1$, we aim to obtain the fully-turbulent scalar field.
We need to note that the spatial resolution required for the passive scalars
with $Sc\geq 1$ is defined based on $\eta_{B}=\eta\,Sc^{-1/2}$ [23] and in
this example the spatial resolution for the velocity field is sufficient for a
well-resolved passive scalar. We manage to resolve the passive scalar field
for 25 large-eddy turnover times and the rest of the simulation parameters
remain the same as values reported in the Table 2. Figure 3 shows the records
of scalar variance production over dissipation rate,
$-2\,\langle\phi\,u_{2}\rangle/\langle\chi\rangle$, and the flatness factor
for the scalar gradient along the direction of mean scalar gradient,
$\mathcal{K}_{\phi_{,2}}$. As it is observed, after resolving the passive
scalar field for approximately two large-eddy turnover times,
$-2\,\langle\phi\,u_{2}\rangle/\langle\chi\rangle\sim 1.0$ that means the
equilibrium state for the passive scalar variance is obtained. Moreover, after
approximately three large-eddy turnover times the high-order statistical
moments of the scalar gradient reach to a statically stationary state. For
instance, $\mathcal{S}_{\phi_{,2}}\sim 1.4$, and $\mathcal{K}_{\phi_{,2}}\sim
20.8$ throughout the time-averaging of these statistical moments when
$t/T_{e}\geq 5$. By resolving the passive scalar field through AD equation and
for long enough time after the equilibrium and stationary state, the fully-
turbulent and realistic scalar field is ensured.

##  4 Impact

Current work offers a framework to obtain highly accurate spatio-temporal data
for homogeneous turbulent transport with proper statistical testing from the
recorded quantities. In turbulent transport research, this provides a great
source of high-fidelity data for a variety of innovative contributions. For
instance in large-eddy simulation (LES) of turbulent transport, novel nonlocal
models for the subgrid-scale (SGS) stress and flux terms appearing in the
filtered NS and AD equations heavily depend on the DNS data for such transport
phenomena to compute the exact values of SGS terms to evaluate the model
performance [24]. On the other hand, in the abundance of data and emergence of
the data-driven turbulence models [25], current computational platform would
be a reliable candidate to generate data for training and testing such models
[26, 27, 28, 29, 30, 31, 32]. Moreover, high-Reynolds and well-examined high-
fidelity turbulent transport data from the present platform could be directly
employed in studying the role of coherent turbulent structures and their
effects on the turbulence statistics [33, 1], as well as investigating
topological characteristics of turbulent transport [34, 35], analysis of
extreme events and the internal intermittency [36, 37, 38], and diverse
modeling strategies for the arising closures [39, 40].

##  5 Conclusion

This work presented a computational platform for DNS of homogeneous turbulent
flow and passive scalar transport. This open-source software works based upon
a pseudo-spectral representation of the NS and AD equations on a triply cubic
computational domain with periodic boundary conditions for fluctuating fields.
Using Fourier collocation method, the governing equations are discretized in
space and by employing RK4 scheme the time-stepping is performed. The software
provides a pre-processing step to construct homogeneous and isotropic
divergence-free velocity IC based on prescribed energy spectrum and decompose
it into user-defined partitions. Using artificial forcing scheme, the
dissipated energy is injected to the low wavenumbers so that after long-time
integration, the statistically stationary state is achieved. In order to
examine and identify if the fully-developed turbulent flow is obtained, small-
scale statistical quantities of turbulence in addition to the central moments
of VGT components are computed and recorded. Once the realistic turbulent
velocity field is obtained, the user is able to start resolving a passive
scalar that is transported on the HIT flow while a uniform mean scalar
gradient is imposed. Resolving the scalar fluctuations for long enough time
after reaching to the equilibrium and stationary state provides the fully-
developed turbulent scalar field. Statistical records of scalar gradients in
addition to the records of production and dissipation of scalar variance helps
the user to properly identify when the fully-developed scalar turbulence is
achieved.

## Acknowledgement

This work was financially supported by the MURI/ARO award (W911NF-15-1-0562),
the ARO Young Investigator Program (YIP) award (W911NF-19-1-0444), and the NSF
award (DMS-1923201). The HPC resources and services were provided by the
Institute for Cyber-Enabled Research (ICER) at Michigan State University.

## References

  * [1] A. Akhavan-Safaei, S. H. Seyedi, M. Zayernouri, Anomalous features in internal cylinder flow instabilities subject to uncertain rotational effects, Physics of Fluids 32 (9) (2020) 094107.
  * [2] P. Moin, K. Mahesh, Direct Numerical Simulation: A tool in turbulence research, Annual Review of Fluid Mechanics 30 (1) (1998) 539–578.
  * [3] C. D. Cantwell, D. Moxey, A. Comerford, A. Bolis, G. Rocco, G. Mengaldo, D. De Grazia, S. Yakovlev, J.-E. Lombard, D. Ekelschot, et al., Nektar++: An open-source spectral/hp element framework, Computer physics communications 192 (2015) 205–219.
  * [4] D. Moxey, C. D. Cantwell, Y. Bao, A. Cassinelli, G. Castiglioni, S. Chun, E. Juda, E. Kazemi, K. Lackhove, J. Marcon, et al., Nektar++: Enhancing the capability and application of high-fidelity spectral/hp element methods, Computer Physics Communications 249 (2020) 107110.
  * [5] P. He, A high order finite difference solver for massively parallel simulations of stably stratified turbulent channel flows, Computers & Fluids 127 (2016) 161–173.
  * [6] P. Bartholomew, G. Deskos, R. A. Frantz, F. N. Schuch, E. Lamballais, S. Laizet, Xcompact3D: An open-source framework for solving turbulence problems on a cartesian mesh, SoftwareX 12 (2020) 100550.
  * [7] P. T. Bauman, R. H. Stogner, GRINS: a multiphysics framework based on the libmesh finite element library, SIAM Journal on Scientific Computing 38 (5) (2016) S78–S100.
  * [8] M. Mortensen, H. P. Langtangen, High performance python for direct numerical simulations of turbulent flows, Computer Physics Communications 203 (2016) 53 – 65.
  * [9] T. O. Foundation, Openfoam v8 user guide.
URL https://cfd.direct/openfoam/user-guide

  * [10] Z. Warhaft, Passive scalars in turbulent flows, Annual Review of Fluid Mechanics 32 (1) (2000) 203–240.
  * [11] B. I. Shraiman, E. D. Siggia, Scalar turbulence, Nature 405 (6787) (2000) 639–646.
  * [12] R. Rogallo, Numerical Experiments in Homogeneous Turbulence, NASA technical memorandum, National Aeronautics and Space Administration, 1981.
  * [13] A. G. Lamorgese, D. A. Caughey, S. B. Pope, Direct numerical simulation of homogeneous turbulence with hyperviscosity, Physics of Fluids 17 (1) (2005) 015106\.
  * [14] L. Dalcín, R. Paz, M. Storti, MPI for python, Journal of Parallel and Distributed Computing 65 (9) (2005) 1108–1115.
  * [15] L. Dalcín, R. Paz, M. Storti, J. D’Elía, MPI for python: Performance improvements and MPI-2 extensions, Journal of Parallel and Distributed Computing 68 (5) (2008) 655–662.
  * [16] L. D. Dalcin, R. R. Paz, P. A. Kler, A. Cosimo, Parallel distributed computing using python, Advances in Water Resources 34 (9) (2011) 1124–1139.
  * [17] G. S. Patterson, S. A. Orszag, Spectral calculations of isotropic turbulence: Efficient removal of aliasing interactions, The Physics of Fluids 14 (11) (1971) 2538–2541.
  * [18] N. P. Sullivan, S. Mahalingam, R. M. Kerr, Deterministic forcing of homogeneous, isotropic turbulence, Physics of Fluids 6 (4) (1994) 1612–1614.
  * [19] V. Eswaran, S. Pope, An examination of forcing in direct numerical simulations of turbulence, Computers & Fluids 16 (3) (1988) 257 – 278.
  * [20] K. Alvelius, Random forcing of three-dimensional homogeneous turbulence, Physics of Fluids 11 (7) (1999) 1880–1889.
  * [21] A. Monin, A. Yaglom, Statistical Fluid Mechanics, Volume II: Mechanics of Turbulence, Dover Books on Physics, Dover Publications, 2013.
  * [22] M. R. Overholt, S. B. Pope, Direct numerical simulation of a passive scalar with imposed mean gradient in isotropic turbulence, Physics of Fluids 8 (11) (1996) 3128–3148.
  * [23] G. K. Batchelor, Small-scale variation of convected quantities like temperature in turbulent fluid part 1. general discussion and the case of small conductivity, Journal of Fluid Mechanics 5 (1) (1959) 113–133.
  * [24] M. Samiee, A. Akhavan-Safaei, M. Zayernouri, A fractional subgrid-scale model for turbulent flows: Theoretical formulation and a priori study, Physics of Fluids 32 (5) (2020) 055102.
  * [25] K. Duraisamy, G. Iaccarino, H. Xiao, Turbulence modeling in the age of data, Annual Review of Fluid Mechanics 51 (2019) 357–377.
  * [26] Z. Wang, K. Luo, D. Li, J. Tan, J. Fan, Investigations of data-driven closure for subgrid-scale stress in large-eddy simulation, Physics of Fluids 30 (12) (2018) 125101.
  * [27] A. Beck, D. Flad, C.-D. Munz, Deep neural networks for data-driven LES closure models, Journal of Computational Physics 398 (2019) 108910.
  * [28] G. D. Portwood, B. T. Nadiga, J. A. Saenz, D. Livescu, Interpreting neural network models of residual scalar flux, Journal of Fluid Mechanics 907 (2021) A23.
  * [29] A. T. Mohan, D. Tretiak, M. Chertkov, D. Livescu, Spatio-temporal deep learning models of 3D turbulence with physics informed diagnostics, Journal of Turbulence 0 (0) (2020) 1–41.
  * [30] A. T. Mohan, N. Lubbers, D. Livescu, M. Chertkov, Embedding hard physical constraints in neural network coarse-graining of 3D turbulence, arXiv preprint arXiv:2002.00021.
  * [31] S. Beetham, J. Capecelatro, Formulating turbulence closures using sparse regression with embedded form invariance, Phys. Rev. Fluids 5 (2020) 084611.
  * [32] J. Sirignano, J. F. MacArt, J. B. Freund, DPM: A deep learning PDE augmentation method with application to large-eddy simulation, Journal of Computational Physics (2020) 109811.
  * [33] M. Zayernouri, M. Metzger, Coherent features in the sensitivity field of a planar mixing layer, Physics of Fluids 23 (2) (2011) 025105.
  * [34] M. Gonzalez, Kinematic properties of passive scalar gradient predicted by a stochastic lagrangian model, Physics of Fluids 21 (5) (2009) 055104.
  * [35] C. Meneveau, Lagrangian dynamics and models of the velocity gradient tensor in turbulent flows, Annual Review of Fluid Mechanics 43 (2011) 219–245.
  * [36] P. Yeung, X. Zhai, K. R. Sreenivasan, Extreme events in computational turbulence, Proceedings of the National Academy of Sciences 112 (41) (2015) 12633–12638.
  * [37] M. Farazmand, T. P. Sapsis, A variational approach to probing extreme events in turbulent dynamical systems, Science advances 3 (9) (2017) e1701533.
  * [38] T. P. Sapsis, Statistics of extreme events in fluid flows and waves, Annual Review of Fluid Mechanics 53.
  * [39] R. Das, S. S. Girimaji, On the reynolds number dependence of velocity-gradient structure and dynamics, Journal of Fluid Mechanics 861 (2019) 163–179.
  * [40] R. Das, S. S. Girimaji, Revisiting turbulence small-scale behavior using velocity gradient triple decomposition, New Journal of Physics.