# Determination of multifractal dimensions of complex networks by means of the
sandbox algorithm

Jin-Long Liu1, Zu-Guo Yu1,2, and Vo Anh2
1Hunan Key Laboratory for Computation and Simulation in Science and
Engineering and Key
Laboratory of Intelligent Computing and Information Processing of Ministry of
Education,
Xiangtan University, Xiangtan, Hunan 411105, China.
2School of Mathematical Sciences, Queensland University of Technology, GPO Box
2434,
Brisbane, Q4001, Australia. Corresponding author, email<EMAIL_ADDRESS>
###### Abstract

Complex networks have attracted much attention in diverse areas of science and
technology. Multifractal analysis (MFA) is a useful way to systematically
describe the spatial heterogeneity of both theoretical and experimental
fractal patterns. In this paper, we employ the sandbox (SB) algorithm proposed
by Tél et al. (Physica A, 159 (1989) 155-166), for MFA of complex networks.
First we compare the SB algorithm with two existing algorithms of MFA for
complex networks: the compact-box-burning (CBB) algorithm proposed by Furuya
and Yakubo (Phys. Rev. E, 84 (2011) 036118), and the improved box-counting
(BC) algorithm proposed by Li et al. (J. Stat. Mech.: Theor. Exp., 2014 (2014)
P02020) by calculating the mass exponents $\tau(q)$ of some deterministic
model networks. We make a detailed comparison between the numerical and
theoretical results of these model networks. The comparison results show that
the SB algorithm is the most effective and feasible algorithm to calculate the
mass exponents $\tau(q)$ and to explore the multifractal behavior of complex
networks. Then we apply the SB algorithm to study the multifractal property of
some classic model networks, such as scale-free networks, small-world
networks, and random networks. Our results show that multifractality exists in
scale-free networks, that of small-world networks is not obvious, and it
almost does not exist in random networks.

Key words: Complex network; multifractal analysis; sandbox algorithm; box-
counting algorithm

##  1 Introduction

Many studies have shown that complex networks play an important role in
characterizing complicated dynamic systems in nature and society [1]. This is
because the nodes of a complex network represent the elements, and the edges
represent and simplify the complexity of their interactions so that we can
better focus on the topological relation between two elements in a complex
system. Recently, complex networks have attracted the attention of a lot of
researchers from different fields. Based on the self-similarity of fractal
geometry [2, 3, 4], Song et al. [1] generalized the box-counting algorithm and
used it in the field of complex networks. They found that many complex
networks are self-similar under certain length-scales. The fractal dimension
has been widely used to characterize complex fractal sets [2, 3, 4]. Because
the metric on graphs is not the same as the Euclidian metric on Euclidian
spaces, the box-counting algorithms to calculate the fractal dimension of a
network is much more complicated than the traditional box-counting algorithm
for fractal sets in Euclidian spaces. Song et al. [5] developed some
algorithms to calculate the fractal dimension of complex networks. Then Kim et
al. [6, 7] proposed an improved algorithm to investigate the skeleton of
networks and the fractal scaling in scale-free networks. Zhou et al. [8]
proposed an algorithm based on the edge-covering box counting to explore self-
similarity of complex cellular networks. Later on, a ball-covering approach
[9] and an approach defined by the scaling property of the volume [10, 11]
were proposed for fractal dimensions of complex networks. The features of
topology and statistics [13, 12], the fractality and percolation transition
[14], fractal transition [15] in complex networks, and properties of a scale-
free Koch networks [16, 17, 18, 19, 20] have turned out to be hot topics in
recent years.

As a generalization of fractal analysis, the tool of multifractal analysis
(MFA) may have a better performance on characterizing the complexity of
complex networks in real world. MFA has been widely applied in a variety of
fields such as financial modelling [21, 22], biological systems [23, 24, 25,
26, 27, 28, 29, 30, 31, 32], and geophysical data analysis [33, 34, 35, 36,
37, 38, 39, 40, 41, 42].

In recent years, MFA also has been successfully used in complex networks and
seems more powerful than fractal analysis. Lee and Jung [43] found that MFA is
the best tool to describe the probability distribution of the clustering
coefficient of a complex network. Some algorithms have been proposed to
calculate the mass exponents $\tau(q)$ and to study the multifractal
properties of complex networks [44, 45, 46, 47]. Based on the compact-box-
burning algorithm for fractal analysis of complex networks which is introduced
by Song et al. [5], Furuya and Yakubo [44] proposed a compact-box-burning
(CBB) algorithm for MFA of complex networks and applied it to show that some
networks have multifractal structures. Wang et al. [45] proposed a modified
fixed-size box-counting algorithm to detect the multifractal behavior of some
theoretical and real networks. Li et al. [46] improved the modified fixed-size
box-counting algorithm [45] and used it to investigate the multifractal
properties of a family of fractal networks introduced by Gallos et al. [48].
We call the algorithm in Ref. [46] the improved BC algorithm. Recently, we
adopted the improved BC to study the multifractal properties of the recurrence
networks constructed from fractional Brownian motions [47].

In order to easily obtain the generalized fractal dimensions of real data, Tél
et al. [49] introduced a sandbox algorithm which is originated from the box-
counting algorithm [50]. They [49] pointed out that the sandbox algorithm
gives a better estimation of the generalized fractal dimensions in practical
applications. So far, the sandbox algorithm also has been widely applied in
many fields. For example, Yu et al. [26] used it to perform MFA on the
measures based on the chaos game representation of protein sequences from
complete genomes.

In this article, we employ the sandbox (SB) algorithm proposed by Tél et al.
[49] for MFA of complex networks. First we compare the SB algorithm with the
CBB and improved BC algorithms for MFA of complex networks in detail by
calculating the mass exponents $\tau(q)$ of some deterministic model networks.
We make a detailed comparison between the numerical and theoretical results of
these model networks. Then we apply the SB algorithm to study the multifractal
property of some classic model networks, such as scale-free networks, small-
world networks, and random networks.

##  2 Sandbox algorithm for multifractal analysis of complex networks

It is well known that the fixed-size box-covering algorithm [50] is one of the
most common and important algorithms for multifractal analysis. For a given
measures $\mu$ with support set $E$ in a metric space, we consider the
following partition sum

$Z_{\epsilon}(q)=\sum_{\mu(B)\neq 0}[\mu(B)]^{q},$ (1)

$q\in R$, where the sum runs over all different nonempty boxes $B$ of a given
size $\epsilon$ in a box covering of the support set $E$. From the definition
above, we can easily obtain $Z_{\epsilon}(q)\geq 0$ and $Z_{\epsilon}(0)=1$.
The mass exponents $\tau(q)$ of the measure $\mu$ can be defined as

$\tau(q)=\lim_{\epsilon\rightarrow 0}\frac{\ln Z_{\epsilon}(q)}{\ln\epsilon}.$
(2)

Then the generalized fractal dimensions $D(q)$ of the measure $\mu$ are
defined as

$D_{q}=\frac{\tau(q)}{q-1},~{}~{}~{}~{}~{}~{}\textrm{for}~{}~{}q\neq 1,$ (3)

and

$D_{q}=\lim_{\epsilon\rightarrow
0}\frac{Z_{1,\epsilon}}{\ln\epsilon},~{}~{}~{}\textrm{for}~{}~{}q=1,$ (4)

where $Z_{1,\epsilon}=\sum_{\mu(B)\neq 0}\mu(B)\ln\mu(B)$. The linear
regression of $[\ln Z_{\epsilon}(q)]/(q-1)$ against $\ln\epsilon$ for $q\neq
1$ gives a numerical estimation of the generalized fractal dimensions $D_{q}$,
and similarly a linear regression of $Z_{1,\epsilon}$ against $\ln\epsilon$
for $q=1$. In particular, $D_{0}$ is the box-counting dimension (or fractal
dimension), $D_{1}$ is the information dimension, and $D_{2}$ is the
correlation dimension.

In complex network, the measure $\mu$ of each box can be defined as the ratio
of the number of nodes covered by the box and the total number of nodes in the
entire network. In addition, we can determine the multifractality of complex
network by the shape of $\tau(q)$ or $D(q)$ curve. If $D(q)$ is a constant or
$\tau(q)$ is a straight line, the object is monofractal; on the other hand, if
$D(q)$ or $\tau(q)$ is convex, the object is multifractal.

Before we use the following SB algorithm to perform MFA of a network, we need
to apply the Floyd’s algorithm [51] of Matlab-BGL toolbox [52] to calculate
the shortest-path distance matrix $D$ of this network according to its
adjacency matrix $A$.

The sandbox algorithm proposed by Tél et al. [49] is an extension of the box-
counting algorithm [50]. The main idea of this sandbox algorithm is that we
can randomly select a point on the fractal object as the center of a sandbox
and then count the number of points in the sandbox. The generalized fractal
dimensions $D(q)$ are defined as

$D_{q}=\lim_{r\rightarrow
0}\frac{\ln\langle[M(r)/M(0)]^{q-1}\rangle}{\ln(r/d)}\frac{1}{q-1},~{}~{}~{}~{}q\in
R,$ (5)

where $M(r)$ is the number of points in a sandbox with a radius of $r$, $M(0)$
is the total number of points in the fractal object. The brackets
$\langle\cdot\rangle$ mean to take statistical average over randomly chosen
centers of the sandboxes. As a matter of fact, the above equation can be
rewritten as

$\ln(\langle[M(r)]^{q-1}\rangle)\ \propto\ D(q)(q-1)\ln(r/d)+(q-1)\ln(M_{0}).$
(6)

So, in practice, we often estimate numerically the generalized fractal
dimensions $D(q)$ by performing a linear regression of
$\ln(\langle[M(r)]^{q-1}\rangle)$ against $(q-1)\ln(r/d)$; and estimate
numerically the mass exponents $\tau(q)$ by performing a linear regression of
$\ln(\langle[M(r)]^{q-1}\rangle)$ against $\ln(r/d)$. In a complex network, we
can randomly choose a node of a network as the center of a sandbox. $M(r)$ and
$M(0)$ represent the number of nodes in the sandbox of radius $r$ and the size
of the network, respectively. The SB algorithm for MFA of complex networks can
be described as follows.

  1. (i)

Initially, make sure all nodes in the entire network are not selected as a
center of a sandbox.

  2. (ii)

Set the radius $r$ of the sandbox which will be used to cover the nodes in the
range $r\in[1,d]$, where $d$ is the diameter of the network.

  3. (iii)

Rearrange the nodes of the entire network into random order. More
specifically, in a random order, nodes which will be selected as the center of
a sandbox are randomly arrayed.

  4. (iv)

According to the size $N$ of networks, choose the first 1000 nodes in a random
order as the center of 1000 sandboxes, then search all the neighbor nodes by
radius $r$ from the center of each sandbox.

  5. (v)

Count the number of nodes in each sandbox of radius $r$, denote the number of
nodes in each sandbox as $M(r)$.

  6. (vi)

