aff1]Physikalisches Institut, Albert-Ludwigs-Universität Freiburg, Hermann-
Herder-Str. 3, 79104 Freiburg, Germany aff2]EUCOR Centre for Quantum Science
and Quantum Computing, Hermann-Herder-Str. 3, 79104 Freiburg, Germany
aff3]Institut für Theoretische Physik III, Heinrich-Heine-Universität
Düsseldorf, Universitätsstr. 1, 40225 Düsseldorf, Germany
aff4]Naturwissenschaftlich-Technische Fakultät, Universität Siegen, Walter-
Flex-Str. 3, 57068 Siegen, Germany

# Entanglement characterization using quantum designs

Andreas Ketterer  Nikolai Wyderka  Otfried Gühne  [  [  [  [

###### Abstract

We present in detail a statistical approach for the reference-frame-
independent detection and characterization of multipartite entanglement based
on moments of randomly measured correlation functions. We start by discussing
how the corresponding moments can be evaluated with designs, linking methods
from group and entanglement theory. Then, we illustrate the strengths of the
presented framework with a focus on the multipartite scenario. We discuss a
condition for characterizing genuine multipartite entanglement for three
qubits, and we prove criteria that allow for a discrimination of $W$-type
entanglement for an arbitrary number of qubits.

##  1 Introduction

The experimental detection of multipartite entanglement usually requires a
number of appropriately chosen local quantum measurements which are aligned
with respect to a previously shared common reference frame [1, 2]. The latter,
however, can be a challenging prerequisite for photonic free-space quantum
communication over distances of several hundreds of kilometers [3, 4], which
is currently in the process of being extended to space involving satellites
orbiting the earth [5, 6, 7, 8]. Here, due to the motion, distance and number
of involved satellites, the issue of sharing classical reference frames
becomes particularly challenging, making the development of alternative
detection strategies desirable.

In recent years, there has been a number of proposals of experimental
protocols that avoid the need of sharing classical reference frames. One
possibility is to encode logical qubits into rotational invariant subspaces of
combined degrees of freedom of photons, i.e. their polarization and transverse
degrees of freedom [9, 10]. While the latter procedure provides one with a
complete experimental toolbox for alignment free quantum communication, one
can also find experimentally less demanding strategies that allow for the
reference frame independent certification of entanglement. For instance, one
can use entanglement criteria that are invariant under local unitary (LU)
transformations, commonly termed as reference-frame-independent [11, 12, 13,
14, 15, 16, 17, 18]. This type of entanglement criteria requires that the
experimenters are capable of measuring a fixed set of local observables, but
completely avoids the need of aligning measurements among different parties.

In several recent works it has been shown how to go beyond such procedures by
relaxing also the assumption of being able to measure a fixed set of local
observables and instead allow only for local measurements with settings drawn
uniformly at random [19, 20, 21, 22, 23, 24, 25, 27, 26]. The common idea of
these approaches is to measure a certain correlation function, and average the
result over random local unitaries applied to the state. Clearly, in this
scenario one has to resort to statistical tools based on the moments of the
resulting probability distribution of correlations in order to infer the
nonlocal properties of the underlying quantum states (see also Fig. 1).
Furthermore, apart from its reference-frame-independent nature, such protocols
are advantageous for the characterization of large multipartite systems where
a complete reconstruction of the underlying quantum state becomes practically
impossible due to the required measurement resources.

So far, however, most of the aforementioned approaches have been focused on
the lowest statistical moments only [19, 20, 23, 24, 25, 27], and applied them
to the detection of either the absence of full-separability [19, 20], or bi-
separability with respect to a predetermined bipartition [23, 24, 25, 27].
Though the latter approaches have proven useful for experimental
implementations in trapped ion experiments [24, 27], it is important to note
that the knowledge of all moments corresponding to different sectors, i.e.,
subsets of the involved parties, gives more insight into the entanglement
properties of the underlying state. For the case of second moments, this
follows from the relation between the respective moments and so-called sector
lengths which have been independently under investigation in the context of
entanglement detection [28, 29, 30]. Similar insights have been reached in
Ref. [26] where the detection and characterization of multipartite
entanglement based on second moments was studied experimentally with entangled
photons.

Figure 1: Example for entanglement detection using the statistics of random
measurements. For a given two-qubit state, one can measure the expectation
value
$E({\boldsymbol{u}}_{1},{\boldsymbol{u}}_{2})=\langle\sigma_{{\boldsymbol{u}}_{1}}\otimes\sigma_{\boldsymbol{u}_{2}}\rangle$
in randomly chosen local directions ${\boldsymbol{u}}_{1}$ and
${\boldsymbol{u}_{2}}$. The figure shows a Histogram of the observed counts
for $10^{4}$ measured values of $E({\boldsymbol{u}}_{1},{\boldsymbol{u}}_{2})$
with uniformly sampled directions ${\boldsymbol{u}}_{1}$ and
${\boldsymbol{u}}_{2}$. Colors indicate three different two qubit states:
product state (blue), the two-body marginal state of the tripartite $W$-state
(orange) and a specific Werner state (green). The distributions share the same
mean value and, in fact, also the same variance (as one has
$\mathcal{R}^{(2)}=1/9$ in Eq. (4)), but the forth moment $\mathcal{R}^{(4)}$
differs and allows to detect the entanglement of the Werner state and the
marginal state.

Moreover, it has recently been demonstrated that a combination of statistical
moments beyond the second order can lead to further improvements in terms of
entanglement detection [21]. In particular, it was shown that novel reference-
frame-independent criteria for the detection and also characterization of
multipartite entanglement can be derived by expressing the respective moments
in terms of spherical designs, i.e., pseudo-random processes allowing to mimic
uniform averages over the sphere. This insight also has experimental
consequences because the evaluation of the respective moments in terms of
spherical designs is exact and, for multipartite systems involving a small or
intermediate sized number of parties, does not exploit too many measurement
resources. Lastly, it was also shown that a quantification of bipartite
entanglement in terms of the respective moments is possible as well [22].

In the present contribution we will further investigate the potential of this
framework. After recalling for completeness some of the results from Ref.
[21], we first discuss the characterization of multipartite entanglement
classes based on the first two non-vanishing moments in small multipartite
systems consisting of three and four qubits. Finally, we prove two novel
criteria enabling the discrimination of $W$-type entangled mixed states for an
arbitrary number of parties.

The paper is organized as follows. In Sec. 2 we introduce the necessary
theoretical tools from Ref. [21], i.e., the moments of random correlation
functions and the concept of unitary and spherical designs (see Sec. 2.1 and
2.2, respectively). These are then used in Sec. 3 to evaluate moments of
random correlation functions of different order. Section 4 then focuses on the
application of the introduced framework to few-qubit system discussing the
characterization of multipartite entanglement based on the first two non-
vanishing moments. Lastly, the criterion enabling the discrimination of
$W$-type entanglement in multi-qubit systems is introduced in Sec. 5. Finally,
we conclude in Sec. 6 and give a short outlook.

##  2 Theoretical framework

###  2.1 Moments of random correlations

To set the stage we consider a system of $N$ $d$-dimensional quantum systems
(qudits) prepared in the initial state $\rho$. Subsequently, each of the
qudits is measured in a randomly drawn basis

$\displaystyle\left\\{\left(\rvert u_{n}^{(0)}\rangle:=U_{n}\rvert
0_{n}\rangle,\rvert u_{n}^{(1)}\rangle:=U_{n}\rvert 1_{n}\rangle,\ldots,\rvert
u_{n}^{(d-1)}\rangle:=U_{n}\rvert(d-1)_{n}\rangle\right)\right\\}_{n=1,\ldots,N},$
(1)

each specified by a random unitary transformation $U_{n}$ picked uniformly
from the unitary group $\mathcal{U}(d)$, i.e., according to the Haar measure.
One round of such random measurements yields the corresponding correlation
function $\langle U_{1}\mathcal{O}U_{1}^{\dagger}\otimes\ldots\otimes
U_{N}\mathcal{O}U_{N}^{\dagger}\rangle$, where $\mathcal{O}$ describes an
arbitrary qudit observable diagonal in the computational basis $\\{\rvert
0_{n}\rangle,\ldots,\rvert(d-1)_{n}\rangle\\}$, and $\langle...\rangle$
denotes the expectation value with respect to the quantum state $\rho$. We
note that in general the choice of $\mathcal{O}$ is relevant if one considers
systems of local dimensions $d$. However, it turns out that in the case $d=2$,
which is the main focus of Secs. 4 and 5, any local qubit observable will
suffice111This is a direct consequence of the isomorphism between
SU$(2)/\mathbb{Z}_{2}$ and SO$(3)$..

Further on, as the unitary transformations $U_{n}$ are chosen randomly, a
single set of random measurement settings will not give much insight into the
nonlocal properties of the initial state $\rho$. In order to achieve this we
have to perform several rounds of random measurements and seek a statistical
treatment in terms of the moments of the randomly measured correlation
functions. As we have assumed that each correlation function is characterized
by a set of Haar random unitaries $\\{U_{n}\\}_{n=1,\ldots,N}$, we can define
these moments as follows [21]:

$\displaystyle\mathcal{R}^{(t)}$
$\displaystyle=\int_{\mathcal{U}(d)}d\eta(U_{1})\ldots\int_{\mathcal{U}(d)}d\eta(U_{N})\langle
U_{1}\mathcal{O}U_{1}^{\dagger}\otimes\ldots\otimes
U_{N}\mathcal{O}U_{N}^{\dagger}\rangle^{t},$ (2)

where $t$ is a positive integer, and $\eta$ the Haar measure on the unitary
group $\mathcal{U}(d)$. It is important to note that one can gain more
information about $\rho$ by considering also moments of smaller qubit sectors,
i.e., the respective reduced states, as has been investigated for $t=2$ in
Refs. [28, 24, 25, 27, 29, 30, 26]. However, in the remainder of this
manuscript we will focus on the characterization of multipartite entanglement
based on full $N$-qubit moments (2).

Further on, we note that the random observables $U\mathcal{O}U^{\dagger}$,
with $U\in\mathcal{U}(d)$, can be parametrized by $d(d-1)$ angles. This
follows directly from the fact that any $U\in\mathcal{U}(d)$ can be decomposed
as $U=e^{i\varphi}Z_{1}XZ_{2}$, where $\varphi$ is a global phase, $Z_{1}$ and
$Z_{2}$ are diagonal unitary matrices with $[Z_{1}]_{11}=[Z_{2}]_{11}=1$ and
$X$ is a unitary matrix with $\sum_{i}[X]_{ij}=\sum_{i}[X]_{ji}=1$, for all
$j$ [31]. We thus have that
$U\mathcal{O}U^{\dagger}=Z_{1}X\mathcal{O}X^{\dagger}Z_{1}^{\dagger}$ and,
since the matrices $A$ are isomorphic to $\mathcal{U}(d-1)$, a simple count of
parameters leads to $d(d-1)$. Hence, in the case of two-level systems (qubits)
we find that the local measurement settings are characterized by only two
angles corresponding to the spherical coordinates fixing a direction on the
Bloch sphere $S^{2}$. In other words, we can associate to each random basis
$(\rvert u_{n}^{(0)}\rangle:=U_{n}\rvert 0_{n}\rangle,\rvert
u_{n}^{(1)}\rangle:=U_{n}\rvert 1_{n}\rangle)$, a direction
$\boldsymbol{u}_{n}\in S^{2}$, defined by the components
$[\boldsymbol{u}_{n}]_{i}=\mathrm{tr}[\sigma_{\boldsymbol{u}_{n}}\sigma_{i}]/2$,
where $\sigma_{i}$, with $i=x,y,z$, denote the usual Pauli matrices and
$\sigma_{\boldsymbol{u}_{n}}=U_{n}\sigma_{z}U_{n}^{\dagger}$. Equation (2) for
the $t$-th moment thus becomes

$\displaystyle\mathcal{R}^{(t)}$
$\displaystyle=\int_{\mathcal{U}(2)}d\eta(U_{1})\ldots\int_{\mathcal{U}(2)}d\eta(U_{N})\langle
U_{1}\sigma_{z}U_{1}^{\dagger}\otimes\ldots\otimes
U_{N}\sigma_{z}U_{N}^{\dagger}\rangle^{t}$ (3)
$\displaystyle=\frac{1}{(4\pi)^{N}}\int_{S^{2}}d\boldsymbol{u}_{1}\ldots\int_{S^{2}}d\boldsymbol{u}_{N}E(\boldsymbol{u}_{1},\ldots,\boldsymbol{u}_{N})^{t},$
(4)

where we defined
$E(\boldsymbol{u}_{1},\ldots,\boldsymbol{u}_{N}):=\langle\sigma_{\boldsymbol{u}_{1}}\otimes\ldots\otimes\sigma_{\boldsymbol{u}_{N}}\rangle$,
and $d\boldsymbol{u}_{i}=\sin{\theta_{i}}d\theta_{i}d\phi_{i}$ denotes the
uniform measure on the Bloch sphere $S^{2}$. In particular, it is easy to see
from Eq. (4) that all odd moments are zero due to the symmetry of the
correlation functions $E(\boldsymbol{u}_{1},\ldots,\boldsymbol{u}_{N})$ with
respect to a reflection on the Bloch sphere:
$E(\boldsymbol{u}_{1},\ldots,-\boldsymbol{u}_{i},\ldots,\boldsymbol{u}_{N})=-E(\boldsymbol{u}_{1},\ldots,\boldsymbol{u}_{N})$.

In the following, we show that the moments $\mathcal{R}^{(t)}$ can be
calculated using unitary $t$-designs (or spherical $t$-designs for $d=2$),
rather than averaging over the whole unitary group.

###  2.2 Designs

Generally speaking quantum designs are pseudo-random processes that allow to
mimic uniform averages over some group if one is only interested in moments up
to some finite degree. Depending on the specific choice of the group such
processes are either referred to as unitary or spherical designs. In the
following, we will give short introduction to the two of them.

####  2.2.1 Unitary designs

Let us denote by $\mathrm{Hom}(r,s)$ the set of all homogeneous polynomials
$P_{r,s}(U)$, with support on the space of unitary matrices $\mathcal{U}(d)$,
that is of degree at most $r$ and $s$, respectively, in each of the matrix
elements of $U$ and their complex conjugates. For example, the polynomial
$P_{r=1,s=2}(U)=U^{\dagger}VUVU^{\dagger}$ is of degree $r=1$ and $s=2$ in the
entries of the matrix $U$. With this we arrive at the following definition
[32]: A unitary $t$-design is a set of unitary matrices
$\\{U_{k}|k=1,\ldots,K^{(t)}\\}\subset\mathcal{U}(d)$, with cardinality
$K^{(t)}$, such that

$\displaystyle\frac{1}{K^{(t)}}\sum_{k=1}^{K^{(t)}}P_{t^{\prime},t^{\prime}}(U_{k})=\int_{\mathcal{U}(d)}P_{t^{\prime},t^{\prime}}(U)d\eta(U),$
(5)

for all homogeneous polynomials
$P_{t^{\prime},t^{\prime}}\in\mathrm{Hom}(t^{\prime},t^{\prime})$, with
$t^{\prime}\leq t$, and where $\eta(U)$ denotes the normalized Haar measure on
$\mathcal{U}(d)$.

We note that, while the existence of unitary designs has been proven [33], no
universal strategy for their construction in case of an arbitrarily given $t$
exists. This fact led to the study of approximate unitary designs for which
the property (5) is accordingly relaxed [34, 35, 36]. However, in the
remainder of this manuscript we will restrict ourselves to particular cases in
which exact design are known. A prominent example of an exact unitary
$t$-design is given by the multi-qubit Clifford group consisting of all
unitary matrices mapping the multi-qubit Pauli group onto itself. The latter
has been shown to constitute a unitary $3$-design [37], however, it fails to
be a unitary $4$-design [38]. In the case of a single qubit the Clifford group
has $24$ elements which can be generated from the Hadamard gate $H$ and the
phase gate $S=e^{i\frac{\pi}{4}\sigma_{z}}$. Furthermore, in Ref. [39] the
existence of a qubit $5$-design of one qubit was noted. The latter is given by
the unitary representation of the special linear group $SL(2,\mathbb{F}_{5})$
of invertible $2\times 2$ matrices over the finite field $\mathbb{F}_{5}$ with
five elements. In App. A we shortly outline how to generate this design.

####  2.2.2 Spherical designs

Figure 2: Plot of the Bloch vectors of various spherical designs, for $d=2$,
and their corresponding polyhedra. The components of the Bloch vector are
denoted as $\langle\sigma_{i}\rangle$, with the Pauli matrices $\sigma_{x}$,
$\sigma_{y}$, and $\sigma_{z}$. (a) $3$-design with $K_{\text{octa}}=6$
vertices forming an octahedron. (b) $5$-design with $K_{\text{icosa}}=12$
vertices forming an icosahedron. (c) $5$-design with $K_{\text{icosi}}=30$
vertices forming an icosidodecahedron. (d) $7$-design with
$K_{7-\text{design}}=24$ vertices forming a deformed snub cube (red points).
For comparison, a regular snub cube only forms a $3$-design (blue polyhedron).

In Sec. 2.1 we saw that in the case of systems consisting of qubits the
evaluation of the moments $\mathcal{R}^{(t)}$ boils down to local integrations
over the spheres $S^{2}$. In this case, instead of using a unitary $t$-design
to evaluate the respective moments, we can resort to the concept of spherical
designs. In general, a spherical $t$-design in dimension three consist of a
finite set of points $\\{\boldsymbol{u}_{i}|i=1,\ldots,L^{(t)}\\}\subset
S^{2}$ fulfilling the property

$\displaystyle\frac{1}{L^{(t)}}\sum_{k=1}^{L^{(t)}}P_{t^{\prime}}(\boldsymbol{u}_{k})=\frac{1}{4\pi}\int_{S^{2}}d\boldsymbol{u}\
P_{t^{\prime}}(\boldsymbol{u}),$ (6)

for all homogeneous polynomials $P_{t^{\prime}}:S^{2}\rightarrow\mathbb{R}$,
with $t^{\prime}\leq t$. It thus suffices to resort to spherical $t$-designs
as long as one is interested in calculating averages of polynomials of degree
at most $t$ over the Bloch sphere $S^{2}$.

One way to generate spherical designs is to extract them from unitary designs.
For instance, by applying the elements of the single-qubit Clifford group to
one of the Pauli matrices, e.g. $\sigma_{z}$, we are left with the following
set of inequivalent operators
$\\{\pm\sigma_{x},\pm\sigma_{y},\pm\sigma_{z}\\}$. The latter correspond to
the following set of unit vectors $\\{\pm\boldsymbol{e}_{i}|i=x,y,z\\}$ which
form a spherical $3$-design (see Fig. 2(a)). Similarly, one can generate a
spherical $5$-design from the $60$ element unitary $5$-design
$SL(2,\mathbb{F}_{5})$. To do so, we calculate again all inequivalent
directions on the Bloch sphere originating from the operators
$\sigma_{{\boldsymbol{u}}^{(k)}}=U^{(k)}\sigma_{z}{U^{(k)}}^{\dagger}$, for
all elements $U^{(k)}$ of the unitary $5$-design. The result is a set of 30
vertices on the Bloch sphere forming an icosidodecahedron (see Fig. 2(c)).

In general, however, spherical designs are easier to find than unitary designs
because it is easier to search them directly by checking the relation (6) for
sets of vertices on the sphere $S^{2}$. Such a search was carried out in Ref.
[40] where a number of spherical designs on the $2$-sphere $S^{2}$ for $t$’s
up to 20 and consisting of up to 100 elements were found. In Fig. 2(b) and (d)
we present two more examples of such spherical designs. The $5$-designs
presented in Fig. 2(b) corresponds to an icosahedron with
$K_{\text{icosa}}=12$ vertices. Furthermore, the $7$-design with
$K_{7-\text{design}}=24$ vertices, shown in Fig. 2(d), represents a deformed
snub cube, i.e., a regular snub cube with slightly smaller square faces and
slightly larger triangular faces [40]. Interestingly, the icosidodecahedron,
presented in Fig. 2(c), and the regular snub cube, presented in blue in Fig.
2(d), constitute only spherical $5$\- and $3$-designs, respectively.

Lastly, we emphasize that for systems of larger local dimensions it is
difficult to resort to spherical designs instead of unitary ones for the
evaluation of the moments $\mathcal{R}^{(t)}$. This is due to the fact that
the parametrization of the space of observables of the form
$U\mathcal{O}U^{\dagger}$ (see Sec. 2.1), is not in one to one correspondence
with points on a generalized Bloch sphere $S^{d^{2}-1}$. It rather forms a
submanifold of $S^{d^{2}-1}$ characterized by $d(d-1)$ parameters. Hence,
generalized spherical designs on $S^{d^{2}-1}$ will not be useful in this
case.

##  3 Random moments from designs

We first note that the $t$-th power of the correlation function $\langle
U_{1}\mathcal{O}U_{1}^{\dagger}\otimes\ldots\otimes
U_{N}\mathcal{O}U_{N}^{\dagger}\rangle$ is a polynomial of degree $t$ in the
entries of the local random unitary matrices $U_{n}\in\mathcal{U}(d)$, and
their complex conjugates. Hence, in order to evaluate the moments
$\mathcal{R}^{(t^{\prime})}$, with $t^{\prime}\leq t$, it suffices to average
$\langle U_{1}\mathcal{O}U_{1}^{\dagger}\otimes\ldots\otimes
U_{N}\mathcal{O}U_{N}^{\dagger}\rangle^{t^{\prime}}$ locally over a respective
unitary $t$-design instead over the whole unitary group $\mathcal{U}(d)$. This
leads to the result

$\displaystyle\mathcal{R}^{(t^{\prime})}$
$\displaystyle=\frac{1}{(K^{(t)})^{N}}\sum_{k_{1},\ldots,k_{N}=1}^{K^{(t)}}\langle
U_{k_{1}}\mathcal{O}U_{k_{1}}^{\dagger}\otimes\ldots\otimes
U_{k_{N}}\mathcal{O}U_{k_{N}}^{\dagger}\rangle^{t^{\prime}},$ (7)

for $t^{\prime}\leq t$, and with a unitary $t$-design
$\\{U_{n}\\}_{k=1}^{K^{(t)}}$. Equation (7) is a general formula that allows
one to calculate the moments (4) up to a certain $t$ provided a unitary
$t$-design can be found. In the remainder of this paper we will focus on the
qubit case ($d=2$) for which Eq. (7) takes a particularly simple form. In this
case Eq. (7) becomes

$\displaystyle\mathcal{R}^{(t)}$
$\displaystyle=\frac{1}{(L^{(t)})^{N}}\sum_{k_{1},\ldots,k_{N}=1}^{L^{(t)}}\langle\sigma_{\boldsymbol{u}_{k_{1}}}\otimes\ldots\otimes\sigma_{\boldsymbol{u}_{k_{N}}}\rangle^{t},$
(8)

where $L^{(t)}\leq K^{(t)}$ denotes the count of the remaining non-equivalent
measurement directions $\\{\boldsymbol{u}_{i}\\}_{i=1,\ldots,L^{(t)}}$, after
dropping those for which $j$’s exist with
$\sigma_{\boldsymbol{u}_{k_{i}}}=\sigma_{\boldsymbol{u}_{k_{j}}}$. Hence, in
order to calculate the moments $\mathcal{R}^{(t)}$ for systems of qubits it is
enough to average
$\left[E(\boldsymbol{u}_{1},\ldots,\boldsymbol{u}_{N})\right]^{t}$ over a
finite number $L^{(t)}$ of nonequivalent Bloch sphere directions
$\\{\boldsymbol{u}_{i}|i=1,\ldots,L^{(t)}\\}\subset S^{2}$ which themselves
form a spherical $t$-design.

We explicitly evaluate the moments (8) using the $6$ element spherical
$3$-design and $12$ element spherical $5$-design, presented in Sec. 2.2,
yielding:

$\displaystyle\mathcal{R}^{(2)}$
$\displaystyle=\frac{1}{3^{N}}\sum_{i_{1},\ldots,i_{N}=x,y,z}E(\boldsymbol{e}_{i_{1}},\ldots\boldsymbol{e}_{i_{N}})^{2},$
(9) $\displaystyle\mathcal{R}^{(4)}$
$\displaystyle=\frac{1}{6^{N}}\sum_{i_{1},\ldots,i_{N}=1}^{6}E(\boldsymbol{v}_{i_{1}},\ldots\boldsymbol{v}_{i_{N}})^{4},$
(10)

with the upper summation bounds given by $L^{(3)}/2=3$ and $L^{(5)}/2=6$,
respectively. Note that in Eqs. (9) and (10) the number of summands is
$L^{(t)}/2$ because for even $t$ one can drop the respective anti-parallel
settings $-\boldsymbol{e}_{i}$ and $-\boldsymbol{v}_{i}$. Also note that we
always use those spherical designs with the smallest number of elements, in
order to minimize the number of involved measurement settings.

In the remainder of the manuscript we will mainly focus on entanglement
criteria involving the first two non-vanishing moments, i.e.,
$\mathcal{R}^{(2)}$ and $\mathcal{R}^{(4)}$. However, we emphasize that is
generally possible to evaluate also higher order moments using higher order
designs. For instance, with the deformed snub cube spherical $7$-design
$\\{\boldsymbol{w}_{i}|i=1,\ldots,L^{(7)}=24\\}$ (see Fig. 2(d)) we obtain the
following formula for the sixth moment

$\displaystyle\mathcal{R}^{(6)}=\frac{1}{24^{N}}\sum_{i_{1},\ldots,i_{N}=1}^{24}E(\boldsymbol{w}_{i_{1}},\ldots\boldsymbol{w}_{i_{N}})^{6}.$
(11)

Note that the number of summands in Eq. (11) cannot be reduced to $L^{(7)}/2$,
as in Eqs. (9) and (10), because the deformed snub cube is not point
symmetric.

##  4 Applications to qubit systems

###  4.1 Two qubits

A general two-qubit density matrix can be expressed in the following form

$\displaystyle\rho=\frac{1}{4}\big{[}\mathbbm{1}_{4}+(\boldsymbol{a}\cdot\boldsymbol{\sigma})\otimes\mathbbm{1}_{2}+\mathbbm{1}_{2}\otimes(\boldsymbol{b}\cdot\boldsymbol{\sigma})+\sum_{i,j=x,y,z}c_{i,j}\sigma_{i}\otimes\sigma_{j}\big{]},$
(12)

where $\boldsymbol{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})^{\top}$ denotes
the Pauli spin operator, $\boldsymbol{a}$ and $\boldsymbol{b}$ the local Bloch
vectors of the reduced subsystems of each individual qubit, and
$C=\\{c_{i,j}\\}_{i,j=x,y,z}=\\{\text{tr}[\rho\sigma_{i}\otimes\sigma_{j}]\\}_{i,j=x,y,z}$
the correlation matrix. To start, we focus first on so-called Bell diagonal
states, i.e. states that are diagonal in the Bell basis and whose density
matrix can be expressed in the following way
$\rho_{\text{BD}}=\frac{1}{4}\big{[}\mathbbm{1}_{4}+\sum_{j=x,y,z}c_{j}\sigma_{j}\otimes\sigma_{j}\big{]}$,
with real parameters $c_{j}$, such that $0\leq|c_{j}|\leq 1$, and the
corresponding eigenvalues $\lambda_{1,2}=(1\mp c_{x}\mp c_{y}-c_{z})/4$ and
$\lambda_{3,4}=(1\pm c_{x}\mp c_{y}+c_{z})/4$. Bell diagonal states are
separable if and only if $|c_{x}|+|c_{y}|+|c_{z}|\leq 1$ [41].

Figure 3: (a) Representation of the set of Bell diagonal states in the space
spanned by the moments $\mathcal{R}^{(2)}$ and $\mathcal{R}^{(4)}$ obtained
from the analytical solution of the system of Eqs. (15)-(17). Labeled black
circles indicate the maximally mixed state (A), the pure product states (B)
and the Bell states (C). (b) Zoom into plot (a) in the range
$0\leq\mathcal{R}^{(2)}\leq 1/3^{2}$ highlighting separable (blue solid lines)
and entangled (red dotted lines) Bell diagonal states. Note that for
$\mathcal{R}^{(2)}\leq 1/3^{3}$ all states are separable, and for
$1/3^{3}\leq\mathcal{R}^{(2)}\leq 1/3^{2}$ separable and entangled states have
a non-zero overlap (striped region). The white dashed curve indicates the
entanglement criteria resulting form Eqs. (15)-(17).

Direct evaluation of Eq. (9) and (10) for Bell diagonal states yields:

$\displaystyle\mathcal{R}^{(2)}$
$\displaystyle=\frac{1}{9}(c_{x}^{2}+c_{y}^{2}+c_{z}^{2}),$ (13)
$\displaystyle\mathcal{R}^{(4)}$
$\displaystyle=\frac{2}{75}(c_{x}^{4}+c_{y}^{4}+c_{z}^{4})+\frac{27}{25}(\mathcal{R}^{(2)})^{2},$
(14)

which are polynomials of degree smaller or equal than four in the coefficients
$c_{i}$, with $i=x,y,z$. In order to determine the borders of the set of Bell
diagonal states in the space spanned by the moments $\mathcal{R}^{(2)}$ and
$\mathcal{R}^{(4)}$ we have to maximize (minimize), e.g., $\mathcal{R}^{(4)}$
over the set of states while keeping $\mathcal{R}^{(2)}$ fixed. The latter
optimization has been carried out in Refs. [21, 22] and is presented in Fig.
3(a). Similarly, we can obtain the borders of the set of separable states by
additionally imposing the condition $|c_{x}|+|c_{y}|+|c_{z}|\leq 1$. It turns
out that solving this optimization problem is equivalent to finding the
solution of the system of equations:

$\displaystyle 1$ $\displaystyle=|c_{x}|+|c_{y}|+|c_{z}|,$ (15)
$\displaystyle\alpha$ $\displaystyle=c_{x}^{2}+c_{y}^{2}+c_{z}^{2},$ (16)
$\displaystyle\beta$ $\displaystyle=c_{x}^{4}+c_{y}^{4}+c_{z}^{4},$ (17)

where we used the transformations $\alpha:=3^{2}\mathcal{R}^{(2)}$ and
$\beta:=15^{2}\mathcal{R}^{(4)}/6-\alpha^{2}/2$. The relevant solutions of
Eqs. (15)-(17) have been derived in Ref. [21] and are presented in detail in
Fig. 3(b). We thus obtain a sufficient entanglement criterion which is a
polynomial function of the moments $\mathcal{R}^{(2)}$ and
$\mathcal{R}^{(4)}$. We note that the latter criterion can be made necessary
and sufficient by including the next higher non-vanishing moment
$\mathcal{R}^{(6)}$ [21, 22].

Further on, in the case of general two-qubit states (12) one can always find
corresponding Bell diagonal states $\rho_{\text{BD}}$ having the same moments.
On the one hand, this is a consequence of the fact that the moments
$\mathcal{R}^{(t)}$ are LU invariant which allows us to diagonalize the
correlation matrix $C=\\{c_{i,j}\\}_{i,j=x,y,z}$. On the other hand, we can
subsequently drop those terms in Eq. (12) involving only the local Bloch
vector components $\boldsymbol{a}$ and $\boldsymbol{b}$ by applying a local
dephasing operation
$\rho\rightarrow\frac{1}{4}(\rho+\sum_{i=x,y,z}\sigma_{i}\otimes\sigma_{i}\rho\sigma_{i}\otimes\sigma_{i})$.
As the latter transformations are completely positive and, in particular,
correspond to the class of LOCC operations, we end up with a resulting Bell
diagonal state that has equal moments as the original two-qubit state but
whose entanglement cannot have increased. Hence, the derived entanglement
criteria for Bell diagonal states yield also sufficient entanglement criteria
for general two-qubit states.

Lastly, we note that the results presented in Fig. 3 can be generalized in
order to account for a quantification of entanglement. In Ref. [22] it was
shown that there exists a one-to-one correspondence between the concurrence of
a Bell diagonal state and its corresponding moments $\mathcal{R}^{(2)}$,
$\mathcal{R}^{(4)}$ and $\mathcal{R}^{(6)}$. Likewise, for general two-qubit
states $\rho$ this one-to-one correspondence turns into a sufficient criterion
allowing to lower bound the concurrence of $\rho$.

###  4.2 Three qubits

We move on now to the more involved multipartite scenarios and start with the
simplest case consisting of three qubits. However, already for three qubits a
complete characterization of the set of states in the space spanned by the
moments $\mathcal{R}^{(2)}$ and $\mathcal{R}^{(4)}$ becomes considerably more
difficult due to the increased Hilbert space dimension and the more involved
structure of the space of multipartite entangled states. Nevertheless, there
have been some formulations of criteria allowing to detect three-qubit
entanglement and also genuine three-qubit entanglement based on the second
moment only [14, 18, 29]. In the following, we will go beyond these results
and discuss the characterization of different classes of genuine multipartite
entanglement in terms of $\mathcal{R}^{(2)}$ and $\mathcal{R}^{(4)}$.

Classes of genuine multipartite entangled states are usually defined through
the concept of stochastic local operations and classical communication (SLOCC)
[44, 45] which form the probabilistic counterpart of LOCC operations [42, 43].
In this framework, two pure $N$-qubit states, $\rvert\Psi\rangle$ and
$\rvert\Phi\rangle$, are equivalent if there exist LOCC operations allowing to
transform them into each other with some finite probability. Mathematically,
SLOCC equivalence implies the existence of invertible operations $A_{i}$, with
$i=1,\ldots,N$, such that [44]

$\displaystyle\rvert\Psi\rangle=A_{1}\otimes\ldots\otimes
A_{N}\rvert\Phi\rangle.$ (18)

The states are then called SLOCC equivalent and the corresponding equivalence
classes are referred to as SLOCC classes. For systems consisting of three
qubits there exist in total six SLOCC classes: the pure product states
$\mathcal{S}^{(3)}$, three classes of bi-separable states
$\mathcal{B}_{1|23}$, $\mathcal{B}_{12|3}$ and $\mathcal{B}_{13|2}$, and two
classes of genuinely multipartite entangled states, the $W$\- and the GHZ-
class [44], referred to as $\mathcal{W}^{(3)}$ and $\mathcal{GHZ}^{(3)}$,
respectively. As their names suggest, $\mathcal{W}^{(3)}$ and
$\mathcal{GHZ}^{(3)}$ consist of those states which are SLOCC equivalent to
the states $\rvert W_{3}\rangle=(\rvert 001\rangle+\rvert 010\rangle+\rvert
100\rangle)/\sqrt{3}$ and $\rvert\text{GHZ}_{3}\rangle=(\rvert
000\rangle+\rvert 111\rangle)/\sqrt{2}$, respectively.

Figure 4: The set of (bi-separable) three-qubit states on the plane spanned by
$\mathcal{R}^{(2)}$ and $\mathcal{R}^{(4)}$. The subset of bi-separable states
is displayed in light blue with thick blue boundaries. The lower bound of the
set of all states is displayed in solid black, and $150,000$ random three-
qubit states are displayed in gray. The conjectured criterion in Eq. (19) is
displayed as a dashed, green line. The red circles label the maximally mixed
state (A), all pure product states (B), bi-separable states of the form
$\rvert\phi\rangle\rvert\text{Bell}\rangle$ (C), the uniform mixture of the
bi-separable states $\frac{1}{\sqrt{2}}\rvert 0\rangle(\rvert 00\rangle+\rvert
11\rangle)$ and $\frac{1}{\sqrt{2}}\rvert 1\rangle(\rvert 01\rangle+\rvert
10\rangle)$ (D), the three-qubit $W$-state (E) and the GHZ-state (F).

For mixed states one can define similar classes by taking the convex hull of
the corresponding SLOCC classes [45]. In this way, we obtain the usual set of
separable mixed states $\text{Conv}(\mathcal{S}^{(3)})$, the set of all bi-
separable mixed states $\text{Conv}(\mathcal{B}_{\text{bi-sep}}^{(3)})$, with
$\mathcal{B}_{\text{bi-sep}}^{(3)}:=\mathcal{B}_{1|23}\cup
B_{12|3}\cup\mathcal{B}_{13|2}$, and consequently
$\text{Conv}(\mathcal{W}^{(3)})$ and $\text{Conv}(\mathcal{GHZ}^{(3)})$, where
the convex hull is defined as $\text{Conv}(X):=\\{\sum_{i}p_{i}x_{i}|x_{i}\in
X,p_{i}\geq 0,\sum_{i}p_{i}=1\\}$. This definition leads to the onion like
structure of the set of three-qubit mixed states as presented in Ref. [45]. We
note that in other contexts one might define mixed state analogs of
$\mathcal{W}^{(3)}$ and $\mathcal{GHZ}^{(3)}$ as
$\text{Conv}(\mathcal{W}^{(3)})\setminus\text{Conv}(\mathcal{B}_{\text{bi-
sep}}^{(3)})$ and
$\text{Conv}(\mathcal{GHZ}^{(3)})\setminus\text{Conv}(\mathcal{B}_{\text{ bi-
sep}}^{(3)})$, respectively, in order to ensure that they only contain states
which are genuinely multipartite entangled.

Figure 5: Representations of different three-qubit SLOCC classes in the space
spanned by the moments $\mathcal{R}^{(2)}$ and $\mathcal{R}^{(4)}$. Plotted
points correspond to the moments of randomly generated quantum states from the
class of separable states (a), the pure $W$-class $\mathcal{W}^{(3)}$ (b), its
convex hull $\text{Conv}(\mathcal{W}^{(3)})$ (c), and the whole state space of
three-qubit density matrices (d). Labeled red circles indicate the maximally
mixed state (A), all pure product states (B), bi-separable states of the form
$\rvert\phi\rangle\rvert\text{Bell}\rangle$ (C), the three-qubit $W$\- (D) and
GHZ-state (E) (see App. B for details on the sampling of the random states).
While the black solid lines connecting (B), (C) and (D) enclose the set
$\mathcal{W}^{(3)}$, its mixed extension $\text{Conv}(\mathcal{W}^{(3)})$ is
indicated by black dashed lines.

As genuine multipartite entanglement constitutes a necessary ingredient for
many applications [46, 47], we start by distinguishing the set of bi-separable
states $\text{Conv}(\mathcal{B}_{\text{bi-sep}}^{(3)})$ from the set of
genuinely multipartite entangled states in
$\text{Conv}(\mathcal{GHZ}^{(3)})\setminus\text{Conv}(\mathcal{B}_{\text{bi-
sep}}^{(3)})$ and
$\text{Conv}(\mathcal{W}^{(3)})\setminus\text{Conv}(\mathcal{B}_{\text{bi-
sep}}^{(3)})$ in order to find new entanglement criteria based on higher
moments. To that end, we numerically optimize the moments $\mathcal{R}^{(2)}$
and $\mathcal{R}^{(4)}$ for mixed bi-separable states and compare the result
to an optimization over the set of all states. The results are depicted in
Figure 4 and show a clear difference between the two sets, allowing for
improved entanglement detection compared to any criterion based on
$\mathcal{R}^{(2)}$ alone. The boundary between the two sets can be
approximated by a quadratic function, which is displayed in the figure as
well. This yields the following conjectured entanglement criterion for three-
qubit states:

$\mathcal{R}^{(4)}_{\rho_{\text{bisep}}}\geq\frac{1}{425}[972(\mathcal{R}^{(2)}_{\rho_{\text{bisep}}})^{2}+90\mathcal{R}^{(2)}-5].$
(19)

We note that it is possible to prove a similar bound on the set of pure bi-
separable states. To do so, one simply assumes an arbitrary product of single-
and two-qubit states
$\rvert\text{Bisep}\rangle=\rvert\phi\rangle\rvert\Psi\rangle$ for which the
respective fourth moment factorizes to
$\mathcal{R}^{(4)}_{\rvert\text{Bisep}\rangle}=\frac{1}{5}\mathcal{R}^{(4)}_{\rvert\Psi\rangle}$.
Subsequently, one can apply the results of Sec. 4.1 in order to obtain the
criterion $\mathcal{R}^{(4)}_{\rvert\text{Bisep}\rangle}\geq
5^{4}\mathcal{R}^{(2)}_{\rvert\text{Bisep}\rangle}/3^{5}$. The latter leads to
a slightly tighter curve than the one corresponding to Eq. (19), shown in Fig.
4, thus confirming the conjecture for pure separable states.

Next, we focus on the remaining SLOCC classes of fully separable states and
the two genuinely multipartite entangled ones. Figure 5 contains the results
of a numerical analysis aiming at identifying these classes through a mixture
of numerical optimizations and the generation of random states. For instance,
Fig. 5(b) shows a numerical estimation of the borders of the pure $W$-class
$\mathcal{W}^{(3)}$, which is confirmed by more than $10^{5}$ randomly
generated three-qubit states. Both, the numerical estimation of the borders of
the pure $W$-class and the generation of random states is based on the three-
qubit standard form presented in Ref. [48]. The latter allows us to
numerically determine the border of $\mathcal{W}^{(3)}$ by numerically
optimizing the moments with respect to the standard form parameters.
Furthermore, by drawing the same parameters randomly we were able to sample
states from $\mathcal{W}^{(3)}$ (see App. B for details). As the moments are
invariant under LU transformations, this suffices to sample exhaustively
points in the $(\mathcal{R}^{(2)},\mathcal{R}^{(4)})$-plane corresponding to
states in $\mathcal{W}^{(3)}$.

Further on, in Fig. 5(a), (c) and (d) we focus on the mixed SLOCC classes
$\text{Conv}(\mathcal{S}^{(3)})$, $\text{Conv}(\mathcal{W}^{(3)})$ and
$\text{Conv}(\mathcal{GHZ}^{(3)})$, respectively. However, due to the
structure of the convex hull these classes are more difficult to characterize.
We thus only estimated their boundaries roughly by minimizing the moments over
a subset of permutational invariant mixed states: $\text{Conv}(\\{\lvert
000\rangle\langle 000\rvert,\lvert W_{3}\rangle\langle
W_{3}\rvert,\mathbbm{1}_{8}/2^{3}\\})\subset\text{Conv}(\mathcal{W}^{(3)})$,
the results of which are indicated in Fig. 5 by black dashed lines. Moreover,
by sampling states from $\mathcal{S}^{(3)}$, $\mathcal{W}^{(3)}$ and
$\mathcal{GHZ}^{(3)}$, applying random LU transformations and subsequently
mixing them with randomly drawn mixing parameters, we are able to sample
states from the corresponding mixed classes. The latter seem to cover quiet
well the $(\mathcal{R}^{(2)},\mathcal{R}^{(4)})$-plane as presented in Fig. 5.
However, we stress that some areas of the
$(\mathcal{R}^{(2)},\mathcal{R}^{(4)})$-plane are underrepresented by this
procedure because of the low sampling probability of the corresponding states.
For instance, one can find sequences of states in
$\text{Conv}(\mathcal{W}^{(3)})$ which converge towards the uniform mixture of
the bi-separable states $\frac{1}{\sqrt{2}}\rvert 0\rangle(\rvert
00\rangle+\rvert 11\rangle)$ and $\frac{1}{\sqrt{2}}\rvert 1\rangle(\rvert
01\rangle+\rvert 10\rangle)$ (see also Fig. 4). Hence, it is expected that the
upper boundary of the set $\text{Conv}(\mathcal{W}^{(3)})$ shows a similar
structure as that of the mixed bi-separable sates presented in Fig. 4.

In conclusion, the above analysis shows that a discrimination of the mixed
SLOCC classes $\text{Conv}(\mathcal{S}^{(3)})$,
$\text{Conv}(\mathcal{W}^{(3)})$ and $\text{Conv}(\mathcal{GHZ}^{(3)})$ in
terms of $\mathcal{R}^{(2)}$ and $\mathcal{R}^{(4)}$ is possible.

###  4.3 Four qubits

In this section we carry out a similar analysis in the slightly more
complicate case consisting of four qubits. The main aim will be to give the
reader an impression of the growing complexity of the problem of
characterizing SLOCC classes in the space spanned by the moments
$\mathcal{R}^{(2)}$ and $\mathcal{R}^{(4)}$ by discussing important
differences to the results presented in Sec. 4.2.

Figure 6: Representations of different four-qubit SLOCC classes in the space
spanned by the moments $\mathcal{R}^{(2)}$ and $\mathcal{R}^{(4)}$. Plotted
points correspond to the moments of randomly generated quantum states from the
class of separable states (a), the pure $W$-class $\mathcal{W}^{(4)}$ (b), its
convex hull $\text{Conv}(\mathcal{W}^{(4)})$ (c), and the whole state space of
four-qubit density matrices (d). Labeled red circles indicate the maximally
mixed state (A), all pure product states (B), tri-separable states of the form
$\rvert\phi\rangle\rvert\phi^{\prime}\rangle\rvert\text{Bell}\rangle$ (C), bi-
separable states of the form $\rvert\phi\rangle\rvert W\rangle$ (D),
$\rvert\phi\rangle\rvert\text{GHZ}_{3}\rangle$ (E) and
$\rvert\text{Bell}\rangle\rvert\text{Bell}\rangle$ (F), and the four-qubit
$W$\- (G) and GHZ-state (H) (see App. B for details on the sampling of the
random states). Again, the black solid lines enclose the set
$\mathcal{W}^{(4)}$ and its mixed extension $\text{Conv}(\mathcal{W}^{(4)})$
is indicated by black dashed lines.

Figure 6 shows an overview of the characterization of four-qubit SLOCC
classes. From this we first note the remarkable resemblance between the set of
separable three- and four-qubit states (see Fig. 5(a) and 6(a)). In contrast,
the respective four-qubit $W$-classes, $\mathcal{W}^{(4)}$ and
$\text{Conv}(\mathcal{W}^{(4)})$, show some more detailed structures which
correspond to different types of bi- and tri-separable states (see Fig. 6(b)
and (c)). Furthermore, one recognizes that in the four-qubit case the relative
size of the pure, as well as the mixed $W$-classes is smaller in comparison to
the set of all states (see Fig. 6(d)). Lastly, we see that $\mathcal{R}^{(2)}$
is maximized simultaneously by the four qubit GHZ state
$\rvert\text{GHZ}_{4}\rangle$ and the bi-separable state
$\rvert\text{Bell}\rangle\rvert\text{Bell}\rangle$ consisting of a product of
two Bell states, the latter of which even maximizes $\mathcal{R}^{(4)}$. This
last point shows that it is no longer possible to detect genuine multipartite
entanglement based solely on $\mathcal{R}^{(2)}$, as its maximum within the
set $\text{Conv}(\mathcal{B}^{(4)}_{\text{bi-sep}})$ coincides with its
maximum over all four-qubit states. A similar conclusion can be reached for
$\mathcal{R}^{(4)}$, as its value for the state
$\rvert\text{Bell}\rangle\rvert\text{Bell}\rangle$ is larger than that for
$\rvert\text{GHZ}_{4}\rangle$. If the detection of genuine multipartite
entanglement is possible for more than three qubits by suitably combining
moments of different orders will be subject of future investigations.

In order to understand some of the above observations, it is important to note
that the structure of multipartite entanglement classes of four qubits is
already considerably more complicated than for three qubits. In fact, in the
case of four qubits we are dealing with infinitely many SLOCC equivalence
classes [49]. Moreover, the onion like structure of the mixed SLOCC classes
given in the case of three qubits is no longer present [48, 45]. Hence, it is
generally possible to find states that are not contained in the $W$-class and
at the same time are not genuinely multipartite entangled.

In the next section we will push further in this direction and investigate the
discrimination of $W$-class states for an arbitrary number of qubits.

##  5 Discrimination of $W$-class states

In this Section we derive a criterion that allows for a discrimination of
mixed $W$-class states for an arbitrary number $N$ of qubits based on the
second moment $\mathcal{R}^{(2)}$. In particular, we find that

$\displaystyle\mathcal{R}^{(2)}_{\rho}\leq\frac{5-\frac{4}{N}}{3^{N}}=:\chi^{(N)},$
(20)

for all $\rho\in\text{Conv}(\mathcal{W}^{(N)})$, with equality for the pure
$W$-state

$\displaystyle\rvert W_{N}\rangle=\frac{1}{\sqrt{N}}\left(\rvert 10\ldots
0\rangle+\rvert 010\ldots 0\rangle+\ldots+\rvert 0\ldots 01\rangle\right).$
(21)

Hence, any multi-qubit state whose second moment violates the criterion (20)
cannot belong to the mixed $W$-class. While we have addressed such an indirect
characterization of the class $\text{Conv}(\mathcal{W}^{(N)})$ already by
numerical means for up to ten qubits in Ref. [21], we provide here an
analytical proof of the criterion (20) for arbitrary $N$. The basic proof idea
is to maximize $\mathcal{R}^{(2)}$ first over all pure states contained in
$\mathcal{W}^{(N)}$ and then conclude by convexity that the same value also
yields the maximum over the class $\text{Conv}(\mathcal{W}^{(N)})$.
Furthermore, the maximum of $\mathcal{R}_{\rho}^{(2)}$, with
$\rho\in\mathcal{W}^{(N)}$, can be determined via a rather simple geometric
argument, as the value of $\mathcal{R}_{\rho}^{(2)}$ can be associated with a
section of the surface area of a square with side length $1$ (see Fig. 8). For
further details on the proof of criterion (20) see App. C.

Figure 7: Plot of the threshold value $p^{*}$ for the detection of the noisy
GHZ state $\rho_{\text{GHZ}}(p,N)$ in terms of the criteria
$\mathcal{R}^{(2)}\leq\chi^{(N)}$ (purple circles),
$\text{Lin}(\mathcal{R}^{(2)},\mathcal{R}^{(4)})\leq 0$ (yellow squares) and
$\langle W^{(N)}_{\text{GHZ}}\rangle\geq 0$ (cyan diamonds), respectively, as
a function of the number $N$ of qubits.

In order to investigate the performance of criterion (20) we apply it to a
$N$-qubit GHZ state mixed with white noise, i.e.,
$\rho_{\text{GHZ}}(p,N):=p\mathbbm{1}/2^{N}+(1-p)\lvert\text{GHZ}_{N}\rangle\langle\text{GHZ}_{N}\rvert$,
where $\rvert\text{GHZ}_{N}\rangle=\big{(}\rvert 0\rangle^{\otimes N}+\rvert
1\rangle^{\otimes N}\big{)}/\sqrt{2}$. As
$\mathcal{R}^{(2)}_{\mathbbm{1}/2^{N}}=0$, for all $N$, it is easy to see that
the second moment of the noisy GHZ-state $\rho_{\text{GHZ}}(p,N)$ is given by

$\displaystyle\mathcal{R}^{(2)}=\mathcal{R}^{(2)}_{\rvert\text{GHZ}_{N}\rangle}(1-p)^{2},$
(22)

where

$\displaystyle\mathcal{R}^{(2)}_{\rvert\text{GHZ}_{N}\rangle}=\begin{cases}2^{N-1}/3^{N},&N\
\mathrm{odd},\\\ (2^{N-1}+1)/3^{N},&N\ \mathrm{even}.\end{cases}$ (23)

In particular, it has been proven recently in Refs. [20, 30, 29] that the GHZ-
value (23) also yields the maximum of the second moment $\mathcal{R}^{(2)}$.
Hence, the noisy GHZ state $\rho_{\text{GHZ}}(p,N)$ approaches the minimum and
maximum of $\mathcal{R}^{(t)}$ for the maximally mixed state and the GHZ
state, respectively. With Eq. (22) in hand we can easily determine the noise
threshold $p^{*}$ up to which the GHZ state violates the criterion (20):

$\displaystyle
p^{*}=\begin{cases}1-2^{\frac{1}{2}-\frac{N}{2}}\sqrt{5-\frac{4}{N}},&N\
\mathrm{odd},\\\ 1-\frac{\sqrt{10N-8}}{\sqrt{\left(2^{N}+2\right)N}},&N\
\mathrm{even}.\end{cases}$ (24)

As comparison, we also calculate the threshold value $\tilde{p}^{*}$ obtained
from the SLOCC witness found in Refs. [45, 50, 51]. The latter reads
$W^{(N)}_{\text{GHZ}}:=\lambda\mathbbm{1}_{2^{N}}-\lvert\text{GHZ}_{N}\rangle\langle\text{GHZ}_{N}\rvert$,
with $\lambda=3/4$ for $N=3$ and $\lambda=1/2$ for $N\geq 4$, and leads to

$\displaystyle\tilde{p}^{*}=\begin{cases}\frac{2}{7},&N=3,\\\
\frac{2^{N-1}}{2^{N}-1},&N\geq 4.\end{cases}$ (25)

Equations (24) and (25) are plotted in Fig. 7. While the SLOCC witness is more
noise robust for low number of qubits, our criterion (20) performs better for
$N>5$. In this respect, one has to keep in mind that the SLOCC witness
$W^{(N)}_{\text{GHZ}}$ also detects the genuine multipartite entanglement of
the states which partly explains the discrepancies between the performances
given in Eqs. (24) and (25) [1]. Asymptotically, for $N\rightarrow\infty$, we
find that the threshold value $p^{*}$ approaches $1$ which is a direct
consequence of the fact that the bound (20) approaches zero in the same limit.
Equation (24) thus demonstrates the decreasing relevance of the $W$-class in
large multipartite systems.

Hence, if one is interested in the $W$-class alone it is already quiet
exhaustive to consider only the second moment $\mathcal{R}^{(2)}$. We further
investigate this issue by considering a simple generalizaiton of the criterion
(20) which includes also the fourth moment $\mathcal{R}^{(4)}$. The idea is to
define a linear criterion in the $(\mathcal{R}^{(2)},\mathcal{R}^{(4)})$-plane
that passes through the $W$-state $\rvert W_{N}\rangle$ and the bi-seprable
state $\rvert\phi\rangle\rvert W_{N-1}\rangle$. In App. C.2 we prove that this
line indeed forms another possible criterion, denoted as
$\text{Lin}(\mathcal{R}^{(2)},\mathcal{R}^{(4)})\leq 0$, allowing to
discriminate states of the $W$-class $\text{Conv}(\mathcal{W}^{(N)})$. As for
Eq. (20) we further study the performance of this criterion by applying it to
the noisy GHZ-state $\rho_{\text{GHZ}}(p,N)$ and include the results into Fig.
7. We find that in this particular case the generalized criterion
$\text{Lin}(\mathcal{R}^{(2)},\mathcal{R}^{(4)})\leq 0$ does not provide an
advantage over that of Eq. (20). Nevertheless, it allows one to detect other
classes of states which are not in reach of criterion (20). For instance, the
states $\rvert 0\rangle^{\otimes(N-3)}\rvert\text{GHZ}_{3}\rangle$ yield the
second moment $\mathcal{R}^{(2)}=4/3^{N}$ which is smaller or equal to the
bound $\chi^{(N)}$, for all $N\geq 4$ (see Fig. 6(d) for the case $N=4$).
However, it violates the generalized criterion
$\text{Lin}(\mathcal{R}^{(2)},\mathcal{R}^{(4)})\leq 0$, for all $N\geq 3$.
Hence, depending on the class of states under consideration, including the
fourth moment can provide an advantage in detecting states outside of the
W-class $\text{Conv}(\mathcal{W}^{(N)})$.

##  6 Conclusion

In conclusion, we showed how to characterize multipartite entangled states in
a reference frame independent manner through moments of correlation functions
obtained from locally randomized measurements. As the latter correspond
mathematically to uniform averages over polynomials of finite degree, we were
able to exploit the concept of pseudo-random processes in order to evaluate
them exactly for variable system sizes. In particular, we showed, in the case
of qubit systems, how to express the moments in terms of specific examples of
spherical designs and discussed the dependency of the involved number of
measurements on the respective order of the moments. As application of the
introduced framework, we presented an analysis of the structure of different
classes of multipartite entangled states of few qubit systems using the first
two non-vanishing moments. Lastly, we introduced two novel criteria that allow
for a discrimination of $W$-type entangled states in terms of the second and
fourth moment for an arbitrary number of qubits.

The presented framework for the characterization of multipartite entanglement
based on randomized measurements promises advantages, in particular, in the
limit of large system sizes where conventional methods, such as state
tomography, become impractical. In this limit one can attempt to evaluate the
moments (2) approximately using a finite number of randomly measured
correlation functions. As such an approach will consequently lead to
additional statistical errors, it will be subject of future investigations to
develop analytical methods allowing to predict such errors in advance.

## Acknowledgements

We thank Cornelia Spee and Jasmin Meinecke for fruitful discussions. This work
was supported by the ERC (Consolidator Grant 683107/TempoQ), and the DFG. NW
acknowledges support by the QuantERA grant QuICHE and the BMBF. AK
acknowledges support by the Georg H. Endress foundation. The article
processing charge was funded by the Baden-Wuerttemberg Ministry of Science,
Research and Art and the University of Freiburg in the funding programme Open
Access Publishing.

##  Appendix A Qubit unitary $5$-design $SL(2,\mathbb{F}_{5})$

First, we extract a set of generators of $SL(2,\mathbb{F}_{5})$ from the GAP
character library using the package REPSN, as outlined in [39], leading to:

$\displaystyle\begin{pmatrix}-1&0\\\
0&-1\end{pmatrix},\begin{pmatrix}-\omega^{11}-\omega^{14}&\omega^{6}+\omega^{9}\\\
-\omega-\omega^{2}-\omega^{4}-\omega^{7}-\omega^{8}-\omega^{13}&\omega^{11}+\omega^{14}\end{pmatrix},$
$\displaystyle\begin{pmatrix}\omega^{10}&\omega^{11}+\omega^{14}\\\
-\omega^{2}-\omega^{8}&-\omega^{10}\end{pmatrix},\begin{pmatrix}0&\omega^{5}\\\
-\omega^{10}&-\omega^{3}-\omega^{17}\end{pmatrix},$ (26)

with $\omega=e^{i2\pi/15}$. Next, with the matrices (26) we generate the $120$
group elements of $SL(2,\mathbb{F}_{5})=\\{S_{k}|k\in\\{1,\ldots,120\\}\\}$
from which we can then extract an appropriate unitary representation using the
transformation:

$\displaystyle U_{k}=\sqrt{P}S_{k}\sqrt{P}^{-1},\ \ \text{with}\ \
P:=\sum_{k=1}^{120}S_{k}^{\dagger}S_{k}>0.$ (27)

Further on, after eliminating those matrices that are equal up to a global
phase $e^{i\phi}$, we end up with a set of $60$ unitary matrices representing
the corresponding unitary $5$-design.

##  Appendix B Generation of random states

In Sec. 4.2 and 4.3 we used standard forms of multi-qubit states in order to
sample random states from different SLOCC classes. Standard forms are
multipartite generalizations of the Schmidt decomposition, e.g., a three qubit
standard form was introduced in [48]:

$\displaystyle\rvert\Psi_{\text{s.f.}}\rangle=\lambda_{0}\rvert
000\rangle+\lambda_{1}e^{i\phi}\rvert 100\rangle+\lambda_{2}\rvert
101\rangle+\lambda_{3}\rvert 110\rangle+\lambda_{4}\rvert 111\rangle,$ (28)

where $\lambda_{i}\geq 0$, with $i=0,\ldots,4$, and $0\leq\phi\leq\pi$. Any
pure three-qubit state can be represented by Eq. (28) plus an appropriate LU
transformation. For a generalization of Eq. (28) to $N$-partite systems of
local dimensions $d$ we refer the reader to Ref. [52]. More generally,
standard forms allow for a characterization of different SLOCC classes. For
instance, with $\lambda_{4}=\phi=0$ and $\lambda_{i}>0$, for $i=0,1,2,3$, Eq.
(28) becomes a standard form describing the three-qubit $W$-class
$\mathcal{W}^{(3)}$. A similar standard form for the $N$-qubit $W$-class
$\mathcal{W}^{(N)}$ has been presented in Ref. [53]:

$\displaystyle\rvert W_{N}(\vec{x})\rangle=\sqrt{x_{0}}\rvert 0\ldots
0\rangle+\sqrt{x_{1}}\rvert 10\ldots 0\rangle+\ldots+\sqrt{x_{N-1}}\rvert
0\ldots 010\rangle+\sqrt{x_{N}}\rvert 0\ldots 01\rangle.$ (29)

where $\vec{x}=(x_{0},\ldots,x_{N})^{\text{T}}$, and
$|\vec{x}|_{1}={\sum_{i=0}^{N}x_{i}}=1$. Hence, we can sample random pure
states from each of the mentioned classes by drawing randomly the respective
real parameters and apply additional random LU transformations. The latter
step can be skipped in our case as we are interested in the distribution
respective moments $\mathcal{R}^{(t)}$ of the states which are LU invariant.
In Fig. 5(b) and 6(b) of the main text we present randomly sampled three- and
four-qubit $W$-states in the $(\mathcal{R}^{(2)},\mathcal{R}^{(4)})$-plane.

In order to sample mixed quantum states different strategies have been studied
in the literature. One way is to draw a random pure state from a higher
dimensional Hilbert space and subsequently trace over the extended dimensions
[54]. This approach leads to random states which are not homogeneously enough
distributed in the space spanned by the moments $\mathcal{R}^{(2)}$ and
$\mathcal{R}^{(4)}$. Alternatively, we can use the above methods for sampling
random pure states, apply additional random LU transformations and
subsequently mix them with randomly drawn mixing parameters, e.g., for
$\text{Conv}(\mathcal{W}^{(N)})$ we sample
$\rho_{W}=\sum_{\alpha=1}^{2^{N}}p_{\alpha}\left(\bigotimes_{i=1}^{N}U_{i}^{(\alpha)}\right)\lvert
W_{N}(\vec{x}_{\alpha})\rangle\langle
W_{N}(\vec{x}_{\alpha})\rvert\left(\bigotimes_{i=1}^{N}{U_{i}^{(\alpha)}}\right)^{\dagger}$,
with randomly drawn $x_{i}^{(\alpha)}$’s and $p_{\alpha}$’s with
$\sum_{\alpha=1}^{2^{N}}p_{\alpha}=1$. Here, in contrast to the case of pure
states, we cannot simply skip the step of applying additional random LU
transformations as the process of mixing the individual pure states can in
general lead to a larger class of mixed states. However, we found numerically
that mixing of pure standard form states already exhausts quite well the
numerically estimated lower boundaries shown in Fig. 5 and 6. Nevertheless,
the outlined methods to sample mixed states have of course their limits, as
explained in the main text.

##  Appendix C Proof of the criteria presented in Sec. 5

###  C.1 Criterion based on $\mathcal{R}^{(2)}$

In the following, we prove the statement:

$\displaystyle\mathcal{R}^{(2)}_{\rho}\leq\frac{5-\frac{4}{N}}{3^{N}}:=\chi^{(N)},$
(30)

for all $\rho\in\text{Conv}(\mathcal{W}^{(N)})$, with equality for the pure
$W$-state (21).

###### Proof.

In order to prove the bound presented in Eq. (20) we first have to determine
the maximum value of $\mathcal{R}^{(2)}$ for all pure $W$-states, i.e.,
$\chi^{(N)}:=\max_{\rvert\Psi\rangle\in\mathcal{W}^{(N)}}{\mathcal{R}^{(2)}_{\rvert\Psi\rangle}}$.
Given $\chi^{(N)}$, we can subsequently conclude by convexity of
$\mathcal{R}^{(2)}$ that for all
$\rho_{\mathcal{W}}\in\text{Conv}(\mathcal{W}^{(N)})$ the following inequality
holds:

$\displaystyle\mathcal{R}^{(2)}_{\rho_{W}}\leq\sum_{\alpha}p_{\alpha}\mathcal{R}^{(2)}_{\rvert
W_{N}\rangle}\leq\max\limits_{\rho^{\prime}\in\mathcal{W}^{(N)}}\mathcal{R}^{(2)}_{\rvert
W_{N}\rangle}:=\chi^{(N)}.$ (31) $\displaystyle\mathcal{R}^{(2)}_{\rvert
W_{N}\rangle}=\frac{1}{3^{N}}\left[1+8\sum_{i<j;i,j=1}^{N}1/N^{2}\right]=\frac{1}{3^{N}}\left(5-\frac{4}{N}\right).$
(32)

Figure 8: Geometric representation of the value of $\mathcal{R}^{(2)}_{\rvert
W_{N}(\underline{\vec{x}})\rangle}$, for $N=6$, as the blue colored area
within a square of side length $1$ for three different cases: with
$\sum_{i=1}^{N}x_{i}<1$ (a), with $\sum_{i=1}^{N}x_{i}=1$ (b), and with
maximal $\mathcal{R}^{(2)}_{\rvert W_{N}(\underline{\vec{x}})\rangle}$ (c).

In order to calculate the maximum $\chi^{(N)}$, we first evaluate
$\mathcal{R}^{(2)}$ with respect to the standard form (29) and subsequently
optimize it with respect to the involved parameters $(x_{0},\ldots,x_{N})$.
This is sufficient because, on the one hand, for every state
$\rvert\Psi\rangle\in\mathcal{W}^{(N)}$ one can find an appropriate vector
$\vec{x}^{\prime}$ and a LU transformation $U_{1}^{\prime}\otimes\ldots\otimes
U_{N}^{\prime}$ such that
$\rvert\Psi\rangle=(U_{1}^{\prime}\otimes\ldots\otimes U_{N}^{\prime})\rvert
W_{N}(\vec{x}^{*})\rangle$ and, on the other hand, the moments
$\mathcal{R}^{(t)}$ are invariant under LU transformations. Evaluating
$\mathcal{R}^{(t)}$ in terms of Eq. (29) thus yields

$\displaystyle\mathcal{R}^{(2)}_{\rvert
W_{N}(\vec{x})\rangle}=\frac{1}{3^{N}}\left[|\vec{x}|_{2}^{2}+2x_{0}\sum_{i=1}^{N}x_{i}+10\sum_{i<j;i,j=1}^{N}x_{i}x_{j}\right],$
(33)

where $\vec{x}=(x_{0},\ldots,x_{N})^{\text{T}}$, and
$|\vec{x}|_{2}=\sqrt{\sum_{i=0}^{N}x_{i}^{2}}$. Next, we substitute the
constraint $x_{0}=1-\sum_{i=1}^{N}x_{i}$, yielding

$\displaystyle\mathcal{R}^{(2)}_{\rvert
W_{N}(\underline{\vec{x}})\rangle}=\frac{1}{3^{N}}\left[1+8\sum_{i<j;i,j=1}^{N}x_{i}x_{j}\right],$
(34)

where $\underline{\vec{x}}=(x_{1},\ldots,x_{N})^{\text{T}}$. From here it
remains to be checked if the maximum of Eq. (34) is contained in the interior
of the parameter space or lies on its boundaries determined by the constraints
$1\geq\sum_{i=1}^{N}x_{i}$, and $x_{i}\geq 0$, for all $i=1,\ldots,N$. This
question can already be answered in terms of a simple geometric picture. To
see this, we rewrite Eq. (34) as follows:

$\displaystyle\mathcal{R}^{(2)}_{\rvert
W_{N}(\underline{\vec{x}})\rangle}=\frac{1}{3^{N}}\left\\{1+4\left[\left(\sum_{i=1}^{N}x_{i}\right)^{2}-\sum_{i=1}^{N}x_{i}^{2}\right]\right\\},$
(35)

and note that the $x_{i}$-dependent part in the rectangular brackets
determines the size of the blue area of the square with length $1$ presented
in Fig. 8. From this geometric picture one recognizes that in order to reach
the maximum value of $\mathcal{R}^{(2)}_{\rvert
W_{N}(\underline{\vec{x}})\rangle}$, the vector $\underline{\vec{x}}$ has be
on the boundary determined by $\sum_{i=1}^{N}x_{i}=1$. This is clear because
if one takes an exemplary value of $\underline{\vec{x}}$ with
$\sum_{i=1}^{N}x_{i}<1$, as shown in Fig. 8(a), one can simply increase the
value of $\mathcal{R}^{(2)}_{\rvert W_{N}(\underline{\vec{x}})\rangle}$ by
increasing one of the components $x_{i}$. Hence, if all the components $x_{i}$
are non-zero and the maximum is on the boundary $\sum_{i=1}^{N}x_{i}=1$, we
are left with

$\displaystyle\mathcal{R}^{(2)}_{\rvert
W_{N}(\underline{\vec{x}})\rangle}=\frac{1}{3^{N}}\left\\{1+4\left[1-\sum_{i=1}^{N}x_{i}^{2}\right]\right\\},$
(36)

which is maximal whenever $\sum_{i=1}^{N}x_{i}^{2}$ is minimal, i.e., if
$x_{i}=1/N$ for all $i=1,\ldots,N$ (see Fig. 8(c)). This solution corresponds
exactly to the $N$-qubit $W$-state $\rvert W_{N}\rangle$ and leads to the
following value of the second moment

In order to prove that Eq. (32) is the true maximum, it remains to be shown
that $\mathcal{R}^{(2)}$ takes no larger value on the other boundaries. The
remaining boundaries are those where $1=\sum_{i=1}^{N}\underline{x}_{i}$ and
simultaneously $k$ of the $N$ remaining $\underline{x}_{i}$, with
$i=1,\ldots,k$, and all permutations thereof, are equal to zero. In all those
cases we can use a similar geometric argument to find that the maxima are
reached whenever the remaining $N-k$ components read $x_{i}=1/(N-k)$, leading
to

$\displaystyle\mathcal{R}^{(2)}_{\rvert
W_{N}({\underline{\vec{x}}(k)})\rangle}=\frac{1}{3^{N}}\left[1+8\sum_{i<j;i,j=1+k}^{N}\frac{1}{(N-k)^{2}}\right]=\frac{1}{3^{N}}\left(5-\frac{4}{N-k}\right).$
(37)

where $\underline{\vec{x}}(k)$ denote all vectors with $k$ vanishing
components and $N-k$ components equal to $1/(N-k)$. Hence, we have
$\mathcal{R}^{(2)}_{\rvert
W_{N}({\underline{\vec{x}}(k)})\rangle}<\mathcal{R}^{(2)}_{\rvert
W_{N}\rangle}$, for all $k=1,\ldots,N-1$, which shows that the true maximum is
taken by the $W$-state yielding Eq. (32). ∎

###  C.2 Criterion based on a linear combination of $\mathcal{R}^{(2)}$ and
$\mathcal{R}^{(4)}$

We want to show that the linear bound connecting the points $\rvert
W_{N}\rangle$ and $\rvert\phi\rangle\rvert W_{N-1}\rangle$ in the
$(\mathcal{R}^{(2)},\mathcal{R}^{(4)})$-plane provides a bound for the
detection of non-$W$-class states. To do so, it is useful to note the
following relations:

$\displaystyle\mathcal{R}^{(4)}_{\rvert W_{N}\rangle}$
$\displaystyle=\frac{83N^{3}+216N^{2}-176N-96}{27N^{3}5^{N}},$ (38)

and

$\displaystyle\mathcal{R}^{(2)}_{\rvert\phi\rangle\rvert W_{N-1}\rangle}$
$\displaystyle=\frac{1}{3}\mathcal{R}^{(2)}_{\rvert W_{N-1}\rangle},$ (39)
$\displaystyle\mathcal{R}^{(4)}_{\rvert\phi\rangle\rvert W_{N-1}\rangle}$
$\displaystyle=\frac{1}{5}\mathcal{R}^{(4)}_{\rvert W_{N-1}\rangle}.$ (40)

A line passing through the $W$-state $\rvert W_{N}\rangle$ and the bi-seprable
state $\rvert\phi\rangle\rvert W_{N-1}\rangle$ in the
$(\mathcal{R}^{(2)},\mathcal{R}^{(4)})$-plane can then be defined as follows

$\displaystyle\tilde{R}_{4}=m\tilde{R}_{2}-b,$ (41)

with the rescaled moments $\tilde{R}_{t}=(t-1)^{N}\mathcal{R}^{(t)}$, and

$\displaystyle m$
$\displaystyle=\frac{-54N^{4}+196N^{3}-114N^{2}-28N+24}{27N^{2}(N-1)^{2}},$
(42) $\displaystyle b$
$\displaystyle=\frac{353N^{4}-1146N^{3}+829N^{2}+156N-216}{27N^{2}(N-1)^{2}}.$
(43)

Lastly, we need to prove that for all states in the $W$-class
$\text{Conv}(\mathcal{W}^{(N)})$ the expression
$\text{Lin}(\tilde{R}_{2}/3^{N},\tilde{R}_{4}/5^{N}):=\tilde{R}_{4}-m\tilde{R}_{2}-b$
is non-positive. In order to do so, we make use of two facts: first, $m<0$ for
all $N\geq 2$. Thus, the expression $\tilde{R}_{4}-m\tilde{R}_{2}-b$ combines
the two convex functions $\tilde{R}_{2}$ and $\tilde{R}_{4}$ with positive
coefficients, yielding again a convex function. Hence, in order to prove the
criterion it suffices to maximize $\tilde{R}_{4}-m\tilde{R}_{2}-b$ over pure
states in the $W$-class, i.e., $\mathcal{W}^{(N)}$. Second, the expression is
still invariant under local unitary transformations, thus we can optimize over
pure states in the standard form (29) instead. For these states, the second
moment is given in Eq. (34) and the fourth moment reads

$\displaystyle\tilde{R}_{4}$
$\displaystyle=1+\frac{16}{3}\sum_{i<j}x_{i}x_{j}+\frac{128}{3}\sum_{i<j}x_{i}^{2}x_{j}^{2}-\frac{448}{9}\sum_{i<j<k}x_{i}x_{j}x_{k}$
(44) $\displaystyle\phantom{=}+64\sum_{i\neq j,k\wedge
j<k}x_{i}^{2}x_{j}x_{k}+\frac{1664}{9}\sum_{i<j<k<l}x_{i}x_{j}x_{k}x_{l}.$
(45)

where we have already eliminated $x_{0}$ using the normalization constraint
and all sums are assumed to iterate over $1,\ldots,N$. Due to the symmetry of
the target function and the constraints, we assume in the following that
$x_{1}\leq x_{2}\leq\ldots\leq x_{N}$. In order to maximize the expression
$\tilde{R}_{4}-m\tilde{R}_{2}-b$, we first show that for fixed
$x_{2},\ldots,x_{N}$, the optimum lies at the boundary of the set of allowed
values for $x_{1}$, namely ${x}_{1}=\frac{1}{N}$ or ${x}_{1}=0$. To do so, we
closer investigate the partial derivative of $\tilde{R}_{4}-m\tilde{R}_{2}-b$
with respect to $x_{1}$ which shows that it does not vanish in the region
$[0,1/N]$ and thus implies that the optimum lies at the boundary.

In the first case, $x_{1}=x_{2}=\ldots=x_{N}=\frac{1}{N}$ due to the ordering.
In the second case, inserting $x_{1}=0$ into $\tilde{R}_{4}-m\tilde{R}_{2}-b$
yields the same expression for the case of $N-1$ particles, except for the
slope of $m(N)$, which still depends on the larger $N$. However, the
optimization of the remaining variables follows exactly the same steps,
yielding for $x_{2}$ the two candidates $x_{2}=0$ and $x_{2}=\frac{1}{N-1}$.
Thus, we end up with $N+1$ candidate points for the optimum: (1)
$x_{1}=\ldots=x_{N}=\frac{1}{N}$, (2) $x_{1}=0,\
x_{2}=\ldots=x_{N}=\frac{1}{N-1}$, …, (N+1) $x_{1}=\ldots x_{N}=0$. Inserting
all the candidates into $\tilde{R}_{4}-m\tilde{R}_{2}-b$ yields zero for cases
(1) and (2) as expected, as these cases correspond to choosing $\rvert
W_{N}\rangle$ and $\rvert 0\rangle\rvert W_{N-1}\rangle$, and something
negative for the remaining ones, thus proving the bound.

## References

  * [1] O. Gühne and G. Tóth, http://dx.doi.org/10.1016/j.physrep.2009.02.004Phys. Rep. 474, 1 (2009).
  * [2] S. D. Bartlett, T. Rudolph, and R. W. Spekkens, https://doi.org/10.1103/RevModPhys.79.555Rev. Mod. Phys. 79, 555 (2007).
  * [3] R. Ursin, et al., https://doi.org/10.1038/nphys629Nat. Phys. 3, 481 (2007).
  * [4] F. Flamini, N. Spagnolo, and F. Sciarrino, http://dx.doi.org/10.1088/1361-6633/aad5b2 Rep. Prog. Phys. 82, 016001 (2018).
  * [5] J. G. Rarity, P. R. Tapster, P. M. Gorman, and P. Knight, http://dx.doi.org/10.1088/1367-2630/4/1/382New J. Phys. 4, 82 (2002).
  * [6] M. Aspelmeyer, et al., http://dx.doi.org/10.1126/science.1085593Science 301, 621 (2003).
  * [7] P. Villoresi, et al., http://dx.doi.org/10.1088/1367-2630/10/3/033038 New J. Phys. 10, 033038 (2008).
  * [8] C. Bonato, A. Tomaello, V. Da Deppo, G. Naletto, and P. Villoresi, http://dx.doi.org/10.1088/1367-2630/11/4/045017New J. Phys. 11, 045017 (2009).
  * [9] L. Aolita and S. P. Walborn, https://doi.org/10.1103/PhysRevLett.98.100501Phys. Rev. Lett. 98, 100501 (2007).
  * [10] V. D’Ambrosio, E. Nagali, S. P. Walborn, L. Aolita, S. Slussarenko, L. Marrucci, and F. Sciarrino, http://dx.doi.org/10.1038/ncomms1951Nat. Comm. 3, 961 (2012).
  * [11] H. Aschauer, J. Calsamiglia, M. Hein, and H. J. Briegel, https://doi.org/10.26421/QIC4.5Quantum Inf. Comput. 4, 383 (2004).
  * [12] J. I. de Vicente, https://doi.org/10.26421/QIC4.5Quantum Inf. Comput. 7, 624 (2007).
  * [13] J. I. de Vicente, https://doi.org/10.1088/1751-8113/41/6/065309J. Phys. A: Math. Theor. 41, 065309 (2008).
  * [14] J. I. de Vicente and M. Huber, https://doi.org/10.1103/PhysRevA.84.062306 Phys. Rev. A 84, 062306 (2011).
  * [15] P. Badziag, C. Brukner, W. Laskowski, T. Paterek, and M. Żukowski, https://doi.org/10.1103/PhysRevLett.100.140403Phys. Rev. Lett. 100, 140403 (2008).
  * [16] W. Laskowski, M. Markiewicz, T. Paterek, and M. Żukowski, https://doi.org/10.1103/PhysRevA.84.062305Phys. Rev. A 84, 062305 (2011).
  * [17] T. Lawson, A. Pappa, B. Bourdoncle, I. Kerenidis, D. Markham, and E. Diamanti, https://doi.org/10.1103/PhysRevA.90.042336Phys. Rev. A 90, 042336 (2014).
  * [18] C. Klöckl and M. Huber, https://doi.org/10.1103/PhysRevA.91.042339Phys. Rev. A 91, 042339 (2015).
  * [19] M. C. Tran, B. Dakić, F. Arnault, W. Laskowski, and T. Paterek, https://doi.org/10.1103/PhysRevA.92.050301 Phys. Rev. A 92, 050301(R) (2015).
  * [20] M. C. Tran, B. Dakić, W. Laskowski, and T. Paterek, https://doi.org/10.1103/PhysRevA.94.042302Phys. Rev. A 94, 042302 (2016).
  * [21] A. Ketterer, N. Wyderka, O. Gühne, http://dx.doi.org/10.1103/PhysRevLett.122.120505Phys. Rev. Lett. 122, 120505 (2019).
  * [22] M. Krebsbach, https://doi.org/10.6094/UNIFR/150706Bachelor thesis, Albert-Ludwigs-Universität Freiburg (2019).
  * [23] A. Elben, B. Vermersch, M. Dalmonte, J. I. Cirac, and P. Zoller, http://dx.doi.org/10.1103/PhysRevLett.120.050406Phys. Rev. Lett. 120, 050406 (2018).
  * [24] T. Brydges, A. Elben, P. Jurcevic, B. Vermersch, C. Maier, B. P. Lanyon, P. Zoller, R. Blatt, C. F. Roos, http://doi.org/10.1126/science.aau4963Science 364, 260 (2019).
  * [25] A. Elben, B. Vermersch, C. F. Roos, P. Zoller, http://dx.doi.org/10.1103/PhysRevA.99.052323Phys. Rev. A 99, 052323 (2019).
  * [26] L. Knips, J. Dziewior, W. Kłobus, W. Laskowski, T. Paterek, P. J. Shadbolt, H. Weinfurter, and J. D. A. Meinecke, https://doi.org/10.1038/s41534-020-0281-5npj Quantum Information 6, 51 (2020).
  * [27] A. Elben, B. Vermersch, R. van Bijnen, C. Kokail, T. Brydges, C. Maier, M. K. Joshi, R. Blatt, C. F. Roos, and P. Zoller, http://dx.doi.org/10.1103/PhysRevLett.124.010504Phys. Rev. Lett. 124, 010504 (2020).
  * [28] N. Wyderka, F. Huber, and O. Gühne, http://dx.doi.org/10.1103/PhysRevA.97.060101Phys. Rev. A 97, 060101 (2018).
  * [29] N. Wyderka and O. Gühne, https://doi.org/10.1088/1751-8121/ab7f0aJ. Phys. A: Math. Theor. 53, 345302.
  * [30] C. Eltschka and J. Siewert, https://doi.org/10.22331/q-2020-02-10-229 Quantum 4, 229 (2020).
  * [31] M. Idel and M. M. Wolf, https://doi.org/10.1016/j.laa.2014.12.031 Lin. Alg. Appl. 471, 76 (2015).
  * [32] C. Dankert, M.Sc. thesis, University of Waterloo, (2005); also available as e-print https://arxiv.org/abs/quant-ph/0512217quant-ph/0512217.
  * [33] P. D. Seymour, T. Zaslavsky, https://doi.org/10.1016/0001-8708(84)90022-7 Advances in Mathematics 52, 213 (1984).
  * [34] F. G. S. L. Brandão, A. W. Harrow, and M. Horodecki, http://doi.org/10.1103/PhysRevLett.116.170502 Phys. Rev. Lett. 116, 170502 (2016).
  * [35] Y. Nakata, C. Hirche, M. Koashi, A. Winter, https://doi.org/10.1103/PhysRevX.7.021006 Phys. Rev. X 7, 021006 (2017).
  * [36] J. Haferkamp, F. Montealegre-Mora, M. Heinrich, J. Eisert, D. Gross, and I. Roth, https://arxiv.org/abs/2002.09524arXiv:2002.09524.
  * [37] Z. Webb, https://doi.org/10.26421/QIC16.15-16Quantum Inf. Comput. 16, 1379 (2016).
  * [38] H. Zhu, R. Kueng, M. Grassl, and D. Gross, https://arxiv.org/abs/1609.08172arXiv:1609.08172.
  * [39] D. Gross, K. Audenaert, and J. Eisert, https://doi.org/10.1063/1.2716992 J. Math. Phys. 48, 052104 (2007).
  * [40] R. H. Hardin and N. J. A. Sloane, https://doi.org/10.1007/BF02711518Discrete & Computational Geometry 15, 429 (1996).
  * [41] R. Horodecki and M. Horodecki, https://doi.org/10.1103/PhysRevA.54.1838Phys. Rev. A 54, 1838 (1996).
  * [42] M. A. Nielsen and I. L. Chuang, Quantum Computation and Quantum Information https://doi.org/10.1017/CBO9780511976667(Cambridge University Press, New York, 2000).
  * [43] R. Horodecki, P. Horodecki, M. Horodecki, and K. Horodecki, http://doi.org/10.1103/RevModPhys.81.865Rev. Mod. Phys. 81, 865 (2009).
  * [44] W. Dür, G. Vidal, and J. I. Cirac, http://doi.org/10.1103/PhysRevA.62.062314 Phys. Rev. A 62, 062314 (2000).
  * [45] A. Acín, D. Bruß, M. Lewenstein, and A. Sanpera, https://doi.org/10.1103/PhysRevLett.87.040401Phys. Rev. Lett. 87, 040401 (2001).
  * [46] R. Raussendorf and H. J. Briegel, http://doi.org/10.1103/PhysRevLett.86.5188Phys. Rev. Lett. 86, 5188 (2001).
  * [47] R. Cleve, D. Gottesman, and H.-K. Lo, https://doi.org/10.1103/PhysRevLett.83.648Phys. Rev. Lett. 83, 648 (1999).
  * [48] A. Acín, A. Andrianov, L. Costa, E. Jané, J. I. Latorre, and R. Tarrach, https://doi.org/10.1103/PhysRevLett.85.1560Phys. Rev. Lett. 85, 1560 (2000).
  * [49] C. Spee, J. I. de Vicente, and B. Kraus, https://doi.org/10.1063/1.4946895J. Math. Phys. 57, 052201 (2016).
  * [50] T. Bastin, P. Mathonet, and E. Solano, https://doi.org/10.1103/PhysRevA.91.022310Phys. Rev. A 91, 022310 (2015).
  * [51] C. Ritz, C. Spee, and O. Gühne, http://dx.doi.org/10.1088/1751-8121/ab2f54J. Phys. A: Math. Theor. 52, 335302 (2019).
  * [52] H. A. Carteret, A. Higuchi, and A. Sudbery, https://doi.org/10.1063/1.1319516J. Math. Phys. 41, 7932 (2000).
  * [53] S. Kıntaş and S. Turgut, https://doi.org/10.1063/1.3481573J. Math. Phys. 51, 092202 (2010).
  * [54] J. Maziero, https://doi.org/10.1007/s13538-015-0367-2Braz. J. Phys. 45, 575 (2015).
  * [55] D. S. Bernstein, Matrix Mathematics: Theory, Facts, and Formulas - Second Edition, https://press.princeton.edu/books/paperback/9780691140391/matrix-mathematics(Princeton University Press, Princeton, 2009).


# RMPCDMD: Simulations of colloids with coarse-grained hydrodynamics, chemical
reactions and external fields

Pierre de Buyl  Instituut voor Theoretische Fysica, KU Leuven B-3001, Belgium
Mu-Jie Huang  Chemical Physics Theory Group, Department of Chemistry,
University of Toronto, Toronto, Ontario M5S 3H6, Canada  Laurens Deprez
Instituut voor Theoretische Fysica, KU Leuven B-3001, Belgium

###### Abstract

The RMPCDMD software package performs hybrid Molecular Dynamics simulations,
coupling Multiparticle Collision Dynamics to model the solvent and Molecular
Dynamics to model suspended colloids, including hydrodynamics, thermal
fluctuations, and chemically active solvent particles and catalytic colloids.
The main usage of RMPCDMD is the simulation of chemically powered nanomotors,
but other setups are considered: colloids in the presence of a thermal
gradients or forced flows. RMPCDMD is developed in Fortran 2008 with OpenMP
for multithreaded operation and uses the HDF5-based H5MD file format for
storing data. RMPCDMD comes with documentation and a tutorial for the
simulation of chemically powered nanomotors.

## Keywords

Molecular Dynamics ; Colloids ; Hydrodynamics ; Nanomotors ; Fortran ; Python

##  I Overview

###  I.1 Introduction

Molecular Dynamics (MD) is a computational technique to study condensed matter
systems by numerically solving Newton’s equations of motion for a large number
of particles frenkel_smit_2002 . Often, due to the number of constituents in a
system, so-called coarse-grained or multiscale strategies are used to reduce
the computational cost of a simulation.

One simulation strategy for modelling colloids embedded in a solvent is to
couple the Multiparticle Collision Dynamics (MPCD) algorithm
malevanets_kapral_mpcd_1999  for the solvent with MD for the colloids
malevanets_kapral_mpcd_2000 . It differs from traditional MD by the use of a
coarse-grained representation of the solvent’s structure and dynamics and by
the possible occurrence of chemical species transformations in the course of
the simulation. In the latter case, the extension Reactive MPCD (RMPCD)
provides bulk chemical kinetics that reduce to reaction–diffusion equations
rohlf_et_al_rmpcd_2008 . The name RMPCDMD results from the concatenation of
the different algorithms’ acronyms: Reactive Multiparticle Collision Dynamics
and Molecular Dynamics.

The focus of RMPCDMD lies on the modeling of chemically active colloids (e.g.
nanomotors ruckner_kapral_prl_2007 ) or of colloids subjected to external
gradients (e.g. flows prohm_inertial_2012  or temperature gradients
lusebrink_thermophoresis_2012 ). A RMPCDMD simulation consists in
$N_{\mathrm{solvent}}$ solvent particles and $N_{\mathrm{colloids}}$ colloidal
particles with a mass, position, and velocity, placed in a cuboid domain.
Colloid–solvent and colloid–colloid interactions follow Lennard–Jones type
potentials, but there is no interaction potential acting among the solvent
particles. As $N_{\mathrm{solvent}}\gg N_{\mathrm{colloids}}$, this
considerably speeds up a simulation in comparison to a full MD simulation. The
time evolution of colloidal and solvent particles follows Newton’s equations
of motion, which are solved iteratively using the velocity-Verlet algorithm
allen_tildesley_1987 ; frenkel_smit_2002 . The output from a simulation is a
trajectory representing the coordinates of all particles in the course of
time. Typically, solvent data is preaveraged during the simulation as its
storage requires significant amounts of space.

###  I.2 Implementation and architecture

####  I.2.1 Development

RMPCDMD is implemented in Fortran, following the 2008 version of the standard,
and uses OpenMP for multi-core execution. Trajectories are stored in the H5MD
file format h5md_cpc_2014 , and Python programs are provided to illustrate
data analysis. The module fortran_h5md fortran_h5md_github  provides helper
routines to write H5MD files from Fortran. H5MD files are HDF5 files
(https://www.hdfgroup.org/HDF5/) with an internal organization tailored for
molecular simulations h5md_cpc_2014 .

In designing RMPCDMD, we tried to adopt practices for scientific computing
that help to deliver reproducible results and facilitate the extension of the
code and its maintenance. We also made special efforts on the usability of the
software. The design choices are the results of the authors’ experience on
scientific programming, with further motivation by the article on “Best
Practices for Scientific Computing” by Wilson et al best_practices_2014 :

  * •

The algorithms are implemented in Fortran modules in the src/ directory and
shared between the simulations found in programs/. Table 1 summarizes their
content.

  * •

Subroutines and functions are limited in size, so that understanding any
single code unit remains reasonable.

  * •

Data is encapsulated in derived types.

  * •

No global variable is used, to clarify to the reader what data is passed to
subroutines.

  * •

The programming style strives for legibility, conciseness, and efficiency.

  * •

All simulation parameters are given from a text configuration file, there is
no hard-coded setting in the code that would require re-compilation between
runs.

  * •

The code is tracked with the git revision control software git_web  and
published on GitHub (https://github.com/).

  * •

The build procedure and partial testing are automatically executed on the
Travis CI platform (https://travis-ci.org/) for continuous integration.

Module name |  Functionality
---|---
common |  Routines of general use: histogramming, parsing of command-line arguments, timing and minimum-image convention function.
hilbert |  Computation of the compact Hilbert index of a lattice cell and vice-versa.
cell_system |  Cell data structure.
particle_system |  Particles data structure and sorting routine.
particle_system_io |  Data structures to facilitate the storage of thermodynamic and particle data.
interaction |  Computation of Lennard-Jones force and energy.
neighbor_list |  Data structure and routines to update the neighbor list.
mpcd |  Algorithms for MPCD streaming and collision, and for bulk reactions.
md |  Algorithms for velocity-Verlet and rigid body MD.
Table 1: A summary of the functionality provided by the Fortran modules of
RMPCDMD found in the src/ directory.

Instead of building a generic program that could handle the logic of different
geometries, boundary conditions, and colloid types, specific simulation
programs have been implemented. Each program calls the algorithms from the
modules and resembles a statically compiled “simulation script”. The amount of
administrative code (reading parameters, defining storage, etc.) is however
larger than what could be achieved in a higher-level language such as Python.

CMake is used to build RMPCDMD and run automated tests in test/ that cover
some of the base algorithms such as particle sorting or the histogram routine.
The Fortran module fortran_tester fortran_tester_github  provide convenience
routines for testing (e.g. equality assertions). Automated testing using
actual simulations is difficult and very much time intensive, and we rely on
it only for a small part of the code. In simulations without a thermostat or
wall interactions, the total momentum and energy of the system must be
conserved, which is verified for a short run of the single_dimer_pbc
simulation as part of the continuous integration testing. Further validation
of the code is done via comparison of physical results: the velocity
distribution of particles in equilibrium, the viscosity of the fluid, etc.

####  I.2.2 User perspective

For the end-user, the interface comes as a single command-line tool, rmpcdmd,
that drives the simulations and provides access to utility programs. A typical
command is

⬇

rmpcdmd run single_dimer_pbc dimer.parameters dimer.h5 auto

where single_dimer_pbc is the simulation program, dimer.parameters is the
configuration file, dimer.h5 the output datafile, and auto is the setting for
the random-number seed. While the individual programs can be run directly, the
command-line tool rmpcdmd aims to provide a single file to place in the
program “PATH” on the user’s computer and to automate some actions. The
commands available are:

run

execute one of the simulation programs. The setting of the environment
variable OMP_NUM_THREADS, and the start and end time of the run are shown.
Bash’s time provides the “user time” actually used by the program.

seeder

print out a signed 64-bit integer value generated by the computer’s
/dev/urandom device.

plot

display an observable from a simulation file.

timers

print out (or display) the timing information from a simulation file.

A benefit of using separate simulation programs is that the configuration file
for a simulation only contains directly relevant parameters. All simulation
programs share the same command-line syntax, read parameters from a text file
using the ParseText Fortran library ParseText_github , read the seed for the
RNG from the command line and output the trajectory to a H5MD file. An example
configuration file for single_dimer_pbc is listed in Appendix A.

A generic plotting utility is provided and can be invoked as

    rmpcdmd plot dimer.h5 --obs temperature

to plot the temperature, for instance. This plotting program and the other
analysis programs that come with RMPCDMD are written in Python and rely on
NumPy for storing and processing numerical data, SciPy scipy-web  for
numerical algorithms, Matplotlib matplotlib_2007  to produce 2D figures, h5py
collette_python_hdf5_2014  to read HDF5 files and Mayavi mayavi_2011  to
visualize 3D data.

###  I.3 Documentation

A documentation is found in the directory doc/. It is published on the web,
http://lab.pdebuyl.be/rmpcdmd/ and also serves as the homepage of RMPCDMD. The
documentation is built with the Sphinx documentation generator sphinx_doc_web
and is written in reStructuredText. An extension to Sphinx is bundled with
RMPCDMD to provide (i) automatic links from the documentation to the source
code documentation that is annotated with Doxygen
(http://www.stack.nl/~dimitri/doxygen/) and (ii) automatic inclusion of the
Doxygen header of Fortran programs. Specific instructions for the installation
and execution of RMPCDMD are provided, as well as general informations on the
algorithms are presented. The description of the simulation programs is
generated automatically from the beginning of the corresponding .f90 files and
includes all parameters for the simulations.

A tutorial on the simulation of chemically powered nanomotors, based on
RMPCMD, was given at the recent workshop “Modeling of chemically powered
nanomotors” held at the KU Leuven in April 2016
(http://lab.pdebuyl.be/2016-nanomotor/). This tutorial is reproduced in the
documentation of RMPCDMD and represents a unique resource in the field of
nanomotor simulation.

###  I.4 Algorithms

For completeness, we state here exactly what algorithms are implemented in
RMPCDMD, with the corresponding reference to the literature.

RMPCDMD implements the original MPCD collision rule
malevanets_kapral_mpcd_1999  and the Anderson thermostat noguchi_epl_2007 .
Random shifting of the spatial grid is always performed ihle_kroll_srd_2001
to ensure Galilean invariance. Flows are generated by applying a constant
acceleration to solvent particles allahyarov_gompper_mpcd_flows_2002 ;
whitmer_luitjen_2010 . Stick boundary conditions at the wall use the bounce-
back rule and virtual particles during the collision step lamura_mpcd_epl_2001
.

Molecular Dynamics is performed with the velocity-Verlet algorithm
allen_tildesley_1987 ; malevanets_kapral_mpcd_2000 ; frenkel_smit_2002 . The
interaction potential, both for solvent–colloid and colloid–colloid
interactions, are purely repulsive cut-off 12–6 Lennard-Jones potentials of
the form

$V(r)=4\epsilon\left(\left(\frac{\sigma}{r}\right)^{12}-\left(\frac{\sigma}{r}\right)^{6}+\frac{1}{4}\right)\textrm{
for }r<2^{1/6}\sigma$

where the energy parameter $\epsilon$ and the scale parameter $\sigma$ depend
on the species considered. Verlet and cell lists are used to compute the
solvent–colloid forces and energies frenkel_smit_2002 . The sorting of solvent
particles in cells is the same as for the MPCD collision step. The rigid body
dynamics uses RATTLE andersen_rattle_1983  in its iterative implementation for
the Janus particle and exactly for the dimer motor.

Spatial sorting of the solvent particles is based on compact Hilbert indices
hamilton_compact_hilbert_tr . Random numbers are generated with the Threefry
2x64 RNG random123 , allowing each OpenMP thread to generate random numbers
independently. The RNG is provided by the random_module library
random_module_github , where it is implemented in C99. The choices are
inspired by the nano-dimer colberg_nanodimer_web  simulation code, which
implements these algorithms in OpenCL C for large-scale parallel processors.

To obtain a consistent measure of the temperature $T$ even in the presence of
flows, we compute $T$ by averaging the center-of-mass reference frame
temperature over all MPCD cells as

$T=\frac{1}{3N_{c}}\sum_{\xi}\frac{1}{N_{\xi}-1}\sum_{i}m_{i}\left(v_{i}-v_{\xi}\right)^{2}$

where $N_{c}$ is the number of cells, the variable $\xi$ represents a cell,
$N_{\xi}$ the number of particles in a cell, $v_{\xi}$ the center-of-mass
velocity of the cell, and $m_{i}$ and $v_{i}$ the mass and velocity of
particle.

###  I.5 Applications of RMPCDMD

The prototypical use of the MPCD algorithm to study chemically active
particles is the dimer nanomotor, whose model was introduced by Rückner and
Kapral ruckner_kapral_prl_2007 . It consists of two rigidly linked spheres
evolving in a solvent. Chemical reactions are catalyzed by one of the spheres
and create a local concentration gradient, effectively propelling the dimer
along its symmetry axis. This simulation is implemented in the program
single_dimer_pbc.

To facilitate usage of RMPCDMD by newcomers, ready-made simulation
“experiments” are provided in experiments/. There, makefiles allow the user to
simply type

    make simulation

to start a simulation. The terminology that we have chosen reflects our
opinion that computer simulations can be viewed as numerical experiments in
which one prepares a setup and conducts a study to understand a certain
physical phenomenon. The ability to modify a given experimental choice (i.e.
thermostating or not, surface properties, etc.) to assess its effect is
therefore critical. By providing several chemical models and collision rules,
RMPCDMD provides such a framework.

Four programs are provided to analyze the data of this experiment.
plot_velocity.py displays the laboratory reference frame velocity or the
directed velocity of the dimer (with the option \--directed).
plot_histogram.py displays the cylindrical shell histogram of the reaction
product concentration. plot_msd.py displays the mean-squared displacement of
the dimer’s center-of-mass position. This quantity is a practical interest in
comparison to experimental studies howse_et_al_prl_2007 .
view_last_simulation_frame.py displays a 3D representation of the dimer motor
using the Mayavi visualization library mayavi_2011 , either with its full
trajectory as a line (with option \--unwrap) or with the product solvent
particles B (with option \--show-B). The two first analysis programs are
modeled after the results by Rückner and Kapral ruckner_kapral_prl_2007 .

While the aims of single_dimer_pbc are to reproduce Ref.
ruckner_kapral_prl_2007  and to provide a starting point for other
developments, the other simulation programs in RMPCDMD are related to forward-
looking research projects that will benefit from a robust codebase. These
other setups are: poiseuille_flow generates a Poiseuille flow between parallel
plates, with no colloid. chemotactic_cell generates a Poiseuille flow with a
two-species chemical gradient perpendicular to the flow and an embedded
spherical colloid or a dimer nanomotor. single_janus_pbc implements a
composite Janus nanomotor de_buyl_kapral_nanoscale_2013 .
single_sphere_thermo_trap generates a thermal gradient between two plates,
with an embedded spherical colloid.

##  II Availability

Operating System

RMPCDMD has been tested on Linux using the gcc/gfortran and icc/ifort
compilers and on OS X (El Capitan) using gcc/gfortran (installed via
MacPorts).

Programming Language

Fortran 2008 for the simulation code, C99 for the Random Number Generator,
bash and make for the execution and Python (with NumPy, SciPy, matplotlib and
h5py) for analysis.

Additional system requirements

The largest RAM usage in RMPCDMD comes from the solvent coordinates and
requires about $N_{\textrm{solvent}}\times$250 bytes of storage. The example
simulations in the experiments/ directory require between 100MB and 200MB or
RAM.

Dependencies

The CMake and Make tools, and the HDF5 library. For data analysis, NumPy,
SciPy, matplotlib and h5py.

## List of contributors

Pierre de Buyl

Project creator, code, documentation, build system, tutorial.

Laurens Deprez

Code: RATTLE for the dimer, angle-dependent MPCD, boundary conditions and
forcing for flows, Lennard-Jones 9-3 walls, chemotactic cell setup.

Mu-Jie Huang

Tutorial: Janus particle.

Peter Colberg

Initial OpenMP parallelization, continuous integration testing and build
system.

## Software location

### Archive

Name

Zenodo

Persistent identifier

DOI

Licence

BSD 3-clause

Publisher

Pierre de Buyl

Version published

0.1

Date published

TODO

### Code repository

Name

GitHub

Persistent identifier

https://github.com/pdebuyl-lab/RMPCDMD

Licence

BSD 3-clause

Date published

TODO

##  III Reuse potential

Besides nano-dimer colberg_nanodimer_web  (by Peter Colberg, available under
the MIT license), RMPCDMD is the only open-source software to provide support
for chemically-powered nanomotor simulation. RMPCDMD has, in addition to nano-
dimer, support for external fields or rigid-body dynamics for assemblies
larger than the dimer. nano-dimer, on the other hand, is based on OpenCL C and
on the Message Passing Interface for distributed parallel execution, and is
able to handle many-motors simulations that are out of reach of RMPCDMD for
performance reasons. Given the research activity on the topic of active
colloids, we believe that providing a reference implementation of a chemically
active MPCD fluid can serve as a starting point for further modeling, by us or
by other groups.

RMPCMD has not been used for published work yet. It is in use for ongoing
research projects at the Instituut voor Theoretische Fysica of the KU Leuven.

## Acknowledgments

The authors thank Raymond Kapral for discussions and for motivating this line
of research, and Peter Colberg for his contributions to the code, his comments
on the manuscript and for interesting discussions. The authors thank all
authors of the open-source projects that they rely on for their research. PdB
is a postdoctoral fellow of the Research Foundation-Flanders (FWO).

## Competing interests

The authors declare that they have no competing interests.

## References

  * (1) Frenkel, D and Smit, B (2002).  _Understanding molecular simulation: From algorithms to applications_.  2nd edition (Academic Press, San Diego, CA).
  * (2) Malevanets, A and Kapral, R (1999).  Mesoscopic model for solvent dynamics.  _J. Chem. Phys._ 110: 8605.  doi:10.1063/1.478857.
  * (3) Malevanets, A and Kapral, R (2000).  Solute molecular dynamics in a mesoscale solvent.  _J. Chem. Phys._ 112: 7260.  doi:10.1063/1.481289.
  * (4) Rohlf, K, Fraser, S and Kapral, R (2008).  Reactive multiparticle collision dynamics.  _Comput. Phys. Commun._ 179: 132.  doi:10.1016/j.cpc.2008.01.027.
  * (5) Rückner, G and Kapral, R (2007).  Chemically powered nanodimers.  _Phys. Rev. Lett._ 98: 150603.  doi:10.1103/PhysRevLett.98.150603.
  * (6) Prohm, C, Gierlak, M and Stark, H (2012).  Inertial microfluidics with multi-particle collision dynamics.  _Eur. Phys. J. E_ 35(8): 80.  ISSN 1292-8941.  doi:10.1140/epje/i2012-12080-3.  URL http://dx.doi.org/10.1140/epje/i2012-12080-3.
  * (7) Lüsebrink, D, Yang, M and Ripoll, M (2012).  Thermophoresis of colloids by mesoscale simulations.  _J. Phys. Cond. Matt._ 24(28): 284132.  doi:10.1088/0953-8984/24/28/284132.
  * (8) Allen, MP and Tildesley, DJ (1987).  _Computer Simulation of Liquids_ (Clarendon Press, Oxford).
  * (9) de Buyl, P, Colberg, PH and Höfling, F (2014).  H5MD: A structured, efficient, and portable file format for molecular data.  _Comp. Phys. Commun._ 185: 1546–1553.  doi:10.1016/j.cpc.2014.01.018.
  * (10) de Buyl, P (2015-2016).  fortran_h5md.  https://github.com/pdebuyl/fortran_h5md.
  * (11) Wilson, G, Aruliah, DA, Brown, CT, Hong, NPC, Davis, M, Guy, RT, Haddock, SHD, Huff, KD, Mitchell, IM, Plumbley, MD, White, BWEP and Wilson, P (2014).  Best practices for scientific computing.  _PLoS Biol_ 12: e1001745.  doi:10.1371/journal.pbio.1001745.
  * (12) Git.  URL https://git-scm.com/.  https://git-scm.com/.
  * (13) de Buyl, P (2015-2016).  fortran_tester.  https://github.com/pdebuyl/fortran_tester.
  * (14) de Buyl, P (2011-2016).  Parsetext.  https://github.com/pdebuyl/ParseText.
  * (15) Jones, E, Oliphant, T, Peterson, P _et al._ (2001–).  SciPy: Open source scientific tools for Python.  URL http://www.scipy.org/.
  * (16) Hunter, JD (2007).  Matplotlib: A 2d graphics environment.  _Computing In Science & Engineering_ 9(3): 90–95.  doi:10.1109/MCSE.2007.55.
  * (17) Collette, A (2013).  _Python and HDF5_ (O’Reilly).
  * (18) Ramachandran, P and Varoquaux, G (2011).  Mayavi: 3d visualization of scientific data.  _Comput. Sci. Eng._ 13(2): 40–51.  doi:10.1109/MCSE.2011.35.
  * (19) Brandl, G (2007-2016).  Sphinx.  URL http://www.sphinx-doc.org/.  http://www.sphinx-doc.org/.
  * (20) Noguchi, H, Kikuchi, N and Gompper, G (2007).  Particle-based mesoscale hydrodynamic techniques.  _EPL (Europhysics Letters)_ 78(1): 10005.  doi:10.1209/0295-5075/78/10005.
  * (21) Ihle, T and Kroll, DM (2001).  Stochastic rotation dynamics: A galilean-invariant mesoscopic model for fluid flow.  _Phys. Rev. E_ 63: 20201.  doi:10.1103/PhysRevE.63.020201.
  * (22) Allahyarov, E and Gompper, G (2002).  Mesoscopic solvent simulations: Multiparticle-collision dynamics of three-dimensional flows.  _Phys. Rev. E_ 66: 36702.  doi:10.1103/PhysRevE.66.036702.
  * (23) Whitmer, JK and Luijten, E (2010).  Fluid–solid boundary conditions for multiparticle collision dynamics.  _J. Phys. Condens. Matt._ 22(10): 104106.  doi:10.1088/0953-8984/22/10/104106.
  * (24) Lamura, A, Gompper, G, Ihle, T and Kroll, DM (2001).  Multi-particle collision dynamics: Flow around a circular and a square cylinder.  _EPL (Europhysics Letters)_ 56: 319.  doi:10.1209/epl/i2001-00522-9.
  * (25) Andersen, HC (1983).  Rattle: A ”velocity” version of the shake algorithm for molecular dynamics calculations.  _J. Comp. Phys._ 52: 24.  doi:10.1016/0021-9991(83)90014-1.
  * (26) Hamilton, CH (2006).  Compact hilbert indices.  Technical Report CS-2006-07, Dalhousie University.  URL https://www.cs.dal.ca/research/techreports/cs-2006-07.
  * (27) Salmon, JK, Moraes, MA, Dror, RO and Shaw, DE (2011).  Parallel random numbers: As easy as 1, 2, 3.  In _Proceedings of 2011 International Conference for High Performance Computing, Networking, Storage and Analysis_ , SC ’11, pages 16:1–16:12 (ACM, New York, NY, USA).  ISBN 978-1-4503-0771-0.  doi:10.1145/2063384.2063405.  URL http://doi.acm.org/10.1145/2063384.2063405.
  * (28) de Buyl, P (2015-2016).  random_module.  https://github.com/pdebuyl/random_module.
  * (29) Colberg, PH (2013-2015).  nano-dimer.  URL https://colberg.org/nano-dimer/.  https://colberg.org/nano-dimer/.
  * (30) Howse, JR, Jones, RAL, Ryan, AJ, Gough, T, Vafabakhsh, R and Golestanian, R (2007).  Self-motile colloidal particles: From directed propulsion to random walk.  _Phys. Rev. Lett._ 99: 48102.  doi:10.1103/PhysRevLett.99.048102.
  * (31) de Buyl, P and Kapral, R (2013).  Phoretic self-propulsion: a mesoscopic description of reaction dynamics that powers motion.  _Nanoscale_ 5: 1337–1344.  doi:10.1039/c2nr33711h.

##  Appendix A Example configuration file

An example configuration file for the program single_dimer_pbc is shown below.

⬇

# physical parameters

T = 1

L = 32 32 32

rho = 10

tau = 1

probability = 1.0

bulk_rmpcd = T

bulk_rate = 0.1

# simulation parameters

N_MD = 100

N_loop = 500

# interaction parameters

sigma_N = 3.0

sigma_C = 3.0

d = 6.5

epsilon_N = 1.0 0.2

epsilon_C = 1.0 1.0

epsilon_N_N = 1.0

epsilon_N_C = 1.0

epsilon_C_C = 1.0