Calculate the statistical average $\langle[M(r)]^{q-1}\rangle$ of
$[M(r)]^{q-1}$ over all 1000 sandboxes of radius $r$.

  7. (vii)

For different values of $r$, repeat steps (ii) to (vi) to calculate the
statistical average $\langle[M(r)]^{q-1}\rangle$ and then use
$\langle[M(r)]^{q-1}\rangle$ for linear regression.

We need to choose an appropriate range of $r\in[r_{min},r_{max}]$, then
calculate the generalized fractal dimensions $D(q)$ and the mass exponents
$\tau(q)$ in this scaling range. In our calculation, we perform a linear
regression of $\ln(\langle[M(r)]^{q-1}\rangle)$ against $\ln(r/d)$ and then
choose the slope as an approximation of the mass exponents $\tau(q)$ (the
process for estimating the generalized fractal dimensions $D(q)$ is similar).

For the improved BC and CBB algorithms, we need to cover the entire network by
repeating a large number of same steps and then to find the minimum possible
number of boxes by performing many realizations. Then we can choose an
appropriate range of $r\in[r_{min},r_{max}]$ to calculate the mass exponents
$\tau(q)$ by performing a linear regression of $\ln Z_{r/d}(q)$ against
$\ln(r/d)$. In the process of finding a covering of the network, the two
existing algorithms require that each node in the network cannot be covered by
more than one box at the same time. In addition, for the CBB algorithm, we
have to randomly select many nodes $p$ from the candidate set $C$ to form a
compact box and then repeat these steps until the entire network is covered.
Because of these limitations, these two algorithms must take a large amount of
CPU time and memory resources to record information of the nodes which have
been covered in the previous steps. From the above descriptions of SB
algorithm, however, we find that a big difference from the improved BC and CBB
algorithms is that SB algorithm only requires to randomly choose some nodes as
the center of a sandbox and then to count the number of nodes in each sandbox,
hence we don’t need so many sandboxes to cover the entire network. Moreover,
we focus on the number of nodes of each sandbox, so we also don’t need to know
whether or not the nodes in the sandbox have been covered by other sandboxes.
For the SB algorithm, therefore, we only consume a very small amount of CPU
time and memory resources. In this sense, the SB algorithm can be considered
to be the most effective and feasible algorithm for MFA of complex networks.

##  3 Algorithm comparison

Now we compare the SB algorithm with the CBB and improved BC algorithms for
MFA of complex networks [44, 46] in detail by calculating the mass exponents
$\tau(q)$ of some deterministic model networks. We make a detailed comparison
between the numerical and theoretical results of these model networks.

###  3.1 The model networks

In the past two decades, many network models have been introduced to study and
simulate the topological and fractal properties and the growth mechanism of
many complex dynamical systems in real world. Watts and Strogatz [53], Newman
and Watts [54] proposed the WS and NW small-world network models to explain
the small-world character of many real complex networks respectively. In order
to reveal the generating principle of power law distributions, Barabási et al.
[55] proposed a scale-free network model (BA model) based on the growth and
preferential attachment characteristics of real networks.

In 2002, Dorogovtsev et al. [56] introduced the simple deterministic graphs to
model scale-free networks. They pointed out that the family of deterministic
networks (DGM networks) are pseudo-fractals. In order to understand the self-
similarity of complex networks better, Rozenfeld et al. [57] generalized the
DGM network to a family of scale-free networks, namely $(u,v)$-flowers and
$(u,v)$-trees. The DGM network is a special case of the $(1,v)$-flowers. For
$u=1$, the networks are self-similar only in the weak sense. For $v\geq u>1$,
the networks are self-similar and possess the well-defined fractal dimension.
In 2006, Song et al. [58] proposed the minimal model to study the evolved law
of complex networks and then simulate the emergence of self-similarity and
small-world properties of these networks. Later on, Gallos et al. [48]
introduced a generalized version of the minimal model proposed by Song et al.
[58].

Based on the above network models, some researchers have studied analytically
and numerically the fractal and multifractal properties of networks generated
from these models. Song et al. [58] and Gallos et al. [48] gave an analytical
formula to calculate the fractal dimensions of the minimal model and its
generalized version respectively. Then Rozenfeld et al. [57] also put forward
a theoretical framework for computing the degree exponent and the fractal
dimension of the $(u,v)$-flower network. In addition, Furuya and Yakubo [44]
analytically and numerically investigated the multifractal properties of
several deterministic, stochastic, and real-world fractal scale-free networks.
They gave an analytical formula of the mass exponents $\tau(q)$ for some class
of fractal scale-free model networks by a mean-field approximation. They
define the mass exponents $\tau(q)$ of a complex network as

$\tau(q)=\left\\{\begin{array}[]{cc}(q-1)D_{f},~{}~{}~{}~{}~{}q<\gamma-1,\\\
qD_{f}\frac{\gamma-2}{\gamma-1},~{}~{}~{}~{}~{}~{}~{}q\geq\gamma-1,\end{array}\right.$
(7)

where $D_{f}$ is the fractal dimension, $\gamma$ is the degree exponent of the
complex network.

For the $(u,v)$-flower, we start with a cycle graph consisting of $u+v$ nodes
and $u+v$ links. Then we can obtain the $(u,v)$-flower network of generation
$n+1$ by replacing each link in $n$th generation network by two parallel paths
with length (number of links of the path) $u$ and $v$ respectively. In Fig. 1,
we show how to construct the $(2,2)$-flower network as an example. Rozenfeld
et al. [57] gave the analytical formulas of degree exponent $\gamma$ and
fractal dimension $D_{f}$ of $(u,v)$-flower. The degree exponent $\gamma$ of
the $n$th generation $(u,v)$-flower network is given by

$\gamma=1+\frac{\ln(u+v)}{\ln 2},$ (8)

and the fractal dimension $D_{f}$ is presented by

$D_{f}=\frac{\ln(u+v)}{\ln u},~{}~{}u>1,v\geq u.$ (9)

By substituting Eqs. (8) and (9) into Eq. (7), the mass exponents $\tau(q)$ of
the $(u,v)$-flower network can be written as

$\tau(q)=\left\\{\begin{array}[]{cc}(q-1)\frac{\ln(u+v)}{\ln
u},~{}~{}~{}~{}~{}q<\frac{\ln(u+v)}{\ln 2},\\\ q\frac{\ln((u+v)/2)}{\ln
u},~{}~{}~{}~{}~{}~{}~{}~{}q\geq\frac{\ln(u+v)}{\ln 2}.\end{array}\right.$
(10)

Figure 1: Construction of a $(u,v)$-flower network. Example of the
$(u,v)$-flower network model of generations $n=1,2,3$ with parameters $u=2$
and $v=2$.

The minimal model proposed by Song et al. [58] is a probabilistic combination
of two different connectivity modes: Mode I with probability $e$ and Mode II
with probability $1-e$. Mode I means that all the old connections generated in
the previous generation are remained; Mode II means that we remove all the old
connections generated in the previous generation and add a new edge to connect
two new generated nodes. Before using Modes I and II, we attach $m$ new nodes
to each endpoint of each edge $l$ in the network of the current generation. A
remarkable advantage of this stochastic combination of the two different
growth modes is that its level of fractality can be controlled by the
probability $e$. Song et al. [58] pointed out that the minimal model is a pure
fractal network when the probability $e=0$, and a pure small-world network
when $e=1$. In this paper, we only consider the minimal model with probability
$e=0$. We start with a star structure as in Ref. [58]. Then we obtain the
minimal model of generation $n+1$ by adding $mk_{i}$ new nodes to each node
$i$ with degree $k_{i}$ of generation $n$, where $m$ is a given parameter. In
addition, we adopt the growth Mode II to replace all the old connections in
the previous generation. In Fig. 2, we show how to construct the pure fractal
network with parameters $m=2$ and $e=0$ as an example. As pointed out by Song
et al. [58], the degree exponent $\gamma$ is

$\gamma=1+\frac{\ln(2m+1)}{\ln m},$ (11)

and the fractal dimension $D_{f}$ is

$D_{f}=\frac{\ln(2m+1)}{\ln 3}.$ (12)

By substituting Eqs. (11) and (12) into Eq. (7), the mass exponents $\tau(q)$
of this minimal model can be written as

$\tau(q)=\left\\{\begin{array}[]{cc}(q-1)\frac{\ln(2m+1)}{\ln
3},~{}~{}~{}~{}~{}q<\frac{\ln(2m+1)}{\ln m},\\\ q\frac{\ln((2m+1)/m)}{\ln
3},~{}~{}~{}~{}~{}~{}~{}q\geq\frac{\ln(2m+1)}{\ln m}.\end{array}\right.$ (13)

Figure 2: Construction of a pure fractal network. Example of the minimal model
network of generations $n=0,1$ with parameters $m=2$ and $e=0$.

Based on the minimal model, Gallos et al. [48] proposed a generalized version
of this model. Here we start with two nodes and one edge between them as in
Ref. [48]. Then we obtain the network of next generation by attaching $m$ new
nodes to each endpoint of each edge $l$ in the network of the current
generation. With probability $e$, each edge $l$ of the current generation is
remained and $x-1$ $(x\leq m)$ new edges are added to connect pairs of new
nodes attached to the endpoints of $l$. Otherwise, we remove each edge $l$ in
the network of the current generation and add $x$ $(x\leq m)$ new edges to
connect pairs of nodes attached to the endpoints of $l$. In this paper, we
only consider the network with probability $e=0$. In Fig. 3, we show how to
construct the generalized network with parameters $m=2$, $x=2$, and $e=0$, as
an example. In the case of probability $e=0$, from Refs. [46, 14], the degree
exponent $\gamma$ is

$\gamma=1+\frac{\ln(2m+x)}{\ln m},$ (14)

and the fractal dimension $D_{f}$ is

$D_{f}=\frac{\ln(2m+x)}{\ln 3}.$ (15)

By substituting Eqs. (14) and (15) into Eq. (7), we can obtain the mass
exponents $\tau(q)$ of the generalized version of the minimal model, which can
be written as

$\tau(q)=\left\\{\begin{array}[]{cc}(q-1)\frac{\ln(2m+x)}{\ln
3},~{}~{}~{}~{}~{}q<\frac{\ln(2m+x)}{\ln m},\\\ q\frac{\ln((2m+x)/m)}{\ln
3},~{}~{}~{}~{}~{}~{}~{}q\geq\frac{\ln(2m+x)}{\ln m}.\end{array}\right.$ (16)

Figure 3: Construction of a generalized network of the minimal model. Example
of the generalized network of generations $n=0,1,2$ with parameters $m=2$,
$x=2$, and $e=0$.

In the following, we generate networks using the three models and numerically
study their multifractality by the SB, CBB and improved BC algorithms (in
Section II). Then we give a detailed comparison between the three algorithms
based on the numerical results and theoretical formulas of the mass exponents
$\tau(q)$ for these model networks.

###  3.2 Comparison results

In this work, we set the range of the $q$ values from $-10$ to $10$ with a
step of $1/3$. In order to compare with the results in Ref. [44], we generated
the $(u,v)$-flower network with $u=2$ and $v=2$. Considering the limitation of
the computational capacity of our computer, we only constructed the 7th
generation $(u,v)$-flower network. An important step of calculating the mass
exponents $\tau(q)$ and the generalized fractal dimensions $D(q)$ is to obtain
the appropriate range of $r$ (i.e. $r\in[r_{min},r_{max}]$). As an example, we
show the linear regressions of the $\ln(\langle[M(r)]^{q-1}\rangle)$ vs
$\ln(r/d)$ in the SB algorithm for the 7th generation $(u,v)$-flower network
with $u=2$ and $v=2$ in Fig. 4. We selected the range of $r$ as $[2,20]$ to
fit these data points. From Fig. 4, we can observe the apparent power law
behaviors for the 7th generation $(u,v)$-flower network with $u=2$ and $v=2$.
So we selected the linear fit scaling range of $r$ as $[2,20]$ to calculate
the mass exponents $\tau(q)$. In Fig. 5, we show the mass exponents $\tau(q)$
of the $(u,v)$-flower network calculated by the SB, CBB and improved BC
algorithms. From Fig. 5, we can see that the numerical results obtained by the
three algorithms are consistent with the theoretical results.

For the minimal model network, we started with a star structure of 5 nodes as
in Ref. [58] and then generated the 5th generation network with $m=2$ and
$e=0$. We calculated the mass exponents $\tau(q)$ by the three algorithms. The
numerical results are shown in Fig. 6. From Fig. 6, we can see that the
numerical results obtained by the three algorithms agree well with the
theoretical results.

In addition, we also generated the generalized version of the minimal model
with $m=2$, $x=2$, and $e=0$. Here we only constructed the 5th generation of
the generalized network. The numerical results are shown in Fig. 7. From Fig.
7, we can see that the numerical results obtained by the three algorithms have
good agreement with the theoretical results.

Figure 4: Linear regressions for calculating the mass exponents $\tau(q)$ of
the 7th generation $(u,v)$-flower network with $u=2$ and $v=2$. The result is
calculated by the SB algorithm.

Figure 5: (Color online) Mass exponents $\tau(q)$ for the 7th generation
$(u,v)$-flower network with $u=2$ and $v=2$. Solid line (black line)
represents the mass exponents $\tau(q)$ given by Eq. (10). The curves
indicated by symbols (circles, squares, and diamonds) represent the numerical
estimation of the mass exponent $\tau(q)$ calculated by the three algorithms
for MFA of complex networks, respectively.

Figure 6: (Color online) Mass exponents $\tau(q)$ for the 5th generation
minimal model network with $m=2$ and $e=0$. Solid line (black line) represents
the mass exponent $\tau(q)$ given by Eq. (13). The curves indicated by symbols
(circles, squares, and diamonds) represent the numerical estimation of the
mass exponent $\tau(q)$ calculated by the three algorithms for MFA of complex
networks, respectively.

Figure 7: (Color online) Mass exponents $\tau(q)$ for the 5th generation
network from generalized version of the minimal model with $m=2$, $x=2$, and
$e=0$. Solid line (black line) represents the mass exponent $\tau(q)$ given by
Eq. (16). The curves indicated by symbols (circles, squares, and diamonds)
represent the numerical estimation of the mass exponent $\tau(q)$ calculated
by the three algorithms for MFA of complex networks, respectively.

It is hard to evaluate the performance of the three algorithms only according
to the above three figures. In order to further quantify the performance of
these algorithms, we introduce the relative standard error as in Ref. [29].
Based on the relative standard error $E$, we can determine the goodness of the
numerical results of the mass exponents $\tau(q)$ obtained from the three
algorithms compared with the analytical results for the three deterministic
model networks. The relative standard error $E$ can be defined as

$E=\frac{E_{1}}{E_{2}},$ (17)

where

$E_{1}=\sqrt{\frac{1}{61}\sum_{q}(\tau_{t}(q)-\tau_{n}(q))^{2}},$ (18)

and

$E_{2}=\sqrt{\frac{1}{61}\sum_{q}(\tau_{t}(q)-\bar{\tau_{t}})^{2}},$ (19)

the $\tau_{t}(q)$ and $\tau_{n}(q)$ represent the analytical and numerical
results of the mass exponents $\tau(q)$ respectively; $\bar{\tau_{t}}$ is the
average of the $\tau_{t}(q)$. The goodness of fit is indicated by the result
$E<1$ [29]. We summarize the corresponding relative standard error between the
analytical and numerical results of the mass exponents $\tau(q)$ in Table I.
From Table I, we see that the relative standard errors for these three methods
are all rather small. This result indicates all these three algorithms can
give correct numerical results.

In addition, we compare the consumed CPU time of the three algorithms for MFA
of the networks generated from the three network models. The results are given
in Table II. From Table II, we can CBB algorithm takes a substantial amount of
computation time and memory resources, while SB algorithm consumes the least
amount of CPU time and memory resources. This is to say that the SB algorithm
has an overwhelming advantage in consuming CPU time and memory resources.
Therefore the SB algorithm can be considered as the most effective, feasible,
and accurate algorithm to calculate the mass exponents $\tau(q)$ and then to
study the multifractal properties of complex networks among the three
algorithms.

Table 1: The relative standard error $E$ of the three algorithms for MFA of the networks generated from the three network models. |  $(u,v)$-flower | The minimal model | The generalized model
---|---|---|---
improved box-counting (BC) [46] | 0.01406849 | 0.00533030 | 0.02067627
Compact-box-burning (CBB) [44] | 0.01567521 | 0.02449339 | 0.02245775
Sandbox (SB) | 0.01408235 | 0.00623890 | 0.01299720
Table 2: The CPU time of the three algorithms for MFA of the networks generated from the three network models(unit: s). |  $(u,v)$-flower | The minimal model | The generalized model
---|---|---|---
improved box-counting (BC) [46] | 4069 | 8568 | 1299
Compact-box-burning (CBB) [44] | 817891 | 1003087 | 244206
Sandbox (SB) | 95 | 135 | 44

##  4 Applications

Wang et al. [45] applied the modified fixed-size box-counting algorithm to
study the multifractal properties of some theoretical model networks and real
networks, such as scale-free networks, small-world networks, and random
networks. All of these complex networks have been widely used in various
studies. In this Section, we want to adopt the SB algorithm to detect the
multifractal behavior of these networks.

###  4.1 Scale-free networks

Based on the growth and preferential attachment characteristics of many
complex networks in real world, Barabási et al. [55] proposed a BA model to
explain the generating mechanism of power law distributions. In this paper, we
use the BA model to generate scale-free networks and then investigate their
multifractality.

Here, we start with an initial network with $m=3$ nodes, and its three nodes
are connected each other. At each step, we add one node to this initial
network. Then this new node is connected to $m_{0}=1$ existing nodes with
probability $e$. We denote the probability that the new node is connected to
node $i$ as $e_{i}$. The probability $e_{i}$ is defined as
$e_{i}=\frac{k_{i}}{\sum_{j}k_{j}}$, where $k_{i}$ is the degree of node $i$.

In this paper, we respectively generated 100 scale-free networks with 6000
nodes, 8000 nodes, and 10000 nodes. The SB algorithm is then used to calculate
their average mass exponents $\langle\tau(q)\rangle$ and average generalized
fractal dimensions $\langle D(q)\rangle$, which are shown in Fig. 8. From Fig.
8, we find that the $\langle\tau(q)\rangle$ and $\langle D(q)\rangle$ curves
of scale-free networks are convex. So multifractality exists in these scale-
free networks.

Figure 8: (Color online) Average mass exponents $\langle\tau(q)\rangle$ (upper
panel) and average generalized fractal dimensions $\langle D(q)\rangle$ (lower
panel) for the scale-free networks derived from the SB algorithm. Here the
average is calculated from 100 realizations.

###  4.2 Small-world networks

Based on the random rewiring procedure, Watts and Strogatz [53] introduced the
WS small-world network model which is a transition from a completely regular
network to a completely random graph. Small-world networks not only retain the
high clustering coefficient of regular networks, but also capture the short
average path length of random graphs. Newman and Watts [54] proposed a NW
model which is a modified version of the original WS model. In the NW model,
the shortcuts are added between randomly chosen pairs of nodes, but no
connections are removed from the regular lattice.

The NW model can be described as follows. Firstly, we start with a regular
graph. We consider the nearest-neighbor coupled network containing $N$ nodes.
Each node of the coupled network is connected to its $K$ nearest-neighbors by
undirected edges, where $K$ is an even integer. Secondly, we randomly add some
new connections to the coupled network. With probability $e$, we connect the
pair of nodes chosen randomly.

In this paper, we first generated the three nearest-neighbor coupled networks
containing 6000 nodes, 8000 nodes, and 10000 nodes, respectively. And we set
$K=4$ so that each node of these networks is connected to its 4 nearest-
neighbors. Then we added a connection between pairs of nodes of the three
coupled networks with probability $e=0.0008$, $e=0.0005$, and $e=0.001$,
respectively. For each case, we generated 100 small-world networks using the
NW model. Next, we applied the SB algorithm to calculate their average mass
exponents $\langle\tau(q)\rangle$ and average generalized fractal dimensions
$\langle D(q)\rangle$. The $\langle\tau(q)\rangle$ and $\langle D(q)\rangle$
curves are plotted in Fig. 9. From Fig. 9, we find that the
$\langle\tau(q)\rangle$ and $\langle D(q)\rangle$ curves of small-world
networks are not straight lines, but the fluctuation of these curves is very
small. This means that the multifractality of these small-world networks is
not obvious.

Figure 9: (Color online) Average mass exponents $\langle\tau(q)\rangle$ (upper
panel) and average generalized fractal dimensions $\langle D(q)\rangle$ (lower
panel) for the small-world networks derived from the SB algorithm. Here the
average is calculated from 100 realizations.

###  4.3 Random networks

The Erdős-Rényi (ER) random graph model [59] is one of the classical network
models for generating a completely random network. We start with $N$ isolated
nodes. For every possible pair of nodes, we connect them by an undirected
connection with probability $e$.

In this paper, we considered the three ER random graph models containing 6000
nodes, 8000 nodes, and 10000 nodes. Then we connected all possible node pairs
of the three ER models with probability $e=0.0008$, $e=0.00054$, and
$e=0.0005$, respectively. For each case, we generated 100 random networks by
using the ER model and extract the largest connected component from each
random network. Next, we used the SB algorithm to calculate the average mass
exponents $\langle\tau(q)\rangle$ and average generalized fractal dimensions
$\langle D(q)\rangle$ of these largest connected components. In Fig. 10, we
show the $\langle\tau(q)\rangle$ and $\langle D(q)\rangle$ curves. As we can
see from Fig. 10, the $\langle\tau(q)\rangle$ and $\langle D(q)\rangle$ curves
of random networks are close to straight lines. So this is to say, the
multifractality almost does not exist in these random networks.

Figure 10: (Color online) Average mass exponents $\langle\tau(q)\rangle$
(upper panel) and average generalized fractal dimensions $\langle D(q)\rangle$
(lower panel) for the random networks derived from the SB algorithm. Here the
average is calculated from 100 realizations.

##  5 Conclusion

In this work, we employed the sandbox (SB) algorithm proposed by Tél et al.
[49], for multifractal analysis (MFA) of complex networks.

First we compared the SB algorithm with two existing algorithms of MFA for
complex networks: the compact-box-burning (CBB) algorithm proposed by Furuya
and Yakubo [44], and the improved box-counting (BC) algorithm proposed by Li
et al. [46], by calculating the mass exponents $\tau(q)$ of some deterministic
model networks (networks generated from the $(u,v)$-flower, the minimal model,
and the generalized version of the minimal model). We made a detailed
comparison between the numerical results and the theoretical ones of these
model networks. The comparison results show that the SB algorithm is the most
effective, feasible, and accurate algorithm to calculate the mass exponents
$\tau(q)$ and to explore the multifractal behavior of complex networks.

In addition, we applied the SB algorithm to study the multifractality of some
classic model networks, such as scale-free networks, small-world networks, and
random networks. Our numerical results show that multifractality exists in
scale-free networks, that of small-world networks is not obvious, and it
almost does not exist in random networks.

## Acknowledgments

This project was supported by the Natural Science Foundation of China (Grant
No. 11371016), the Chinese Program for Changjiang Scholars and Innovative
Research Team in University (PCSIRT) (Grant No. IRT1179), the Lotus Scholars
Program of Hunan province of China.

## References

  * [1]
  * [1] C. Song, S. Havlin, and H.A. Makse, Nature(London) 433, 392-395(2005).
  * [2] B.B. Mandelbrot, _The Fractal Geometry of Nature_ (New York: Academic Press, 1983).
  * [3] J. Feder, _Fractals_ (New York: Plenum, 1988).
  * [4] K. Falconer, _Techniques in Fractal Geometry_ (New York: Wiley, 1997).
  * [5] C. Song, L.K. Gallos, S. Havlin, and H.A. Makse, J. Stat. Mech.: Theor. Exp. 3, P03006(2007).
  * [6] J.S. Kim, K.I. Goh, G. Salvi, E. Oh, B. Kahng, and D. Kim, Phys. Rev. E 75, 016110(2007).
  * [7] J.S. Kim, K.I. Goh, B. Kahng, and D. Kim, Chaos 17, 026116(2007).
  * [8] W.X. Zhou, Z.Q. Jiang, and D. Sornette, Physica A 375, 741-752(2007).
  * [9] L. Gao, Y. Hu, and Z. Di, Phys. Rev. E 78, 046109(2008).
  * [10] L. Guo and X. Cai, Chin. Phys. Lett. 26(8), 088901(2009).
  * [11] O. Shanker, Mod.Phys.Lett.B 21, 321-326(2007).
  * [12] G. Bianconi, Chaos 17, 026114(2007).
  * [13] Z.K. Gao and N.D. Jin, Chaos 19, 033137(2009).
  * [14] H.D. Rozenfeld and H.A. Makse, Chem. Eng. Sci. 64, 4572-4575(2009).
  * [15] H.D. Rozenfeld, C. Song, and H. A. Makse, Phys. Rev. Lett. 104, 025701(2010).
  * [16] Z. Zhang, S. Zhou, T. Zou, and G. Chen, J. Stat. Mech.: Theor. Exp. 9, P09008(2008).
  * [17] Z. Zhang, S. Zhou, W. Xie, L. Chen, Y. Lin, and J. Guan, Phys. Rev. E 79, 061113(2009).
  * [18] J.X. Liu and X.M. Kong, Acta Phy. Sin. 59, 2244-2249(2010).
  * [19] Z.Z. Zhang, S.Y. Gao, and W.L. Xie, Chaos 20, 043112(2010).
  * [20] M.F. Dai, X.Y. Li, and L.F. Xi, Chaos 23, 033106(2013).
  * [21] E. Canessa, J. Phys. A: Math. Gen. 33, 3637-3651(2000).
  * [22] V.V. Anh, Q.M. Tieng, and Y.K. Tse, Int. Trans. Oper. Res. 7, 349-363(2000).
  * [23] Z.G. Yu, V. Anh, and K.S. Lau, Physica A 31, 351-361(2001).
  * [24] Z.G. Yu, V. Anh, and K.S. Lau, Phys. Rev. E 64, 031903(2001).
  * [25] Z.G. Yu, V. Anh, and K.S. Lau, Phys. Rev. E 68, 021913(2003).
  * [26] Z.G. Yu, V. Anh, and K.S. Lau, J. Theor. Biol. 226, 341-348(2004).
  * [27] Z.G. Yu, V. Anh, and K.S. Lau, Phys. Rev. E 73 031920(2006).
  * [28] Z.G. Yu, Q.J. Xiao, L. Shi, J.W. Yu, and V. Anh, Chin. Phys. B 19, 068701(2010).
  * [29] V.V. Anh, K.S. Lau, and Z.G. Yu, Phys. Rev. E 66, 031910(2002).
  * [30] L.Q. Zhou, Z.G. Yu, J.Q. Deng, V. Anh, and S.C. Long, J. Theor. Biol. 232 559-567(2005).
  * [31] J.J. Han and W.J. Fu, Chin. Phys. B 19, 010205(2010).
  * [32] S.M. Zhu, Z.G. Yu, and V. Anh, Chin. Phys. B 20, 010505(2011).
  * [33] D. Schertzer and S. Lovejoy, J. Geophys. Res. 92, 9693-9714(1987).
  * [34] S. Lovejoy, M.R. Duncan, and D. Schertzer, J. Geophys. Res. 31D, 26479-426492(1996).
  * [35] S. Lovejoy and D. Schertzer, Computers and Geoscience 36, 1393-1403(2010).
  * [36] S. Lovejoy and D. Schertzer, Computers and Geoscience 36, 1404-1413(2010).
  * [37] J.W. Kantelhardt, E. Koscielny-Bunde, D. Rybski, P. Braun, A. Bunde, and S. Havlin, J. Geophys. Res. 111, D01106(2006).
  * [38] D. Veneziano, A. Langousis, and P. Furcolo, Water Resour. Res. 42, W06D15(2006).
  * [39] V. Venugopal, S.G. Roux, E. Foufoula-Georgiou, and A. Arneodo, Water Resour. Res. 42, W06D14(2006).
  * [40] J.A. Wanliss, V.V. Anh, Z.G. Yu, and S. Watson, J. Geophys. Res. 110, A08214(2005).
  * [41] Z.G. Yu, V. Anh, and R. Eastes, J. Geophys. Res. 114, A05214(2009).
  * [42] Z.G. Yu, V. Anh, Y. Wang, D. Mao, and J. Wanliss, J. Geophys. Res. 115, A10219(2010).
  * [43] C.Y. Lee and S.H. Jung, Phys. Rev. E 73, 066102(2006).
  * [44] S. Furuya and K. Yakubo, Phys. Rev. E 84, 036118(2011).
  * [45] D.L. Wang, Z.G. Yu, and V. Anh, Chin. Phys. B 21(8), 080504(2012).
  * [46] B.G. Li, Z.G. Yu, and Y. Zhou, J. Stat. Mech.: Theor. Exp. 2014, P02020(2014).
  * [47] J.L. Liu, Z.G. Yu, and V. Anh, Phys. Rev. E 89, 032814(2014).
  * [48] L.K. Gallos, C. Song, S. Havlin, and H.A. Makse, Proc. Natl. Acad. Sci. USA. 104(19), 7746-7751(2007).
  * [49] T. Tél, Á. Fülöp, and T. Vicsek, Physica A 159, 155-166(1989).
  * [50] T.C. Halsey, M.H. Jensen, L.P. Kadanoff, I. Procaccia, and B.I. Shraiman, Phys. Rev. A 33, 1141-1151(1986).
  * [51] R.W. Floyd, Commun. ACM 5(6), 345(1962).
  * [52] D.F. Gleich, _A graph library for Matlab based on the boost graph library, http://dgleich.github.com/matlab-bgl._
  * [53] D.J. Watts and S.H. Strogatz, Nature 393(6684), 440-442(1998).
  * [54] M.E.J. Newman and D.J. Watts, Phys. Lett. A 263, 341-346(1999).
  * [55] A.L. Barabási and R. Albert, Science 286(5439), 509-512(1999).
  * [56] S.N. Dorogovtsev, A.V. Goltsev, and J.F.F. Mendes, Phys. Rev. E 65, 066122(2002).
  * [57] H.D. Rozenfeld, S. Havlin, and D. Ben-Avraham, New J. Phys. 9, 175(2007).
  * [58] C. Song, S. Havlin, and H.A. Makse, Nat. Phys. 2, 275-281(2006).
  * [59] P. Erdős and A. Rényi, Publ. Math. Debrecen 6, 290-297(1959).
  * [2]


# Multimode bosonic cat codes with an easily implementable universal gate set

Aurélie Denys<EMAIL_ADDRESS>Anthony Leverrier
<EMAIL_ADDRESS>Inria Paris, France

(August 27, 2024)

###### Abstract

We present a method for designing quantum error correcting codes such that a
specific group of logical operations is implemented using simple physical
operations, provided that this group is a (finite) unitary 1-design. These
physical operations can be transversal gates for qubit codes, or Gaussian
unitaries for bosonic codes. In the latter case, one can exploit this approach
to define multimode extensions of the cat qubit, wherein all single-qubit
Clifford logical gates are obtained from a quadratic Hamiltonian. If a quartic
Hamiltonian is also available, such as a controlled rotation, then it can be
used to implement the CNOT and CS gates, providing a universal gate set.

A main challenge in designing fault-tolerant approaches to quantum computing
is the need to address two seemingly conflicting requirements: protecting
quantum information against various sources of noise, and manipulating the
same quantum information in order to perform a computation. Quantum error
correcting codes offer a solution to the first problem, but good codes tend to
protect information so well that the set of logical gates that can be
performed fault-tolerantly is often very limited. For instance, the Eastin-
Knill theorem puts severe restrictions on the set of transversal gates for a
nontrivial quantum code Eastin and Knill (2009); Faist  _et al._ (2020).

In the case of bosonic codes Terhal  _et al._ (2020); Albert (2022), a number
of single-mode encodings admit interesting logical gate sets that are easily
implementable: Clifford operations for the GKP code Gottesman  _et al._
(2001), Pauli operations and diagonal Clifford gates for cat codes Cochrane
_et al._ (1999); Mirrahimi  _et al._ (2014) and rotation-symmetric codes
Grimsmo  _et al._ (2020). Various additional gadgets can then promote these to
universal gate sets Guillaud and Mirrahimi (2019); Puri  _et al._ (2020);
Grimsmo and Puri (2021). Motivated by a recent construction of a multimode
bosonic code with a nonstandard group of logical gates obtained with Gaussian
unitaries Denys and Leverrier (2023), we ask whether it is possible to
engineer a bosonic code for which the same property holds for a more
interesting logical gate set. We note that a similar idea was investigated by
Gross in the case of encoding a qubit in a spin Gross (2021).

More precisely, in the case of a qubit, and for a subgroup $G$ of $U(2)$, we
wish to find an encoding
$\mathcal{E}:|0\rangle\mapsto|\bar{0}\rangle,|1\rangle\mapsto|\bar{1}\rangle$
such that any logical gate $g\in G$ can be implemented with a Gaussian unitary
$\rho(g)$. A well-known solution to this question is the dual-rail encoding
used for instance in the KLM scheme Knill  _et al._ (2001): the idea is to put
a single photon in two bosonic modes and take the encoding
$|\bar{0}\rangle=|1\rangle|0\rangle$, $|\bar{1}\rangle=|0\rangle|1\rangle$
corresponding to the photon being in the first or second mode. For this
encoding, all single-qubit unitaries are indeed obtained from beamsplitters
and phase-shifts. An issue is that this code is not tolerant to losses since
losing a single photon gives the vacuum state. We are therefore interested in
finding more robust alternatives to the dual-rail encoding. Perhaps
surprisingly, we find that this is indeed possible under relatively mild
assumptions about $G$, essentially that it forms a unitary 1-design, for
instance the Pauli group and the Clifford group. In fact, one can state a
rather general result about the existence of quantum codes allowing to
implement each logical gate $g\in G$ with some physical gate $\rho(g)$, which
can be a Gaussian unitary or a transversal gate.

######  Lemma 1.

Let $\mathcal{H}_{L}:=\mathbbm{C}^{d}$ and $\mathcal{H}_{P}$ be two Hilbert
spaces corresponding to a logical space and a physical space, respectively.
Let $G\leq\mathcal{U}(\mathcal{H}_{L})$ be a subgroup of the unitary group of
$\mathcal{H}_{L}$, and $\rho:G\to\mathcal{U}(\mathcal{H}_{P})$ be a group
homomorphism. Assuming

  1. $(i)$

that there exists a set $\\{g_{i}\\}\subseteq G$ such that
$\\{|i\rangle:=g_{i}|0\rangle\\}$ forms an orthonormal basis of
$\mathcal{H}_{L}$,

  2. $(ii)$

that $G$ is a unitary 1-design,

  3. $(iii)$

that $\Pi_{G}:=\frac{d}{|G|}\sum_{g\in G}\langle
0|g^{\dagger}|0\rangle\rho(g)\neq 0$,

then there exists $|\Phi\rangle\in\mathcal{H}_{P}$ such that the encoding

$\displaystyle\mathcal{E}:\quad\mathcal{H}_{L}$
$\displaystyle\to\mathcal{H}_{P}$ $\displaystyle\quad|i\rangle$
$\displaystyle\mapsto|\bar{i}\rangle:=\mathcal{N}\sum_{g\in G}\langle
0|g^{\dagger}|i\rangle\rho(g)|\Phi\rangle,$ (1)

satisfies

$\displaystyle\mathcal{E}(U|\psi\rangle)=\rho(U)\mathcal{E}(|\psi\rangle),\qquad\forall
U\in G,\quad\forall|\psi\rangle\in\mathcal{H}_{L}.$ (2)

The normalization factor is given by
$\mathcal{N}=\frac{1}{\|\Pi_{G}|\Phi\rangle\|}\frac{d}{|G|}$.

We give the proof in Appendix A. This result is quite general since any
encoding satisfying condition (2) is of the form (1). Indeed, assuming the
existence of a map $\mathcal{E}^{\prime}$ for which (2) holds, define
$|\Phi\rangle=\mathcal{E}^{\prime}(|0\rangle)$, and apply the construction of
the lemma to $|\Phi\rangle$:

$\displaystyle|\bar{i}\rangle$ $\displaystyle=\mathcal{N}\sum_{g\in G}\langle
0|g^{\dagger}|i\rangle\rho(g)\mathcal{E}^{\prime}(|0\rangle)$
$\displaystyle=\mathcal{N}\sum_{g\in G}\mathcal{E}^{\prime}(\langle
0|g^{\dagger}|i\rangle\,g|0\rangle)$
$\displaystyle=\mathcal{N}\frac{|G|}{d}\sum_{k,\ell=0}^{d-1}\mathcal{E}^{\prime}(\langle
0|k\rangle\langle\ell|i\rangle|\ell\rangle\langle k|0\rangle)$
$\displaystyle=\frac{1}{\|\Pi_{G}|\Phi\rangle\|}\mathcal{E}^{\prime}(|i\rangle)=\mathcal{E}^{\prime}(|i\rangle)$

where the second equality exploits the linearity of $\mathcal{E}^{\prime}$ and
(2), and the third equality uses that $G$ is a unitary 1-design. The last
equality results from
$\Pi_{G}\mathcal{E}^{\prime}(|0\rangle)=\mathcal{E}^{\prime}(|0\rangle)$.

Two possible applications are quantum error correcting codes $\llbracket
n,k\rrbracket$ with transversal logical gates, and bosonic codes with Gaussian
logical gates. In the first case, the standard setting is
$\mathcal{H}_{L}=(\mathbbm{C}^{2})^{\otimes k}$ and
$\mathcal{H}_{P}=(\mathbbm{C}^{2})^{\otimes n}$ for $n=mk$ and taking the
group homomorphism $\rho(U)=U^{\otimes m}$, corresponding to a transversal
application of the gate $U$. If one takes $G$ to be the whole unitary group
$\mathcal{U}(\mathcal{H}_{L})$, then the sum in the definition of $\Pi_{G}$
should be replaced by an integral with the Haar measure. The resulting
operator is zero, unless $m=1$, which corresponds to a trivial encoding. This
is consistent with the Eastin-Knill theorem that forbids the existence of a
nontrivial code with a universal gate set. A more interesting group is the
$k$-qubit Clifford group $\mathcal{C}_{k}$, which satisfies the conditions of
the lemma provided $m-1$ is a multiple of 8 (or 4 if the Clifford group does
not contain the 8th root of unity). The idea is that the Clifford group admits
$\langle e^{i\pi/4}\rangle$ as its center, and therefore
$\sum_{g}g^{\dagger}\otimes g^{\otimes n}$ will vanish unless $n-1$ is a
multiple of 8111We thank Markus Heinrich, Richard Küng, and Dominic Wild for
pointing this out.. A word of warning: the codes that we obtain in this way
are not stabilizer codes in general. Similar observations were made in Gross
(2021); Kubischta and Teixeira (2023). As pointed out above, all codes with a
transversal set of gates corresponding to a unitary 1-design can be obtained
in this way. But their error correction capabilities will crucially depend on
the choice of $|\Phi\rangle$, and it is not clear that finding an interesting
initial state is much easier than directly finding a good quantum code with
transversal gates.

The second application, that we explore in more detail now, is concerned with
bosonic codes encoding a qubit. In that case, the set of relevant initial
states is arguably simpler to explore if one focuses for instance on coherent
states. In this setting, we have $\mathcal{H}_{L}=\mathbbm{C}^{2}$ and
$\mathcal{H}_{P}$ is a two-mode Fock space with annihilation operators
$\hat{a},\hat{b}$. The natural homomorphism $\rho$ maps the unitary
$U\in\mathcal{U}(\mathbbm{C}^{2})$ to the Gaussian unitary acting on the
annihilation operators in the following way Weedbrook  _et al._ (2012):
$\rho(U):\left[\begin{smallmatrix}\hat{a}^{\dagger}\\\
\hat{b}^{\dagger}\end{smallmatrix}\right]\mapsto
U\left[\begin{smallmatrix}\hat{a}^{\dagger}\\\
\hat{b}^{\dagger}\end{smallmatrix}\right]$, and mapping a coherent state
$|\vec{\alpha}\rangle$ to $|U\vec{\alpha}\rangle$. A possible choice for the
initial state is an arbitrary coherent state $|\alpha\rangle|\beta\rangle$.
More generally, one could encode a qudit by considering a $d$-mode Fock space,
and the corresponding map $\rho$. Such codes will generalize the one-hot
quantum code Karácsony  _et al._ (2023), which is obtained in the limit of
$\alpha\to 0$.

Three examples of groups that are unitary 1-designs are worth considering.
Taking $U(2)$ (or $SU(2)$), one recovers the encoding of the dual-rail qubit,
which corresponds to no encoding. See Appendix B for details. This is again
expected from the Eastin-Knill theorem. More interestingly, the Pauli and
Clifford groups for one qubit lead to nontrivial codes that we investigate
below. In both cases, the logical states are a superposition of coherent
states lying in a constellation in phase-space. This is reminiscent of quantum
spherical codes Jain  _et al._ (2023). Our codes are not spherical codes in
general since the superpositions are typically far from uniform.

A first example: the Pauli code.— Several versions of the Pauli group exist,
depending on which phases are included: e.g, $\langle X,Z\rangle,\langle
iX,iZ\rangle$ of order 8 and $\langle i,X,Z\rangle$ of order 16. We consider
for instance $\langle X,Z\rangle$. To define the code, we also need to choose
an initial coherent state, and we will see that this choice affects the code
performance as well as the feasibility of the state preparation. For a generic
choice of initial coherent state $|\alpha\rangle|\beta\rangle$, we obtain the
encoding

$\displaystyle|\bar{0}\rangle\propto|c_{1}(\alpha)\rangle|c_{0}(\beta)\rangle,\quad|\bar{1}\rangle\propto|c_{0}(\beta)\rangle|c_{1}(\alpha)\rangle,$
(3)

where we use the notation
$|c_{k}(\alpha)\rangle:=|\alpha\rangle+(-1)^{k}|-\alpha\rangle$ for an
unnormalized two-component single-mode cat state. Note that one recovers the
dual-rail encoding in the limit $\alpha,\beta\to 0$, which suggests that the
code is a finite-energy generalization of the dual-rail qubit. This code
admits a stabilizer $-\rho(-\mathbbm{1})=-e^{i\pi(\hat{n}_{1}+\hat{n}_{2})}$
and jump operators
$\hat{a}^{2}+\hat{b}^{2}-\alpha^{2}-\beta^{2},\hat{a}^{2}\hat{b}^{2}-\alpha^{2}\beta^{2}$,
but is not a stabilizer code unless $\beta=\pm\alpha$. This is because these
operators define a 4-dimensional space
$\mathrm{Span}\\{|\alpha\rangle|\pm\beta\rangle-|-\alpha\rangle|\mp\beta\rangle,|\pm\beta\rangle|\alpha\rangle-|\mp\beta\rangle|-\alpha\rangle\\}$,
and the code is a subspace of dimension 2. By construction, the Pauli
operators can be implemented with Gaussian unitaries. As we will explain
later, one can also obtain a $CZ$ gate between two logical qubits by
performing a controlled rotation ($CROT$) gate, namely
$(-1)^{\hat{n}_{2}\otimes\hat{n}_{4}}$ on two of the four modes, similarly to
what is done for rotation-symmetric bosonic codes Grimsmo  _et al._ (2020).

For this simple code, it is straightforward to simulate the performance for a
pure-loss channel, as was done in Albert  _et al._ (2018) and Denys and
Leverrier (2023): we plot on Fig. 1 the entanglement infidelity for the pure-
loss channel, after the optimal recovery operation. This figure of merit has
the advantage of being efficiently computable provided the constellation size
is not too large, and provides some insight about the protection offered by
the encoding. If the group is $\langle X,Z\rangle$, then the initial state
$|\alpha\rangle|\alpha\rangle$ for $\alpha>0$ has the advantage of yielding a
constellation of minimal size, while $|\alpha\rangle|i\alpha\rangle$ improves
the tolerance to loss. This is consistent with an analysis of the Knill-
Laflamme conditions for the Kraus operators of the pure-loss channel, as
discussed in Appendix C. Taking the variant $\langle i,X,Z\rangle$ of the
Pauli group leads to more complicated states (with a larger constellation) and
appears to degrade the protection against loss.

Figure 1: Entanglement fidelity for the pure-loss channel with loss rate
$\gamma=10^{-2}$ for 3 variants of the Pauli code, depending on the choice of
group $\langle X,Z\rangle$ or $\langle i,X,Z\rangle$ and initial state.

The Clifford code.— Let us now apply the strategy to the single-qubit Clifford
group. Again, there are several variants of this group, and we choose to focus
on $\langle H,S\rangle$, with
$H=\frac{1}{\sqrt{2}}\left[\begin{smallmatrix}1&1\\\
1&-1\end{smallmatrix}\right]$ and $S=\left[\begin{smallmatrix}1&0\\\
0&i\end{smallmatrix}\right]$. This group has order 96. A smaller version of
the group has order 48 (as in Kubischta and Teixeira (2023)), but does not
contain the gate $S$, which will be required later to obtain a universal gate
set. Applying the strategy of Lemma 1 to the initial state
$|\alpha\rangle|\beta\rangle$, we get (see Appendix D for details):

$\displaystyle|\bar{0}\rangle$
$\displaystyle\propto|C_{1}(\alpha)\rangle|C_{0}(\beta)\rangle$
$\displaystyle\quad+\frac{1}{\sqrt{2}}\sum_{\ell=0}^{3}\left|C_{1}\left(\frac{\alpha+i^{\ell}\beta}{\sqrt{2}}\right)\right\rangle\left|C_{0}\left(\frac{\alpha-i^{\ell}\beta}{\sqrt{2}}\right)\right\rangle$
(4)

and $|\bar{1}\rangle$ is obtained by swapping the two modes, that is applying
$\rho(X)$. Here,
$|C_{k}(\gamma)\rangle:=\sum_{\ell=0}^{3}(-i)^{k\ell}|i^{\ell}\gamma\rangle$
is an unnormalized four-component cat state. A possible choice that maximizes
the distance between the states of the constellation is $\beta=\alpha
e^{i\pi/8}$. Again, by construction, all the logical gates from $\langle
H,S\rangle$ are implemented with Gaussian unitaries. To achieve universality,
one requires the ability to perform 2-qubit gates. We do this by applying a
$CROT$ gate $\omega^{\hat{n}_{2}\otimes\hat{n}_{4}}$ on the modes 2 and 4 of a
two-qubit encoding, with $\omega$, an $m$-root of unity. For the Clifford
code, it is sufficient to take $\omega=i$. The following lemma shows that if a
phase gate belongs to $G$, then the corresponding controlled-phase gate can be
implemented with $CROT$.

######  Lemma 2.

If $P_{\omega}=\left[\begin{smallmatrix}1&0\\\
0&\omega\end{smallmatrix}\right]\in G\leq\mathcal{U}(\mathbbm{C}^{2})$, then
$\omega^{\hat{n}_{2}\otimes\hat{n}_{4}}$ implements a $CP_{\omega}$ gate on
the two qubits. In particular, for the Clifford group $\langle H,S\rangle$,
the controlled-rotations $(-1)^{\hat{n}_{2}\otimes\hat{n}_{4}}$ and
$i^{\hat{n}_{2}\otimes\hat{n}_{4}}$ implement a $CZ$ and a $CS$ gate on the
two logical qubits, respectively.

For the Clifford code, we see that $CZ$ gives the whole multiqubit Clifford
group, whereas the $CS$ gate promotes the gate set to a universal one Amy  _et
al._ (2020); Bian and Selinger (2023).

###### Proof.

Let $\omega=e^{i2\pi/m}$. For two arbitrary coherent states, by expanding the
states in the Fock basis, it is straightforward to check that the $CROT$ gate
gives

$\omega^{\hat{n}_{1}\otimes\hat{n}_{2}}|\alpha\rangle|\beta\rangle=\frac{1}{m}\sum_{p,q=0}^{m-1}\omega^{-pq}|\omega^{p}\alpha\rangle|\omega^{q}\beta\rangle.$

If the gate $P_{\omega}$ belongs to the group, then
$\omega^{\hat{n}_{2}\otimes\hat{n}_{4}}|\alpha,\beta\rangle|\gamma,\delta\rangle$
is given by
$\frac{1}{m}\sum_{p,q=0}^{m-1}\omega^{-pq}|\alpha,\omega^{p}\beta\rangle|\gamma,\omega^{q}\delta\rangle$,
that is

$\displaystyle\frac{1}{m}\sum_{p,q=0}^{m-1}\omega^{-pq}\rho(P_{\omega}^{p}\otimes
P_{\omega}^{q})|\alpha,\beta\rangle|\gamma,\delta\rangle.$

In other words, $\omega^{\hat{n}_{2}\otimes\hat{n}_{4}}$ implements a two-
qubit logical gate given by

$\frac{1}{m}\sum_{p,q=0}^{m-1}\omega^{-pq}\begin{bmatrix}1&0\\\
0&\omega^{p}\end{bmatrix}\otimes\begin{bmatrix}1&0\\\
0&\omega^{q}\end{bmatrix}=\mathrm{diag}(1,1,1,\omega)=CP_{\omega}.\qed$

The Clifford code is not a stabilizer code, but it admits some stabilizers:
$i^{\hat{n}_{1}\otimes\hat{n}_{2}}$, or equivalently
$-i\rho(i\mathbbm{1}_{2})$.

State preparation and measurements.— We have already pointed out that our
codes are not stabilizer codes in general, which suggests that state
preparation may be complex. For the Pauli code of (3), we observe that the two
basis states are product states with a state of a two-component cat code in
each mode. Given that such states are routinely prepared and manipulated in
the lab today, we expect the state preparation to be feasible in the near
term. On the other hand, the preparation of the state in (Multimode bosonic
cat codes with an easily implementable universal gate set) for the Clifford
code seems significantly more delicate.

Measuring the states in the (logical) computational basis can be done naively
simply by counting the photons in each of the modes. This is because the
logical states are superpositions of Fock states of the form
$|pm+1\rangle|pn\rangle$, where $p=2$ for the Pauli code and $p=4$ or 8 for
the Clifford code. There again, techniques developed for cat codes may prove
useful since the objective is similar, namely distinguishing logical cat
states.

Discussion and open questions.— We have introduced a general methodology to
design quantum error correcting codes that admit a specific logical group
implementable with simple physical gates, either transversal gates for qubit
codes or Gaussian unitaries for bosonic codes. In the latter case, one can
design a code with a universal gate set consisting of such Gaussian unitaries
together with some controlled-rotation gates, which correspond to quartic
Hamiltonians. While these gates are certainly more challenging to implement,
they do not seem out of reach for circuit QED Blais  _et al._ (2021). On the
other hand, such gates are much more difficult for photonic implementations,
and it would be interesting to understand whether other gadgets can be used to
obtain a universal gate set in that case.

A question that we have not addressed at all is how to perform error
correction. More generally, the question of how to stabilize such states
appears quite daunting since they involve rather large constellations of
coherent states.

###### Acknowledgements.

We thank Mazyar Mirrahimi and Christophe Vuillot for many discussions about
bosonic codes and for feedback on a preliminary version of this manuscript. We
acknowledge support from the Plan France 2030 through the project NISQ2LSQ
ANR-22-PETQ-0006.

## References

  * Eastin and Knill (2009) Bryan Eastin and Emanuel Knill, “Restrictions on transversal encoded quantum gate sets,” Phys. Rev. Lett. 102, 110502 (2009).
  * Faist  _et al._ (2020) Philippe Faist, Sepehr Nezami, Victor V. Albert, Grant Salton, Fernando Pastawski, Patrick Hayden, and John Preskill, “Continuous symmetries and approximate quantum error correction,” Phys. Rev. X 10, 041018 (2020).
  * Terhal  _et al._ (2020) B M Terhal, J Conrad, and C Vuillot, “Towards scalable bosonic quantum error correction,” Quantum Science and Technology 5, 043001 (2020).
  * Albert (2022) Victor V Albert, “Bosonic coding: introduction and use cases,” arXiv preprint arXiv:2211.05714 (2022).
  * Gottesman  _et al._ (2001) Daniel Gottesman, Alexei Kitaev, and John Preskill, “Encoding a qubit in an oscillator,” Phys. Rev. A 64, 012310 (2001).
  * Cochrane  _et al._ (1999) P. T. Cochrane, G. J. Milburn, and W. J. Munro, “Macroscopically distinct quantum-superposition states as a bosonic code for amplitude damping,” Phys. Rev. A 59, 2631–2634 (1999).
  * Mirrahimi  _et al._ (2014) Mazyar Mirrahimi, Zaki Leghtas, Victor V Albert, Steven Touzard, Robert J Schoelkopf, Liang Jiang, and Michel H Devoret, “Dynamically protected cat-qubits: a new paradigm for universal quantum computation,” New Journal of Physics 16, 045014 (2014).
  * Grimsmo  _et al._ (2020) Arne L. Grimsmo, Joshua Combes, and Ben Q. Baragiola, “Quantum computing with rotation-symmetric bosonic codes,” Phys. Rev. X 10, 011058 (2020).
  * Guillaud and Mirrahimi (2019) Jérémie Guillaud and Mazyar Mirrahimi, “Repetition cat qubits for fault-tolerant quantum computation,” Phys. Rev. X 9, 041053 (2019).
  * Puri  _et al._ (2020) Shruti Puri, Lucas St-Jean, Jonathan A. Gross, Alexander Grimm, Nicholas E. Frattini, Pavithran S. Iyer, Anirudh Krishna, Steven Touzard, Liang Jiang, Alexandre Blais, Steven T. Flammia, and S. M. Girvin, “Bias-preserving gates with stabilized cat qubits,” Science Advances 6, eaay5901 (2020).
  * Grimsmo and Puri (2021) Arne L. Grimsmo and Shruti Puri, “Quantum Error Correction with the Gottesman-Kitaev-Preskill Code,” PRX Quantum 2, 020101 (2021).
  * Denys and Leverrier (2023) Aurélie Denys and Anthony Leverrier, “The $2T$-qutrit, a two-mode bosonic qutrit,” Quantum 7, 1032 (2023).
  * Gross (2021) Jonathan A. Gross, “Designing codes around interactions: The case of a spin,” Phys. Rev. Lett. 127, 010504 (2021).
  * Knill  _et al._ (2001) Emanuel Knill, Raymond Laflamme, and Gerald J Milburn, “A scheme for efficient quantum computation with linear optics,” Nature 409, 46–52 (2001).
  * Kubischta and Teixeira (2023) Eric Kubischta and Ian Teixeira, “A nonadditive quantum code with exotic transversal gate set,” arXiv preprint arXiv:2305.07023 (2023).
  * Weedbrook  _et al._ (2012) Christian Weedbrook, Stefano Pirandola, Raúl García-Patrón, Nicolas J. Cerf, Timothy C. Ralph, Jeffrey H. Shapiro, and Seth Lloyd, “Gaussian quantum information,” Rev. Mod. Phys. 84, 621–669 (2012).
  * Karácsony  _et al._ (2023) Márton Karácsony, László Oroszlány, and Zoltán Zimborás, “Efficient qudit based scheme for photonic quantum computing,” arXiv preprint arXiv:2302.07357 (2023).
  * Jain  _et al._ (2023) Shubham P Jain, Joseph T Iosue, Alexander Barg, and Victor V Albert, “Quantum spherical codes,” arXiv preprint arXiv:2302.11593 (2023).
  * Albert  _et al._ (2018) Victor V. Albert, Kyungjoo Noh, Kasper Duivenvoorden, Dylan J. Young, R. T. Brierley, Philip Reinhold, Christophe Vuillot, Linshu Li, Chao Shen, S. M. Girvin, Barbara M. Terhal, and Liang Jiang, “Performance and structure of single-mode bosonic codes,” Phys. Rev. A 97, 032346 (2018).
  * Amy  _et al._ (2020) Matthew Amy, Andrew N. Glaudell, and Neil J. Ross, “Number-Theoretic Characterizations of Some Restricted Clifford+T Circuits,” Quantum 4, 252 (2020).
  * Bian and Selinger (2023) Xiaoning Bian and Peter Selinger, “Generators and relations for 3-qubit Clifford$+$CS operators,” arXiv preprint arXiv:2306.08530 (2023).
  * Blais  _et al._ (2021) Alexandre Blais, Arne L. Grimsmo, S. M. Girvin, and Andreas Wallraff, “Circuit quantum electrodynamics,” Rev. Mod. Phys. 93, 025005 (2021).

## Additional material

##  Appendix A Proof of Lemma 1

###### Proof.

We first check that the operator $\Pi_{G}$ is a projector:

$\displaystyle\Pi_{G}^{2}$ $\displaystyle=\frac{d^{2}}{|G|^{2}}\sum_{g,h\in
G}\langle 0|g^{\dagger}|0\rangle\langle 0|h^{\dagger}|0\rangle\rho(g)\rho(h)$
$\displaystyle=\frac{d^{2}}{|G|^{2}}\sum_{g,h\in G}\langle
0|g^{\dagger}|0\rangle\langle 0|h^{\dagger}|0\rangle\rho(gh)$
$\displaystyle=\frac{d^{2}}{|G|^{2}}\sum_{g,h_{1}\in G}\langle
0|g^{\dagger}|0\rangle\langle 0|h_{1}^{\dagger}g|0\rangle\rho(h_{1})$
$\displaystyle=\frac{d}{|G|}\sum_{g\in G}\sum_{i,j=0}^{d-1}\langle
0|i\rangle\langle j|0\rangle\langle 0|h_{1}^{\dagger}|j\rangle\langle
i|0\rangle\rho(h_{1})$ $\displaystyle=\Pi_{G},$

where the second inequality holds because $\rho$ is a group homorphism, the
third equality follows from the change of variable $h_{1}=gh$, and the fourth
equality relies on the property of a unitary 1-design $\frac{1}{|G|}\sum_{g\in
G}g\otimes g^{\dagger}=\frac{1}{d}\sum_{i,j=0}^{d-1}|i\rangle\langle
j|\otimes|j\rangle\langle i|$.

Provided that this projector is nonzero, one can pick a state
$|\Phi\rangle\in\mathcal{H}_{P}$ such that $\|\Pi_{G}|\Phi\rangle\|\neq 0$ and
observe that
$|\bar{0}\rangle=\frac{1}{\|\Pi_{G}|\Phi\rangle\|}\Pi_{G}|\Phi\rangle$. Let us
compute the overlap $\langle\bar{0}|\rho(g)|\bar{0}\rangle$ for an arbitrary
group element $g\in G$ and show that it equals $\langle 0|g|0\rangle$,

$\displaystyle\langle\bar{0}|\rho(g)|\bar{0}\rangle$
$\displaystyle=\frac{1}{\|\Pi_{G}|\Phi\rangle\|^{2}}\,\langle\Phi|\Pi_{G}\rho(g)\Pi_{G}|\Phi\rangle$
$\displaystyle=\frac{1}{\|\Pi_{G}|\Phi\rangle\|^{2}}\frac{d^{2}}{|G|^{2}}\sum_{h_{1},h_{2}\in
G}\langle 0|h_{1}|0\rangle\langle
0|h_{2}^{\dagger}|0\rangle\langle\Phi|\rho(h_{1}^{\dagger})\rho(g)\rho(h_{2})|\Phi\rangle$
$\displaystyle=\frac{1}{\|\Pi_{G}|\Phi\rangle\|^{2}}\frac{d^{2}}{|G|^{2}}\sum_{h_{1},h_{2}\in
G}\langle 0|h_{1}|0\rangle\langle
0|h_{2}^{\dagger}|0\rangle\langle\Phi|\rho(h_{1}^{\dagger}gh_{2})|\Phi\rangle$
(homomorphism)
$\displaystyle=\frac{1}{\|\Pi_{G}|\Phi\rangle\|^{2}}\frac{d^{2}}{|G|^{2}}\sum_{h_{1},h\in
G}\langle 0|h_{1}|0\rangle\langle
0|h^{\dagger}h_{1}^{\dagger}g|0\rangle\langle\Phi|\rho(h)|\Phi\rangle$
$\displaystyle[h=h_{1}^{\dagger}gh_{2}]$
$\displaystyle=\frac{1}{\|\Pi_{G}|\Phi\rangle\|^{2}}\frac{d}{|G|}\sum_{h\in
G}\langle\Phi|\rho(h)|\Phi\rangle\sum_{i,j=0}^{d-1}\langle 0|i\rangle\langle
j|0\rangle\langle 0|h^{\dagger}|j\rangle\langle i|g|0\rangle$ (unitary
1-design)
$\displaystyle=\frac{1}{\|\Pi_{G}|\Phi\rangle\|^{2}}\frac{d}{|G|}\sum_{h\in
G}\langle\Phi|\rho(h)|\Phi\rangle\sum_{i,j=0}^{d-1}\langle 0|i\rangle\langle
i|g|0\rangle\langle 0|h^{\dagger}|j\rangle\langle j|0\rangle$
$\displaystyle=\frac{1}{\|\Pi_{G}|\Phi\rangle\|^{2}}\frac{d}{|G|}\sum_{h\in
G}\langle\Phi|\rho(h)|\Phi\rangle\langle 0|h^{\dagger}|0\rangle\langle
0|g|0\rangle$
$\displaystyle=\frac{1}{\|\Pi_{G}|\Phi\rangle\|^{2}}\langle\Phi|\Pi_{G}|\Phi\rangle\langle
0|g|0\rangle$ $\displaystyle=\langle 0|g|0\rangle.$

From the choice of $g_{i}$ such that $|i\rangle=g_{i}|0\rangle$, we get

$\displaystyle\rho(g_{i})|\bar{0}\rangle$ $\displaystyle=\mathcal{N}\sum_{g\in
G}\langle 0|g^{\dagger}|0\rangle\rho(g_{i}g)|\Phi\rangle$
$\displaystyle=\mathcal{N}\sum_{g^{\prime}\in G}\langle
0|g^{\prime{\dagger}}g_{i}|0\rangle\rho(g^{\prime})|\Phi\rangle$
$\displaystyle[g^{\prime}=g_{i}g]$
$\displaystyle=\mathcal{N}\sum_{g^{\prime}\in G}\langle
0|g^{\prime{\dagger}}|i\rangle\rho(g^{\prime})|\Phi\rangle=|\bar{i}\rangle.$

This implies
$\langle\bar{i}|\rho(g)|\bar{j}\rangle=\langle\bar{0}|\rho(g_{i}^{\dagger}gg_{j})|\bar{0}\rangle=\langle
0|g_{i}^{\dagger}gg_{j}|0\rangle=\langle i|g|j\rangle$. Taking
$g=\mathbbm{1}$, we obtain that $\\{|\bar{i}\rangle\\}$ is an orthonormal
basis of $\mathcal{E}(\mathcal{H}_{L})$. Finally, for any
$|\psi\rangle\in\mathcal{H}_{L}$, and any $g\in G$, we have

$\displaystyle\rho(g)\mathcal{E}(|\psi\rangle)$
$\displaystyle=\sum_{i,j=0}^{d-1}|\bar{i}\rangle\langle\bar{i}|\rho(g)\mathcal{E}(|j\rangle\langle
j|\psi\rangle)$ $\displaystyle=\sum_{i,j=0}^{d-1}\langle
j|\psi\rangle\langle\bar{i}|\rho(g)|\bar{j}\rangle|\bar{i}\rangle$
$\displaystyle=\sum_{i,j=0}^{d-1}\langle j|\psi\rangle\langle
i|g|j\rangle|\bar{i}\rangle$
$\displaystyle=\sum_{i=0}^{d-1}|\psi\rangle\langle
i|g|\psi\rangle|\bar{i}\rangle$ $\displaystyle=\mathcal{E}(g|\psi\rangle),$

which proves the claim. ∎

##  Appendix B Recovering the dual-rail qubit when $G=SU(2)$

Our goal is to show that

$\Pi_{SU(2)}=2\int_{g\in SU(2)}\langle 0|g^{\dagger}|0\rangle\rho(g)dg$

only has support on the subspace with a single photon in two modes. Here the
integration is over the Haar measure. Let us write
$g=\left[\begin{smallmatrix}g_{00}&g_{01}\\\
g_{10}&g_{11}\end{smallmatrix}\right]$ and compute
$\Pi_{SU(2)}|m\rangle|n\rangle$, where $|m\rangle$ and $|n\rangle$ are Fock
states with $m$ and $n$ photons respectively:

$\displaystyle\Pi_{SU(2)}|m\rangle|n\rangle$ $\displaystyle=2\int_{g\in
SU(2)}dg\;\langle 0|g^{\dagger}|0\rangle\rho(g)|m\rangle|n\rangle$
$\displaystyle=\frac{1}{\sqrt{m!n!}}\int_{g\in SU(2)}dg\;\langle
0|g^{\dagger}|0\rangle\rho(g)\hat{a}^{{\dagger}m}\hat{b}^{{\dagger}n}|0\rangle|0\rangle$
$\displaystyle=\frac{1}{\sqrt{m!n!}}\int_{g\in
SU(2)}dg\;g_{00}^{*}(g_{00}\hat{a}^{\dagger}+g_{01}\hat{b}^{\dagger})^{m}(g_{10}\hat{a}^{\dagger}+g_{11}\hat{b}^{\dagger})^{n}|0\rangle|0\rangle$

The integrand is a homogeneous polynomial of degree $m+n$ in the entries of
$g$ and degree 1 in their complex conjugates. It is well known that this
integral vanishes unless $m+n=1$, which corresponding to the subspace
$\mathrm{Span}(|1\rangle|0\rangle,|0\rangle|1\rangle)$ of the dual-rail qubit
encoding.

##  Appendix C Knill-Laflamme conditions for the pure-loss channel

In this section, we analyse the effect of the pure-loss channel for the Pauli
code of (3), but similar conclusions hold for the Clifford code.

The Kraus operators of the two-mode pure-loss channel $\mathcal{L}_{\gamma}$
with loss rate $\gamma\in[0,1)$ are given by Albert  _et al._ (2018):

$\displaystyle
E_{p_{1},p_{2}}=\left(\frac{\gamma}{1-\gamma}\right)^{(p_{1}+p_{2})/2}\frac{\hat{a}^{p_{1}}\hat{b}^{p_{2}}}{\sqrt{p_{1}!p_{2}!}}(1-\gamma)^{(\hat{n}_{1}+\hat{n}_{2})/2},$

so that the action of the channel $\mathcal{L}_{\gamma}$ is

$\mathcal{L}_{\gamma}(\rho)=\sum_{p_{1},p_{2}=0}^{\infty}E_{p_{1},p_{2}}\rho
E_{p_{1},p_{2}}^{\dagger}.$

A straightforward calculation shows that these Kraus operators attenuate
coherent states: defining $\mu:=\sqrt{1-\gamma}$, we get

$\hat{a}^{p}(1-\gamma)^{\hat{n}/2}|\alpha\rangle=\left(\mu\alpha\right)^{p}e^{-\gamma|\alpha|^{2}/2}|\mu\alpha\rangle.$

We can apply $E_{p_{1},p_{2}}$ to a product of two cat states
$|c_{j}(\alpha_{j})\rangle|c_{1-j}(\alpha_{1-j})\rangle$ for $j\in\\{0,1\\}$,
it gives

$\displaystyle
E_{p_{1},p_{2}}|c_{j}(\alpha_{j})\rangle|c_{j-1}(\alpha_{j-1})\rangle$
$\displaystyle=\left(\frac{\gamma}{1-\gamma}\right)^{(p_{1}+p_{2})/2}\frac{(\mu\alpha_{j})^{p_{1}}(\mu\alpha_{1-j})^{p_{2}}}{\sqrt{p_{1}!p_{2}!}}e^{-\gamma(|\alpha_{j}|^{2}+|\alpha_{1-j}|^{2})/2}|c_{j-p_{1}}(\mu\alpha_{j})\rangle|c_{1-j-p_{2}}(\mu\alpha_{1-j})\rangle$
$\displaystyle=e^{-\gamma(|\alpha_{j}|^{2}+|\alpha_{1-j}|^{2})/2}\frac{\gamma^{(p_{1}+p_{2})/2}\alpha_{j}^{p_{1}}\alpha_{1-j}^{p_{2}}}{\sqrt{p_{1}!p_{2}!}}|c_{j-p_{1}}(\mu\alpha_{j})\rangle|c_{1-j-p_{2}}(\mu\alpha_{1-j})\rangle$
$\displaystyle=f(\alpha_{j},\alpha_{1-j},p_{1},p_{2})|c_{j-p_{1}}(\mu\alpha_{j})\rangle|c_{1-j-p_{2}}(\mu\alpha_{1-j})\rangle$

with indices taken modulo 2.

Focusing on the Pauli code with initial state $|\alpha\rangle|\alpha
e^{i\theta}\rangle$, one can check the Knill-Laflamme conditions for the Kraus
operators of the pure-loss channel by applying the previous expression for
$\alpha_{0}:=\alpha e^{i\theta},\alpha_{1}:=\alpha$,

$\displaystyle\langle\bar{k}|E_{p_{1},p_{2}}^{\dagger}E_{q_{1},q_{2}}|\bar{\ell}\rangle$
$\displaystyle\propto
f(\alpha_{k},\alpha_{1-k},p_{1},p_{2})^{*}f(\alpha_{\ell},\alpha_{1-\ell},q_{1},q_{2})\langle
c_{k-p_{1}}(\mu\alpha_{k})|c_{\ell-q_{1}}(\mu\alpha_{\ell})\rangle\langle
c_{1-k-p_{2}}(\mu\alpha_{1-k})|c_{1-\ell-q_{2}}(\mu\alpha_{1-\ell})\rangle.$

Diagonal terms of the form
$\langle\bar{k}|E_{p_{1},p_{2}}^{\dagger}E_{q_{1},q_{2}}|\bar{k}\rangle$ are
nonzero only if $p_{1}=q_{1}$ and $p_{2}=q_{2}$ since the even and odd cat
states are orthogonal (they have support on the even and odd Fock states,
respectively). In that case, we have

$\displaystyle\langle\bar{k}|E_{p_{1},p_{2}}^{\dagger}E_{p_{1},p_{2}}|\bar{k}\rangle$
$\displaystyle\propto
f(\alpha_{k},\alpha_{1-k},p_{1},p_{2})^{*}f(\alpha_{k},\alpha_{1-k},p_{1},p_{2})$
$\displaystyle=e^{-\gamma(|\alpha_{k}|^{2}+|\alpha_{1-k}|^{2})}\frac{\gamma^{(p_{1}+p_{2})}|\alpha_{k}|^{2p_{1}}|\alpha_{1-k}|^{2p_{2}}}{{p_{1}!p_{2}!}}$

which is independent of $k$, provided that $|\alpha_{0}|=|\alpha_{1}|$.

Moreover, for nondiagonal terms, we observe that if $\theta$ is not a multiple
of $\pi$, then the overlaps between the cat states of amplitude
$\mu\alpha_{0}$ and $\mu\alpha_{1}$ always vanish in the limit of large
energy, $\alpha\to\infty$,

$\lim_{\alpha\to\infty}\langle\bar{0}|E_{p_{1},p_{2}}^{\dagger}E_{q_{1},q_{2}}|\bar{1}\rangle=0.$

This suggests to optimize the choice of the phase $\theta$ so as to maximize
the distance between the constellations of coherent states for the encoded
states $|\bar{0}\rangle$ and $|\bar{1}\rangle$. A similar observation was made
in Jain  _et al._ (2023) where this distance was computed explicitly for many
families of quantum spherical codes. In the case of the Pauli code of (3), we
see that the choices $\theta=0$ and $\theta=\pi/2$ are respectively the worst
and best choices with that respect. This is confirmed numerically, and can
also be seen on Fig. 1.

In the case of the Clifford code of (Multimode bosonic cat codes with an
easily implementable universal gate set), the cat states have either 4 (or 8)
components (depending on whether $e^{i\pi/4}$ is included or not in the
Clifford group), and the Knill-Laflamme conditions are again satisfied in the
limit of large $\alpha$ provided that $\theta$ is not a multiple of $\pi/2$
(or $\pi/4$).

##  Appendix D Expression of the Clifford code for an arbitrary initial state

We consider the group $\langle H,S\rangle$ of order 96 with elements

$\left\\{\begin{bmatrix}i^{k}&0\\\
0&i^{\ell}\end{bmatrix},\quad\begin{bmatrix}0&i^{k}\\\
i^{\ell}&0\end{bmatrix},\quad\frac{1}{\sqrt{2}}\begin{bmatrix}i^{k}&i^{k+\ell}\\\
i^{k+m}&-i^{k+\ell+m}\end{bmatrix},\quad k,\ell,m\in\\{0,1,2,3\\}\right\\}$

and apply the construction of Lemma 1 for an arbitrary initial coherent state
$|\alpha\rangle|\beta\rangle$:

$\displaystyle|\bar{0}\rangle$
$\displaystyle\propto\sum_{k,\ell=0}^{3}i^{-k}|i^{k}\alpha\rangle|i^{\ell}\beta\rangle+\frac{1}{\sqrt{2}}\sum_{k,\ell,m=0}^{3}i^{-k}\left|\frac{i^{k}}{\sqrt{2}}(\alpha+i^{\ell}\beta)\right\rangle\left|\frac{i^{k+m}}{\sqrt{2}}(\alpha-i^{\ell}\beta)\right\rangle$
$\displaystyle\propto|C_{1}(\alpha)\rangle|C_{0}(\beta)\rangle+\frac{1}{\sqrt{2}}\sum_{\ell=0}^{3}\left|C_{1}\left(\frac{\alpha+i^{\ell}\beta}{\sqrt{2}}\right)\right\rangle\left|C_{0}\left(\frac{\alpha-i^{\ell}\beta}{\sqrt{2}}\right)\right\rangle.$

State $|\bar{1}\rangle$ is obtained by applying $\rho(X)$ to
$|\bar{0}\rangle$, which corresponds to swapping the two modes:

$\displaystyle|\bar{1}\rangle$
$\displaystyle\propto|C_{0}(\beta)\rangle|C_{1}(\alpha)\rangle+\frac{1}{\sqrt{2}}\sum_{\ell=0}^{3}\left|C_{0}\left(\frac{\alpha-i^{\ell}\beta}{\sqrt{2}}\right)\right\rangle\left|C_{1}\left(\frac{\alpha+i^{\ell}\beta}{\sqrt{2}}\right)\right\rangle.$