# A Unified Substrate for Body-Brain Co-evolution

Sidney Pontes-Filho1,2,∗, Kathryn Walker3, Elias Najarro3, Stefano
Nichele1,4,5 and Sebastian Risi3

1Department of Computer Science, Oslo Metropolitan University, Oslo (Norway)
2Department of Computer Science, Norwegian University of Science and
Technology, Trondheim (Norway)
3Digital Design Department, IT University of Copenhagen, Copenhagen (Denmark)
4Department of Holistic Systems, Simula Metropolitan Centre for Digital
Engineering, Oslo (Norway)
5Department of Computer Science and Communication, Østfold University College,
Halden (Norway)
∗Corresponding author<EMAIL_ADDRESS>
###### Abstract

The discovery of complex multicellular organism development took millions of
years of evolution. The genome of such a multicellular organism guides the
development of its body from a single cell, including its control system. Our
goal is to imitate this natural process using a single neural cellular
automaton (NCA) as a genome for modular robotic agents. In the introduced
approach, called _Neural Cellular Robot Substrate_ (NCRS), a single NCA guides
the growth of a robot and the cellular activity which controls the robot
during deployment. We also introduce three benchmark environments, which test
the ability of the approach to grow different robot morphologies. In this
paper, NCRSs are trained with covariance matrix adaptation evolution strategy
(CMA-ES), and covariance matrix adaptation MAP-Elites (CMA-ME) for quality
diversity, which we show leads to more diverse robot morphologies with higher
fitness scores. While the NCRS can solve the easier tasks from our benchmark
environments, the success rate reduces when the difficulty of the task
increases. We discuss directions for future work that may facilitate the use
of the NCRS approach for more complex domains.

##  1 Introduction

Multicellular organisms are made of cells that can divide into many, which
specialize in controlling and maintaining the body, sensing the environment,
or protecting from external threats. Such features were acquired by evolution
from the first living cell. After millions of years, colonies of unicellular
organisms appeared and were essential to the development of multicellular
organisms with cellular differentiation [26]. Developmental biologists study
that the growth and specialization of an organism are coordinated by its
genetic code [30].

The field of artificial life tries to create life-like computational models
taking ideas from biological life, such as decentralized and local control
[18]. One of the sub-fields of artificial life, artificial development [13,
7], focuses on modeling or simulating cell division and differentiation. The
techniques applied in artificial development are often based on the indirect
encoding of developmental rules (i.e. analogous to the genome of a biological
organism describing its phenotype). This type of encoding facilitates the
scaling of an organism because the information in the genome is much smaller
than in the resulting phenotype. This property is referred to as genomic
bottleneck [41, 36], and it implies that the genetic code of an organism
compresses the information to grow and maintain its body, and in some species
even complex brains.

One of the simplest computational models of artificial life or dynamical
systems is a cellular automaton (CA) [38]. A CA can be described as a universe
with discrete space and time, which is governed by local rules without any
central control. Such a discrete space is divided into a regular grid of cells
and can possess any number of dimensions. The most commonly studied CAs have
one or two dimensions and their most well-known versions are, respectively,
elementary CA [38] and Conway’s Game of Life [6]. Both have cells with binary
states, but other CA can have many discrete states or continuous ones. In the
1940s, the first CA was introduced by Ulam and von Neumann [35]. Von Neumann
aimed to produce self-replicating machines, and Ulam worked on crystal growth.
In 2002, a CA with rules defined by an artificial neural network was described
[19]. Nowadays, this type of approach is called neural cellular automaton
(NCA). In 2017, Nichele et al. [25] presented an NCA that has developmental
features that were learned through neuroevolution using a method called
compositional pattern-producing network [31]. Recently, Mordvintsev et al.
[23] introduced a differentiable NCA, which possesses growth and regeneration
properties. In their work, an NCA is trained through gradient descent to grow
a colored image from one active "seed" cell.

In evolutionary robotics, co-evolution of morphology and control has the
inherent challenge of optimizing two different features in parallel [2]. It
also presents scalability issues when it deals with modular robots [40]. Our
goal is to implement an approach where the optimization happens in just one
dynamical substrate with local interactions. Here we introduce such a system,
a _Neural Cellular Robot Substrate_ (NCRS), in which a single NCA grows the
morphology of an agent’s body and also controls how that agent interacts with
its environment. The NCA has two phases (Fig. 1). First is the developmental
phase, in which the robot’s body is grown, including where to place its
sensors and actuators. In the following control phase, copies of the same NCA
are running in each cell of the agent, taking into account only local
information from neighboring cells to determine their next state. The
optimization task thus entails figuring out how to transmit information from
the robot’s sensors to its actuators to perform the task at hand.

We also introduce a virtual environment with three benchmark tasks for
evaluating the NCRS’ capacity of designing a robot and then controlling it.
Two benchmarks consist in growing and controlling a robot to approach a light
source (Fig. 1b and Fig. 2(a)). The third task challenges the robot to carry a
ball to a target area. In this benchmark, a second type of rudimentary eye is
added, so the robot can differentiate the ball and the target area (Fig.
2(b)).

The main contribution of this work is the introduction of a single neural
cellular automaton that first grows an agent’s body and then controls it
during deployment. While the solved benchmark domains are relatively simple,
the unified substrate for both body and brain opens up interesting future
research directions, such as opportunities for open-ended evolution [32]. The
source code of this project and videos of the results are available at
https://github.com/sidneyp/neural-cellular-robot-substrate.

Figure 1: Neural Cellular Robot Substrate (NCRS) In the developmental phase
(a), the robot is grown from an initial starting seed, guided by a neural
cellular automaton (c). Once grown, the same neural cellular automaton
determines how the signals propagate through the robot’s morphology during the
control phase (b).

(a) Light chasing with obstacle task (b) Carrying ball to target task

Figure 2: Extensions from the light chasing task. (2(a)) It depicts the
original size of the playing field, which is 60.

##  2 Related work

The co-design of robot bodies and brains has been an active area of research
for decades [22, 29, 16, 37, 10]. Brain and body co-design stands for
producing a control policy and a morphology for a robotic system. For example,
in the work of Lipson and Pollack [20] the same genome directly encodes the
robot’s body and the artificial neural network for control. A method that uses
genetic regulatory networks to define separately a body and an artificial
neural network was introduced by Bongard and Pfeifer [3] and named artificial
ontogeny. The evolved robots are able to locomote and push blocks in noisy
environments. More recent work by Bhatia et al. [2] presents several virtual
environments and also an algorithm for brain and body co-design with separated
description methods for the morphology and control. In comparison with NCRS,
our co-design algorithm consists of only one neural cellular automaton.

The work on NCAs by Mordvintsev et al. [23] is one of the first examples of
self-organizing and self-repair systems that use differentiable models as
rules for cellular automata. Before that, NCA models were typically optimized
with genetic algorithms [25]. After the work on growing NCA, other neural CAs
were introduced, including methods optimized without differentiable
programming. There exist other generative methods for growing 3D artifacts and
functional machines [33], for developing soft robots [15]. Moreover, an NCA
was used as a decentralized classifier of handwritten digits [28].

The developmental phase of our approach is similar to the generative method
with NCA for level design trained with CMA-ME in the work of Earle et al. [8].
Morphology design is also present in other works [14, 34, 17, 5]. The control
phase is based on the NCA for controlling a cart-pole agent introduced by
Variengien et al. [36], but their NCA is trained using a reinforcement
learning algorithm named deep-Q learning and the communication between NCA and
environment happens in predefined cells. Our approach, NCRS, unifies these two
methods by having two phases. The first phase is generative, and the second
one is an agent’s policy.

##  3 Approach: A Unified Substrate

The modular robots grown by the NCA consist of different types of cells such
as sensors, actuators, and connecting tissue. After growth, the robot is
deployed in its particular environment. Importantly, in our approach, the same
NCA controls both the growth of the modular robot (Fig. 1a) and the robot
itself (Fig. 1b). Therefore, it is a unified substrate for body-brain co-
design and is called Neural Cellular Robot Substrate (NCRS). The architecture
of NCRS is illustrated in Fig. 1c. When the growth process is finished, the
channels responsible to define the body modules reflect the robot’s
morphology, then the NCA can observe and act in the environment using the
cells assigned to the specific types of modules, which are sensors, wheels,
and tissue.

The state of a cell is updated considering the eight surrounding neighbors and
itself, then it forms a $3\times 3$ neighborhood. The values of the nine cells
with all the $n$ channels are processed by a trainable convolutional layer
with 30 filters of size $3\times 3\times n$. Followed by a dense layer of 30
neurons and another one with $n$ neurons for the $n$ channels of the neural
CA. After all cells have been computed, the result of this process is added to
the previous state of the neural CA, and then it is clipped to the range of
$[-5,5]$. This update is only valid for the cells that are considered "alive",
which are the ones that have their value in the body channel greater than
$0.1$ and their neighbors. This architecture is very similar to the ones in
self-classifying MNIST [28] and in self-organized control of a cart-pole agent
[36].

The channels have specific roles in the neural CA, as shown in Fig. 3. The
number of channels $n$ differs because of the different number of sensors in
the types of benchmark tasks. The body channel is the one that indicates that
there is a body part in that cell if its value is greater than $0.1$. The
neighbors of a body part are allowed to update their states because they are
considered "growing". The next channel has fixed values and works as a control
flag. When the neural CA is in the developmental phase, all cells in this
channel are set to zero. When it is in the control phase, they are set to one.
The following channels are responsible to define the type of the body part.
The channel with the highest value is the one that specifies the body part. In
the case of a tie, the first channel is selected. The order of those channels
is: body/brain tissue, light/ball sensor, target area sensor (if needed), and
wheel. In this way, it can define a robot as depicted in Fig. 1a. Then, there
are the hidden channels to support the computation in the neural CA. For all
benchmark types, the neural CA contains six hidden channels. Finally, the
input/output channel, which is the one that receives the values from the
sensors and gives the values to the actuators (wheels).

(a) Initial time-step in developmental phase (b) Final time-step of
developmental phase

(c) Final time-step of control phase

Figure 3: Channels of the neural cellular automaton in different stages.

The initial state of the neural CA is a "seed". The middle cell of the grid
has the state set as one in the body channel, and the rest is zero. Fig. 3(a)
illustrates this. After a few time-steps in the developmental phase, all
channels are updated except the control flag channel. This phase lasts for ten
time-steps. The end of the developmental phase is represented in Fig. 3(b).
After development, the control phase starts. In this phase, the benchmark
environment initializes with the developed robot body. For advancing one time-
step in the environment, the NCA takes two time-steps for defining an action
after receiving observations from the sensors. The body and body parts
channels become fixed and their values are defined by the robot body. This is
used to support the neural CA by identifying the cells with body modules, such
as tissues, sensors, and actuators. The cell is assigned the value one to the
body channel if there is a body part and to the specific body parts channel.
Fig. 3(c) shows this assignment for the identification of body parts during
the control phase. The robot designed by this NCA is depicted in Fig. 4(a). At
the start of the control phase, the cellular activity of the hidden and
input/output channels is set to zero. In the input/output channel, only the
input cells are fixed and their values come from the sensors.

In our neural CA, there is no noise. Even all "alive" cells are updated every
time-step. This is done because the stochastic update or any other type of
noise would affect the development of the robot body. After the developmental
phase, the same model could produce different types of robot body.

For our experiments, the neural cellular automaton has a grid of size $5\times
5$. Therefore, it generates a body for an agent with the same size. Since this
is a neural CA, the grid size does not affect the number of trainable
parameters. The light chasing and light chasing with obstacle environments
require just the light sensor. Therefore, the robot can have tissue, light
sensor and wheel. A wheel’s orientation is always vertical during the
initialization of the benchmark environment. The wheel rotates upwards and
downwards relative to the initial angle of the robot. The maximum speeds for
each of those directions are, respectively, +1 and -1. This takes three body
part channels. With one body, one control flag, six hidden, and one
input/output channels, the total number of channels is $12$. In this way, the
number of trainable parameters is 4,572. In the carrying ball to target
environment, the robot needs one additional sensor. Therefore, it adds one
more channel. It results in a neural CA with 4,873 trainable parameters.

##  4 Benchmark environments

To test the capacity of controlling the developed robot, we implemented three
benchmark environments, which are: light chasing (LC), light chasing with
obstacle (LCO), and carrying ball to target (CBT). They are environments where
a modular robot equipped with simple light sensors and wheels can be
evaluated. In those environments, we decided that the size of the playing
field and the distance between the objects are affected by the maximum size
that the robot can have. Thus, the larger the robot can be, the bigger the
playing field. In our experiments, we use a robot and a neural cellular
automaton grid with size $5\times 5$. Because the possible maximum size of the
robot is 5, we chose the size of the playing field to be 60.

The fitness score is calculated using the average score of 12 runs where the
location of the agent, light, ball, and target can differ for each run. The
light or ball has some predefined regions to be initially placed.

The benchmark environments are based on the implementation of the top-down
racing environment in Open AI gym [4]. We use the pybox2d, which is a 2D
physics library in Python.

###  4.1 Light chasing

The light chasing (LC) environment is shown in Fig. 1b. The goal of the agent
is to be closer to the light during the entire simulation. The agent starts in
the middle of the playing field. One light is randomly placed around the
region of one of the four corners of the playing field. The fitness score is
calculated by the average distance between the center of the robot and the
center of the light over all simulation time-steps, and a successful run means
that this distance reached less than 10 times the module size. The activity
$s$ of the agent’s light sensors is calculated as:

$s=e^{-distance/playfield},$ (1)

where the distance between the objects is normalized by the size of the
playing field $playfield$, which is 60. The values of the sensor activity or
fitness score are between 0 and 1, where 1 means no distance. The values
exponentially decay to 0 with an increase in distance.

###  4.2 Light chasing with obstacle

The light chasing with obstacle (LCO) environment is a more difficult version
of the light chasing one (Fig. 2(a)). The robot does not have sensors to
detect the obstacle, thus its morphology plays a bigger role in this
benchmark. The passage width is calculated by the possible maximum size of the
robot. If the robot can have up to $5\times 5$ body parts, then the passage
width would be the size of three body parts. The robot is randomly initialized
at the bottom of the playing field. An obstacle is procedurally generated with
a target passage width and wall roughness. The obstacle has the shape of a
funnel because there are no sensors to it, then this helps the robot to reach
the passage depending on its body. The passage is randomly located on the
horizontal axis and fixed on the vertical axis. The light is at the top and
after the obstacle. The initial light location has four predefined regions on
the horizontal axis, which are left, center-left, center-right and right. The
fitness and success definition are the same as the light chasing task.

###  4.3 Carrying ball to target area

Among the three benchmark environments, the task to carry a ball to a target
area is the most difficult one (Fig. 2(b)). For the control phase, the robot
needs to move towards the ball, and then move to the target area without
losing the ball during the transport. For the developmental phase, the body of
the robot needs to be adequate to push or kick the ball to the target area,
and properly placing the sensors of each type, so it can successfully locate
ball and target area. The agent is located at the bottom in a random
horizontal location. The ball is located in the middle of the vertical axis of
the playing field, but it has the same four predefined regions as the light
chasing with obstacle environment. The target is located at the top and its
location on the horizontal axis is randomly defined. Besides the sensor for
the ball (or light for the other two environments), there is a new sensor type
that calculates the distance to the center of the target area (following (1)).

The fitness score of this environment is the average of the distance between
robot and ball, and the distance between the ball and the center of the target
area. Since they are distances used to calculate the fitness score, they are
normalized using (1). The definition of success in this task means carrying
the ball to the target, so it can have a distance less than ten times the
module size of the robot.

##  5 Training methods

We have chosen to use some derivative-free optimization methods because NCRS
needs some adjustments for using deep reinforcement learning because of the
variable number of inputs and outputs [36]. They are the covariance matrix
adaptation evolution strategy (CMA-ES) [11] and covariance matrix adaptation
MAP-Elites (CMA-ME) [9]. The latter is used to add quality diversity to the
former, broadening the exploration of robot designs. For both training
methods, we use the library CMA-ES/pycma [12]. There are two training methods
and three benchmark tasks. This gives a total of six different combinations.
Because of the computational demands, each of these combinations was trained
only once.

The training process is performed entirely on a CPU. To speed up evaluation
times, robots with a design that would not work properly in the environment
are not simulated. For the two light chasing environments, robots must have at
least one light sensor and two actuators. For the carrying ball to target,
they must have one sensor of each type and two actuators. The fitness scores
of the failed designs are calculated according to the number of correct parts
they have. For each correct body part, the fitness score increases by $0.01$.

To compare the quality diversity of CMA-ES and CMA-ME, we use the percentage
of cells or feature configurations filled, and QD-score. They measure quality
and diversity of the elites [27]. The QD-score is calculated by summing the
fitness score of all elites and dividing it by the total number of possible
feature configurations. Moreover, CMA-ES and CMA-ME have their elites stored,
even though CMA-ES does not use elites during training.

###  5.1 Covariance matrix adaptation evolution strategy

CMA-ES is one of the most effective derivative-free numerical optimization
methods for continuous domains [9]. CMA-ES runs 20,000 generations for all
environments. The initial mean is $0.0$ for all dimensions, and the initial
coordinate-wise standard deviation (step size) is $0.01$. The population size
or the number of solutions acquired to update the covariance matrix is 112.
This number was selected by the number of available threads in the machine
used to train, which contains 56 threads at 2.70GHz.

###  5.2 Covariance matrix adaptation MAP-Elites

CMA-ME is a variant of CMA-ES with the added benefit of quality diversity from
MAP-Elites [24]. The changes to CMA-ES are that there are emitters of CMA-ES
being trained in a cycle. Additionally, a feature map stores one elite for
each possible feature configuration. Because there are invalid body designs,
they do not produce an elite. When there is a successful robot design, the
number of sensors, actuators, and body parts are used as features. If there
are no elites or the current solution is better than the actual elite stored
in the feature map, then the current solution is assigned to its feature
configuration.

We use a slightly modified version of the CMA-ME with improvement emitters
[9]. We only restart an emitter when the number of elites is greater than the
number of emitters and it is stuck for more than 500 generations. Being stuck
means that the emitter could not find a better elite or an elite could not be
placed into an empty feature configuration in the map. When an emitter
restarts, the mean used to initialize the CMA-ES is a random elite in the map.

CMA-ME is executed for 60,000 generations for all environments, except the
light chasing environment with 67,446 generations because we forced it to stop
a longer training and its best fitness score was already better than the one
trained with CMA-ES. The initial mean and the initial coordinate-wise standard
deviation are the same as CMA-ES for all emitters. The population size is 128
because the CMA-ME training was executed in a computer with 128 threads at
2.9GHz.

##  6 Results

The training process took around 2.5 days for optimizing the NCA with CMA-ES.
The evolution with CMA-ME took around 5.5 days. It is important to note that
they do not have the same machine configuration, population size, and maximum
number of generations.

Fig. 4 shows all robot designs with the best fitness scores in regards to
their training method and task. Almost all robots for the LC and CBT tasks
fill the entire $5\times 5$ grid of cells. Those environments do not have any
environmental constraints (any obstacle) for the robot size. Therefore, we
infer that the full grid of modules is easier to design and there are more
computational resources for controlling the robot. Their fitness scores are
shown in Table 1. The results indicate that CMA-ES and CMA-ME can reach almost
the same fitness scores after training. However, CMA-ME has fewer generations
for the 15 emitters (4,000 generations per emitter). It is possible that if we
run 20,000 generations per emitter, CMA-ME could reach a better final
performance than CMA-ES and with more diversity. The history of the maximum
fitness score per generation is depicted in Fig. 5.

Table 1: Best fitness score after training in the tasks of light chasing (LC), light chasing with obstacle (LCO) and carrying ball to target (CBT) | CMA-ES | CMA-ME
---|---|---
LC | 0.58274 | 0.61481
LCO | 0.49295 | 0.47723
CBT | 0.48445 | 0.47884

(a) CMA-ES - LC (b) CMA-ES - LCO (c) CMA-ES - CBT (d) CMA-ME - LC (e) CMA-ME -
LCO (f) CMA-ME - CBT

Figure 4: Robot designs with best fitness scores for the tasks of light
chasing (LC), light chasing with obstacle (LCO) and carrying ball to target
(CBT).

(a) Light chasing (b) Light chasing with obstacle (c) Carrying ball to target

Figure 5: Maximum fitness score through generations.

The elites were saved for both CMA-ES and CMA-ME, then we can compare their
quality diversity. In Table 2, the number of cells filled and QD-scores of all
six methods and tasks combinations are presented. It is noticeable that CMA-ME
provides much more quality diversity because of its bigger number of feature
configurations and its QD-score. We can visualize it in Fig. 6. This shows a
small part of the elites produced for the light chasing tasks with CMA-ES and
CMA-ME. Nevertheless, it confirms those two quality diversity measurements
because more cells are filled, and there are more cells with higher fitness
scores.

Table 2: Elites stored during training for the light chasing (LC), light chasing with obstacle (LCO) and carrying ball to target (CBT) | CMA-ES | CMA-ME
---|---|---
| Cells filled | QD-score | Cells filled | QD-score
LC | 67.58% | 0.29530 | 89.57% | 0.40152
LCO | 17.88% | 0.06069 | 61.80% | 0.19841
CBT | 58.10% | 0.22957 | 93.82% | 0.37996
Figure 6: Selected elites trained in the light chasing environment. Those
modules were selected because they are the most different between CMA-ES and
CMA-ME. Axes and subplots indicate the number of components.

For testing the success of our six trained models, we run 100 times the
simulation and the percentage of success is presented in Table 3. We can
visualize some examples of those simulations in Fig. 7. The trained model with
CMA-ES for the light chasing task got 92% of success rate with $0.58274$
fitness score while the one trained with CMA-ME had 75% of success and
$0.61481$ of fitness score. This means that a higher fitness score does not
indicate a more successful model for reaching the light. This can be observed
in Fig. 7(a)-7(d) for CMA-ES, and Fig. 7(e)-7(h) for CMA-ME. We can see in
Fig. 7(g) that the light is at the top-right corner and the robot goes to the
top-left corner. This explains the 75% success rate of this NCRS because the
light is at the top-right corner in 25 out of the 100 simulations. This model
learned to move faster to the light in the other three corners, but it misses
the one in the top-right corner. For the light chasing with obstacle task, the
reason for the higher success rate of CMA-ES robot is that it is much thinner
than the CMA-ME robot. Therefore, it is easier to pass through the passage. If
we define success in LCO by passing the center of the body through the
passage, then CMA-ES and CMA-ME had a success rate, respectively, of 77% and
45%. The NCRS did not learn to move to the light after passing through the
obstacle. It just moves forward. Because of the difficulty of this task, we
can consider the results for LCO were partially successful in general and
successful in body design. Fig. 7(i)-7(l) and Fig. 7(m)-7(p) show that. The
task of carrying a ball to a target had no successful trained model. The
robots for both training methods just move forward and, by chance, it moves
the ball to target. This can be seen in Fig. 7(q)-7(t) and Fig. 7(u)-7(x).

Fig. 3 shows how the channels progress through time. The hidden channels are
predominantly different in their behavior for the developmental and control
phases. We infer this is mainly due to the control flag channel which
regulates these two phases. We can observe the different patterns that emerged
in their final time-steps. From the initial "seed" state to the state in Fig.
3(b), we can see how the NCA behaves during 10 time-steps of the developmental
phase. In Fig 3(c), we can see the end of the control phase during its 200
time-steps (100 time-steps in the environment). We can still understand its
behavior because the hidden and input/output channels were set to zero at the
beginning of the control phase, and the body, control flag, tissue, sensor,
and actuator channels were fixed according to the morphology of the robot.

Table 3: Testing success percentage over 100 runs for the tasks of light chasing (LC), light chasing with obstacle (LCO) and carrying ball to target (CBT) | CMA-ES | CMA-ME
---|---|---
LC | 92% | 75%
LCO | 20% | 8%
CBT | 1% | 2%

(a) CMA-ES - LC #1 (b) CMA-ES - LC #2 (c) CMA-ES - LC #3 (d) CMA-ES - LC #4

(e) CMA-ME - LC #1 (f) CMA-ME - LC #2 (g) CMA-ME - LC #3 (h) CMA-ME - LC #4

(i) CMA-ES - LCO #1 (j) CMA-ES - LCO #2 (k) CMA-ES - LCO #3 (l) CMA-ES - LCO
#4

(m) CMA-ME - LCO #1 (n) CMA-ME - LCO #2 (o) CMA-ME - LCO #3 (p) CMA-ME - LCO
#4

(q) CMA-ES - CBT #1 (r) CMA-ES - CBT #2 (s) CMA-ES - CBT #3 (t) CMA-ES - CBT
#4

(u) CMA-ME - CBT #1 (v) CMA-ME - CBT #2 (w) CMA-ME - CBT #3 (x) CMA-ME - CBT
#4

Figure 7: Last time-step where the robot is fully visible of the best NCRS
trained with CMA-ES and CMA-ME in the environments for light chasing (LC),
light chasing with obstacle (LCO) and carrying ball to target (CBT).

##  7 Discussion and conclusion

Body-brain co-evolution is a challenging task [2]. In this work, we developed
three benchmark tasks for robot co-design and introduced a novel method by
having a unified substrate as a genome with its own rules. This substrate is a
single neural cellular automaton that works to develop and control a modular
robot. This novelty opens up several possibilities in open-ended evolution
[32], especially because body and brain can co-evolve to the limits of the
capacity of the artificial neural network. Because it defines the local rules
in the CA, NCRS has the advantage of scalability. We also infer that
curriculum learning will be important for complexifying the evolving robot
[1]. For example, the number of body parts and dimensions can increase over
time with the progress of the generations. Evolution in multi-agent
environments may also be applied, such as in PolyWorld [39]. We can also try
to remove the two separated phases into one. Thus, we can observe how
development and control can emerge and the performance the modular robots can
have.

The presented results were successful for the LC task, but our trained models
presented some failures when increasing the difficulty of the tasks. This may
be addressed by adjusting the fitness score to reflect the success conditions,
as well as by applying curriculum learning [1]. In future works, we plan to
apply our method in the Evolution Gym [2], or in a modified version of
VoxCraft [21] for 3D soft robots. Moreover, we aim at training and testing our
approach for self-repair and robustness to noise.

## Acknowledgment

This work was partially funded by the Norwegian Research Council (NFR) through
their IKTPLUSS research and innovation action under the project Socrates
(grant agreement 270961). We thank Henrique Galvan Debarba for his thoughtful
comments about the text. We also thank Joachim Winther Pedersen, Djordje
Grbic, Miguel González Duque, and Rasmus Berg Palm for the helpful discussions
during the implementation of the experiments.

## References

  * Bengio et al. [2009] Yoshua Bengio, Jérôme Louradour, Ronan Collobert, and Jason Weston.  Curriculum learning.  In _Proceedings of the 26th annual international conference on machine learning_ , pages 41–48, 2009.
  * Bhatia et al. [2021] Jagdeep Bhatia, Holly Jackson, Yunsheng Tian, Jie Xu, and Wojciech Matusik.  Evolution gym: A large-scale benchmark for evolving soft robots.  _Advances in Neural Information Processing Systems_ , 34, 2021.
  * Bongard and Pfeifer [2003] Josh C Bongard and Rolf Pfeifer.  Evolving complete agents using artificial ontogeny.  In _Morpho-functional Machines: The new species_ , pages 237–258. Springer, 2003.
  * Brockman et al. [2016] Greg Brockman, Vicki Cheung, Ludwig Pettersson, Jonas Schneider, John Schulman, Jie Tang, and Wojciech Zaremba.  Openai gym, 2016.
  * Brodbeck et al. [2015] Luzius Brodbeck, Simon Hauser, and Fumiya Iida.  Morphological evolution of physical robots through model-free phenotype development.  _PloS one_ , 10(6):e0128444, 2015.
  * Conway et al. [1970] John Conway et al.  The game of life.  _Scientific American_ , 223(4):4, 1970.
  * Doursat et al. [2013] René Doursat, Hiroki Sayama, and Olivier Michel.  A review of morphogenetic engineering.  _Natural Computing_ , 12(4):517–535, 2013.
  * Earle et al. [2021] Sam Earle, Justin Snider, Matthew C Fontaine, Stefanos Nikolaidis, and Julian Togelius.  Illuminating diverse neural cellular automata for level generation.  _arXiv preprint arXiv:2109.05489_ , 2021.
  * Fontaine et al. [2020] Matthew C Fontaine, Julian Togelius, Stefanos Nikolaidis, and Amy K Hoover.  Covariance matrix adaptation for the rapid illumination of behavior space.  In _Proceedings of the 2020 genetic and evolutionary computation conference_ , pages 94–102, 2020.
  * Gupta et al. [2021] Agrim Gupta, Silvio Savarese, Surya Ganguli, and Li Fei-Fei.  Embodied intelligence via learning and evolution.  _Nature communications_ , 12(1):1–12, 2021.
  * Hansen and Ostermeier [1996] Nikolaus Hansen and Andreas Ostermeier.  Adapting arbitrary normal mutation distributions in evolution strategies: The covariance matrix adaptation.  In _Proceedings of IEEE international conference on evolutionary computation_ , pages 312–317. IEEE, 1996.
  * Hansen et al. [2019] Nikolaus Hansen, Youhei Akimoto, and Petr Baudis.  Cma-es/pycma on github.  _Zenodo, doi_ , 10, 2019.
  * Harding and Banzhaf [2009] Simon Harding and Wolfgang Banzhaf.  Artificial development.  In _Organic Computing_ , pages 201–219. Springer, 2009.
  * Hejna III et al. [2021] Donald J Hejna III, Pieter Abbeel, and Lerrel Pinto.  Task-agnostic morphology evolution.  _arXiv preprint arXiv:2102.13100_ , 2021.
  * Horibe et al. [2021] Kazuya Horibe, Kathryn Walker, and Sebastian Risi.  Regenerating soft robots through neural cellular automata.  In _EuroGP_ , pages 36–50, 2021.
  * Komosiński and Ulatowski [1999] Maciej Komosiński and Szymon Ulatowski.  Framsticks: Towards a simulation of a nature-like world, creatures and evolution.  In _European Conference on Artificial Life_ , pages 261–265. Springer, 1999.
  * Kriegman et al. [2018] Sam Kriegman, Nick Cheney, and Josh Bongard.  How morphological development can guide evolution.  _Scientific reports_ , 8(1):1–10, 2018.
  * Langton [2019] Christopher Langton.  _Artificial life: proceedings of an interdisciplinary workshop on the synthesis and simulation of living systems_.  Routledge, 2019.
  * Li and Yeh [2002] Xia Li and Anthony Gar-On Yeh.  Neural-network-based cellular automata for simulating multiple land use changes using gis.  _International Journal of Geographical Information Science_ , 16(4):323–343, 2002.
  * Lipson and Pollack [2000] Hod Lipson and Jordan B Pollack.  Automatic design and manufacture of robotic lifeforms.  _Nature_ , 406(6799):974–978, 2000.
  * Liu et al. [2020] S Liu, D Matthews, S Kriegman, and J Bongard.  Voxcraft-sim, a gpuaccelerated voxel-based physics engine, 2020.
  * Medvet et al. [2021] Eric Medvet, Alberto Bartoli, Federico Pigozzi, and Marco Rochelli.  Biodiversity in evolved voxel-based soft robots.  In _Proceedings of the Genetic and Evolutionary Computation Conference_ , pages 129–137, 2021.
  * Mordvintsev et al. [2020] Alexander Mordvintsev, Ettore Randazzo, Eyvind Niklasson, and Michael Levin.  Growing neural cellular automata.  _Distill_ , 5(2):e23, 2020.
  * Mouret and Clune [2015] Jean-Baptiste Mouret and Jeff Clune.  Illuminating search spaces by mapping elites.  _arXiv preprint arXiv:1504.04909_ , 2015.
  * Nichele et al. [2017] Stefano Nichele, Mathias Berild Ose, Sebastian Risi, and Gunnar Tufte.  Ca-neat: evolved compositional pattern producing networks for cellular automata morphogenesis and replication.  _IEEE Transactions on Cognitive and Developmental Systems_ , 10(3):687–700, 2017.
  * Niklas and Newman [2013] Karl J Niklas and Stuart A Newman.  The origins of multicellular organisms.  _Evolution & development_, 15(1):41–52, 2013\.
  * Pugh et al. [2016] Justin K Pugh, Lisa B Soros, and Kenneth O Stanley.  Quality diversity: A new frontier for evolutionary computation.  _Frontiers in Robotics and AI_ , 3:40, 2016.
  * Randazzo et al. [2020] Ettore Randazzo, Alexander Mordvintsev, Eyvind Niklasson, Michael Levin, and Sam Greydanus.  Self-classifying mnist digits.  _Distill_ , 5(8):e00027–002, 2020.
  * Sims [1994] Karl Sims.  Evolving 3d morphology and behavior by competition.  _Artificial life_ , 1(4):353–372, 1994.
  * Slack and Dale [2021] Jonathan MW Slack and Leslie Dale.  _Essential developmental biology_.  John Wiley & Sons, 2021.
  * Stanley [2007] Kenneth O Stanley.  Compositional pattern producing networks: A novel abstraction of development.  _Genetic programming and evolvable machines_ , 8(2):131–162, 2007.
  * Stanley [2019] Kenneth O Stanley.  Why open-endedness matters.  _Artificial life_ , 25(3):232–235, 2019.
  * Sudhakaran et al. [2021] Shyam Sudhakaran, Djordje Grbic, Siyan Li, Adam Katona, Elias Najarro, Claire Glanois, and Sebastian Risi.  Growing 3d artefacts and functional machines with neural cellular automata.  _arXiv preprint arXiv:2103.08737_ , 2021.
  * Talamini et al. [2021] Jacopo Talamini, Eric Medvet, and Stefano Nichele.  Criticality-driven evolution of adaptable morphologies of voxel-based soft-robots.  _Frontiers in Robotics and AI_ , 8:172, 2021.
  * Topa [2011] Pawel Topa.  Network systems modelled by complex cellular automata paradigm.  _Cellular automata-simplicity behind complexity. InTech, Europe, China_ , pages 259–274, 2011.
  * Variengien et al. [2021] Alexandre Variengien, Sidney Pontes-Filho, Tom Eivind Glover, and Stefano Nichele.  Towards self-organized control: Using neural cellular automata to robustly control a cart-pole agent.  _Innovations in Machine Intelligence_ , 1:1–14, Dec 2021\.  doi: 10.54854/imi2021.01.
  * Veenstra and Glette [2020] Frank Veenstra and Kyrre Glette.  How different encodings affect performance and diversification when evolving the morphology and control of 2d virtual creatures.  In _ALIFE: proceedings of the artificial life conference_ , pages 592–601. MIT Press, 2020.
  * Wolfram [2002] Stephen Wolfram.  _A new kind of science_ , volume 5.  Wolfram media Champaign, IL, 2002.
  * Yaeger et al. [1994] Larry Yaeger et al.  Computational genetics, physiology, metabolism, neural systems, learning, vision, and behavior or poly world: Life in a new context.  In _SANTA FE INSTITUTE STUDIES IN THE SCIENCES OF COMPLEXITY-PROCEEDINGS VOLUME-_ , volume 17, pages 263–263. Citeseer, 1994.
  * Yim et al. [2007] Mark Yim, Wei-Min Shen, Behnam Salemi, Daniela Rus, Mark Moll, Hod Lipson, Eric Klavins, and Gregory S Chirikjian.  Modular self-reconfigurable robot systems [grand challenges of robotics].  _IEEE Robotics & Automation Magazine_, 14(1):43–52, 2007.
  * Zador [2019] Anthony M Zador.  A critique of pure learning and what artificial neural networks can learn from animal brains.  _Nature communications_ , 10(1):1–7, 2019.


# Towards Fully Passive Plug and Play Time-Bin Quantum Key Distribution over
Multi-Mode Channels

Ramy Tannous1, Wilson Wu1, Stéphane Vinet1, Chithrabhanu Perumangatt2, Dogan
Sinar1, Alexander Ling2,3, Thomas Jennewein1 1 Institute for Quantum Computing
and Department of Physics & Astronomy, University of Waterloo, 200 University
Ave W, Waterloo, Ontario, N2L 3G1, Canada  2 Centre for Quantum Technologies,
3 Science Drive 2, National University of Singapore 117543 Singapore  3
Department of Physics, Faculty of Science, National University of Singapore, 2
Science Drive 3, 117551 Singapore

###### Abstract

Phase stabilization of distant quantum time-bin interferometers can be a major
challenge for quantum communication networks, and is typically achieved by
exchanging reference optical signals, which can be particularly challenging
for free-space communication channels. We demonstrate a novel solution that
utilizes reference frame independent time-bin quantum key distribution that
completely avoids the need for active relative phase stabilization of the
remote interferometers. We realized a proof of concept experiment using a
source of hybrid entanglement with polarization and time-bin, and show a
sustained asymptotic secure key rate of greater than
$0.06\text{\,}\mathrm{b}\mathrm{i}\mathrm{t}\mathrm{s}\mathrm{/}\text{coincidence}$
over a $15\text{\,}\mathrm{m}$ multi-mode fiber optical channel. This is
achieved without the use of mode filtering, adaptive optics, active basis
selection, nor active phase alignment. This system can therefore be considered
a passive plug and play self-compensating solution that is suitable for
spatially multi-mode channels, and improves the viability of time-bin quantum
communication links over challenging channels such as airborne, ground to
space, and multi-mode fibers.

Time-bin encoding for quantum key distribution (QKD) continues to be the
encoding of choice for optical-fiber quantum networks [1, 2, 3]. More
recently, there have been demonstrations of time-bin over free space channels
[4], enabled by the development of multi-mode field widened interferometers
[5]. Despite recent efforts in improving the short and long-term phase
stability of these interferometers [6], it remains experimentally challenging.
Active systems are typically required to align the sender and receiver
interferometers’ relative phase which creates additional technical overhead to
successfully perform time-bin QKD [7, 8, 9]. In addition, using active phase
alignment over a fluctuating channel, such as a turbulent atmosphere, is
difficult since a reference laser would suffer from random intensity
fluctuations. Conversely, plug & play QKD systems provide a passive self-
compensating solution and are typically based on retroreflection [10, 11].
However, we make use of a reference frame independent quantum key distribution
(RFI-QKD) protocol that avoids the requirement for active phase locking of the
two interferometers, and is applicable to both prepare-and-measure and
entanglement schemes [12]. The specific benefit of RFI-QKD protocols is that
while the relative phase between the two interferometers is allowed to drift,
the channel purity between the two key exchange parties can still be
determined, provided the phase drift is relatively “slow” to allow for the
sufficient estimation of the channel parameters. Recently, several time-bin
RFI-QKD demonstrations were performed over free space and optical fiber
channels [13, 14, 15]. However, the use of single-mode interferometers and
active basis choice have their disadvantages. The use of single-mode
interferometers over a free space channel can result in significant photon
throughput losses, and thus affect the overall system performance, unless an
adaptive optics system is used which increases the system overhead. In
addition, active basis and phase encoding requires effort in pulse
synchronization that can present experimental challenges for moving platforms
such as planes [16], drones, or satellites [17, 18], where resources may be
limited.

In this work, we demonstrate a fully passive time-bin quantum communication
scheme that has no active basis selection, mode filtering, adaptive optics,
nor interferometer stabilization. To achieve this, we utilize the entanglement
as the source of randomness reducing the required operational overhead of the
implementation, and use interferometers shown to easily perform time-bin QKD
over a free space channel [4]. While the use of the 6-state 4-state time-bin
RFI-QKD removes the need for active interferometer phase stabilization. Due to
the simplicity of our scheme, it can be easily adapted for free space quantum
channels particularly moving platforms. Despite the use of a spatially multi-
mode channel, our system is able to achieve quantum bit error rates (QBER)
below $5\%$ and $7\%$ for the key distillation basis and phase estimation
basis, respectively, while maintaining a continuous asymptotic key rate of
greater than
$0.06\text{\,}\mathrm{b}\mathrm{i}\mathrm{t}\mathrm{s}\mathrm{/}\text{coincidence}$.

We utilize the 6-state 4-state time-bin RFI-QKD protocol that is based on an
entangled photon source (EPS) [19] where one photon is converted from
polarization to time-bin [20]. The state shared between Alice and Bob is thus

$\ket{\psi_{\text{AB}}}=\frac{1}{\sqrt{2}}\left(\ket{H\mathcal{E}}+e^{i\phi_{\text{AB}}}\ket{V\mathcal{L}}\right),$
(1)

where $\phi_{AB}$ is the relative phase been Alice and Bob; $H$, $V$ are the
horizontal and vertical polarizations respectively; and $\mathcal{E}$,
$\mathcal{L}$ are the early and late time-bins respectively. The particular
benefit of this hybrid entangled state is the relative ease of producing a
6-state analyzer for the polarization encoded photons, while a 4-state
analyzer for time-bin photons is an unbalanced interferometer. Alice performs
a 6-state polarization measurement on the first qubit, measuring in the
$Z_{A}\text{, }X_{A}\text{, and }Y_{A}$ bases, while Bob makes a 4-state time-
bin measurement on the other qubit measuring in the in the $Z_{B}\text{ and
}X_{B}$ bases. In our implementation, the computational basis is the
horizontal-vertical polarization basis ($Z_{A}$) for Alice and the early-late
($Z_{B}$) basis for Bob, and are the fixed bases for the RFI protocol that are
used to generate the keymap.111Although we use an entanglement based approach,
this scheme could easily be adapted for a prepare-and-measure type scheme
where the sender produces 6 states while the receiver measures 4 states or
vice-versa.

As shown in [19], a key can be generated despite the phase drifts provided
that strong correlations are maintained in between the polarization
superposition bases and the time-bin superposition basis, which are found in
the $C$-parameter given by

$C=\sqrt{\braket{X_{A}\otimes X_{B}}^{2}+\braket{Y_{A}\otimes X_{B}}^{2}},$
(2)

where $0\leq C\leq 1$ allows the users to estimate the channel quality with
$C=1$ indicating a maximally entangled state or perfect channel [19].
Monitoring the superposition coincidences between the polarization
superposition bases and the time-bin superposition basis ($C$-parameter)
allows the users to estimate the channel quality.

During key exchange, $\phi_{AB}$ can take any value and can drift during key
generation and exchange, provided it does not change too rapidly [21, 22]. The
degree to which the phase can change is a function of the photon detection
rate and timing resolution of the system. To quantify the effects of a rapidly
changing $\phi_{AB}$ we use a Monte Carlo simulation. We consider $N$ signals
sent by the EPS to Alice and Bob that arrive evenly spaced over a time period
$t_{N}$. The instantaneous phase of the entangled state is $\phi(\tau n)$,
where $n\in[0,N]$, while the total phase variation over the measurement time
interval is given by $\phi(t_{N})$. Numerically using the Qutip package [23,
24], we produce $N$ entangled states (1) each with a phase $\phi(\tau n)$.
Then for each $n$ entangled state, we calculate the detection probabilities of
the POVM’s (e.g. $XX_{++}$, $XX_{+-}$, $XX_{-+}$, $XX_{--}$) for each basis
that is involved in calculating the $C$-parameter. Since $\phi(t)$ is
changing, these probabilities should vary with time (and phase). With the
probabilities calculated from every $n$ entangled state, a pseudorandom number
generator determines the detection results of a single photon event. This is
done independently for each event. The result is then counted and statistics
are accumulated for the total $N$ entangled pairs. Once complete, an
experimental expectation value approximation is calculated using equation (2)
of [19] for each basis from which a value for the average $C$-parameter over
the time interval is determined. This process is repeated $i$ times for each
$\phi(t_{N})$ and the average value of the distribution is taken to be the
value of $C(t_{N})$ with error bounds being the standard deviation of the each
distribution.

Figure 1: Experimental setup of the time-bin RFI-QKD proof of principle
demonstration. An entangled photon source (EPS) creates polarization entangled
photon pairs at $785\text{\,}\mathrm{nm}$ (signal) and
$842\text{\,}\mathrm{nm}$ (idler) via spontaneous parametric downconversion.
The idler photon is measured at a local 6-state polarization analyzer, while
the signal photon is converted to time-bin encoding by the polarization to
time-bin converter (PTC) and sent across a $15\text{\,}\mathrm{m}$ multi-mode
fiber (MMF) optic channel to a receiver, the time-bin analyzer (TA), that
performs a 4-state time-bin measurement. A time tagging unit is used to record
detection events (not shown). Both the PTC and TA interferometers have a path
separation of $2.2\text{\,}\mathrm{n}\mathrm{s}$. The TA is a $4f$ field
widened interferometer with multi-mode fiber inputs and outputs. No active
phase compensation is used to maintain a phase reference between the two
interferometers. A piezoelectric actuator is inserted in the TA short path to
induce a phase drift. HWP: half-waveplate, QWP: quarter-waveplate, L: lens,
POL: polarizer, PMF: polarization maintaining fiber, PBS: polarizing beam
splitter, BS: 50:50 beam splitter.

Figure 1 shows a sketch of the experimental setup used for the time-bin RFI-
QKD proof of principle demonstration, where we perform the adapted 6-state
4-state protocol [19]. Our experiment consists of an entangled photon source
[25], a local 6-state polarization measurement, a polarization to time-bin
converting interferometer, a quantum channel, a time-bin analyzing
interferometer, and detection of single photons, including time-tagging.
Coincidence and correlation data is done on a computer during post-processing.
The EPS creates polarization entangled photon pairs in the state
$\ket{\psi_{AB}}=\frac{1}{\sqrt{2}}\left(\ket{H_{A}H_{B}}+e^{i\phi_{AB^{*}}}\ket{V_{A}V_{B}}\right)$,
where $\phi_{AB^{*}}$ is a random phase introduced by the EPS. From each
entangled pair, the idler photon is measured using a 6-state polarization
measurement [19], while the signal photon is converted to time-bin encoding
before being sent across a $15\text{\,}\mathrm{m}$ multi-mode fiber. The
polarization to time-bin converter (PTC) is an unbalanced Michelson
interferometer shown in figure 1. A polarizing beam splitter in combination
with quarter-waveplates are used to maximize the signal sent to the output
port of the Michelson interferometer. The PTC is designed such that the
horizontal polarized light takes the short path $(H\rightarrow\mathcal{E})$
and vertically polarized light takes the long path,
$(V\rightarrow\mathcal{L})$ which creates the entangled state in (1). Thus any
photon in a polarization superposition would create a superposition of both
interferometer paths in the PTC. The polarization information of the photons
is removed at the output of the PTC by placing a polarizer at $None$ with
respect to the horizontal and vertical polarizations. The separation of the
time-bins is given by the path difference of the PTC which is approximately
$2.2\text{\,}\mathrm{ns}$.

By passing through the PTC, the polarized signal photons are converted to
appropriate time-bin states and are then sent across the
$15\text{\,}\mathrm{m}$ multi-mode fiber (Corning InfiniCor 300 62.5/125). The
channel is a graded-index fiber which has a core size of
$62.5\text{\,}\mathrm{\SIUnitSymbolMicro m}$ and a numerical aperture of
$0.275$. The estimated maximum modal dispersion for this fiber is
$0.353\text{\,}\mathrm{ns}\text{\,}{\mathrm{km}}^{-1}$ which amounts to a
negligible $5.2\text{\,}\mathrm{ps}$ over the $15\text{\,}\mathrm{m}$ length.
A field widened multi-mode time-bin analyzer (TA) is used to measure the time-
bin encoded photons after they pass through the quantum channel [4]. The
optical layout of the TA is shown in figure 1. The TA uses a $4f$ imaging
system in each arm. No phase stabilization is used to obtain a phase reference
between the PTC and the TA. In fact, the relative phase of the PTC and TA are
allowed to drift with ambient noise. In some photon exchanges, a piezoelectric
actuator is used to deliberately induce a phase change in the TA to
demonstrate the robustness of the protocol (figure 3). The outputs of the TA
are coupled into multi-mode fibers connected to single photon avalanche
diodes. A time tagging unit is used to record the detection events and
coincidence analysis is done to calculate the correlations in the various
bases. The use of a single time tagging unit is not secure and two independent
time tagging units, one each at Alice and Bob, are needed to properly perform
QKD. However, the single time tagging unit is sufficient for this
demonstration.

The overall number of coincidences measured in our proof of concept experiment
is about $9000\text{\,}\text{coincidence}\mathrm{/}\mathrm{s}$. The EPS has an
intrinsic pair production rate of approximately
$10\text{\times}{10}^{6}\text{\,}\text{pairs}\mathrm{/}\mathrm{s}$. However,
the heralding efficiency of the source $\approx 20\%$, which includes the
losses due to fiber coupling into the PMF, reduces the measured pair rate.
Furthermore, the losses contributed by the PTC and TA, and the detector
efficiencies also reduce the overall number of coincidences measured. The PTC
throughput is measured to be $35.4\%$ including fiber coupling into the output
PMF, however it is intrinsically reduced by $50\%$ due to the polarizer that
eliminates the polarization signature of each time-bin. The
$15\text{\,}\mathrm{m}$ multi-mode channel has a throughput of $86.1\%$, while
the TA has throughput of $67.0\%$, including coupling into the multi-mode
fibers, but not including the detector efficiency. The $6$-state polarization
analyzer at Alice has an approximate throughput of $49.6\%$ across all six
paths, again not including the detector efficiencies. The throughput of the
PTC, TA, quantum channel, and $6$-state analyzer were measured using a laser
and power meter. Experimentally, the total combined efficiencies, including
single photon detector efficiency, are $9.3\%$ for the idler photons measured
by Alice and $1.1\%$ for the signal photons measured by Bob.

The results of our demonstration are shown in figure 2 and figure 3. Here we
describe two experimental scenarios; (i) the system and interferometers are
allowed to drift due to fluctuations from the environment, (ii) a
piezoelectric actuator in the TA is driven to produce a constant phase change
of $0.1\text{\,}\mathrm{rad}\text{\,}{\mathrm{s}}^{-1}$. Despite the lack of
active compensation of the two interferometers’ relative phase, average
visibilities of $88.5\pm 1.5$% and $86.5\pm 1.7$% are observed in the phase-
dependent superposition basis for scenarios (i) and (ii) respectively. This
corresponds to an average QBER in the superposition basis of $5.8\pm 0.8\%$
and $6.8\pm 0.9\%$ for scenarios (i) and (ii) respectively, while the QBER in
the computational basis is $4.2\pm 0.9\%$ and $3.4\%\pm 0.6\%$ for scenarios
(i) and (ii) respectively. The visibility of the superposition basis is
derived from the $C$-parameter and is a combination of the visibility of the
EPS ($95\%$), the PTC interferometer ($99\%$), and the TA interferometer
($93\%$). Combined, the estimated visibility for the superposition basis is
$87\%$ and gives a combined expected QBER of approximately $6.5\%$, which
agrees with the experimentally observed values. Other errors, such as
waveplate errors and imperfections in the beam splitters are not accounted for
but can contribute to the error in our demonstration.

With this $C$-parameter, an asymptotic key rate was calculated following the
analysis of [19]. The number of secret key bits per coincidence detection at
the asymptotic limit is estimated to be
$0.063\text{\,}\mathrm{b}\mathrm{i}\mathrm{t}\mathrm{s}\mathrm{/}\text{coincidence}$
and
$0.063\text{\,}\mathrm{b}\mathrm{i}\mathrm{t}\mathrm{s}\mathrm{/}\text{coincidence}$
((i) and (ii) respectively) compared to the maximum achievable rate for our
system of
$0.167\text{\,}\mathrm{b}\mathrm{i}\mathrm{t}\mathrm{s}\mathrm{/}\text{coincidence}$.
Despite the varying phase, there is no significant variation in QBER and
asymptotic key rate in figures 2(b), 3(b).

(a)

(b)

Figure 2: (a) Measured correlation values and phase as a function of time for
a $60\text{\,}\mathrm{s}$ link (scenario (i)). $\phi$ is the relative phase
difference between the two interferometers. (b) The corresponding asymptotic
key rate and QBER. All errors are shown as shaded regions and are approximated
by Poissonian count errors.

(a)

(b)

Figure 3: (a) Measured correlation values and phase as a function of time for
a $60\text{\,}\mathrm{s}$ demonstration while a piezo actuator is inducing a
phase change ($0.1\text{\,}\mathrm{rad}\text{\,}{\mathrm{s}}^{-1}$) in the TA
(scenario (ii)). (b) The corresponding asymptotic key rate and QBER. All
errors are shown as shaded regions and are approximated by Poissonian count
errors.

Clearly, in figures 2(a), 3(a) the $C$-parameter is shown to be phase
independent and constant throughout the entire photon exchange. This is
achieved despite deliberate phase variation between the PTC and TA (figure
3(a)). We use our Monte Carlo simulation to estimate the effects of a phase
change for our system that measures about $3000\text{\,}\text{coincidences}$
($N=3000$) in the superposition basis during a $1\text{\,}\mathrm{s}$
($t_{N}=1$) measurement time interval (figure 4). The tolerated relative phase
drift is approximately $0.512\text{\,}\mathrm{rad}\text{\,}{\mathrm{s}}^{-1}$
before a $5\%$ reduction in the asymptotic key rate. This tolerance can be
increased by having more coincidences per measurement time interval or by
reducing the measurement time, i.e. a faster measurement system. A longer
measurement time will increase the number of photon detections, but increase
the smearing effect and cause a larger drop in key rate for the same phase
change rate. In figure 4 we compare our Monte Carlo simulation to a modified
version equation (17) from [26]. The modifications are necessary to match the
definition of $C$ used in the 6-state 4-state protocol [19].

Figure 4: Change in the value of the $C$-parameter as a function of phase
variation over a fixed time period ($\phi(t_{N})$). The error bars are the
standard deviation of the Monte Carlo distribution and demonstrate the
statistical nature of the $C$-parameter that is not captured by equation (17)
of [26]. The error bars can be reduced by increasing the number of photons
collected in the time $t_{N}$. Here $i=300$ for each $\phi(t_{N})$.

Furthermore, we perform an indicative finite key rate estimation using the
approach from [26] that analyzed a 6-state entanglement-based RFI protocol.
For our analysis, we make the assumption that the correlations of the
entangled state shared between Alice and Bob are symmetric
($\braket{X_{A}\otimes X_{B}}^{2}=\braket{Y_{A}\otimes Y_{B}}^{2}$ and
$\braket{Y_{A}\otimes X_{B}}^{2}=\braket{X_{A}\otimes Y_{B}}^{2}$) such that
we can take the 6-state 4-state $C$-parameter ($C_{64}$) and relate it to the
6-state $C_{66}$ as $C_{66}=2C^{2}_{64}$. Proper investigation of this
assumption is a topic of further study. In addition, future work will be done
to use modern numerical finite size key analysis [27]. Using equations (2) and
(3) from [26], the expression for the secret-key length secure against
coherent attack for our system given by

$\displaystyle
r_{N}=\frac{n}{N}[1-I_{E}(Q^{\prime},C^{\prime})-n_{Q}fh(Q^{\prime})-\log_{2}\frac{2}{\varepsilon_{\text{EC}}}-2\log_{2}\frac{1}{\varepsilon_{\text{PA}}}-7\sqrt{\frac{\log_{2}\left(2/\bar{\varepsilon}\right)}{n}}$
$\displaystyle-\frac{30}{N}\log_{2}\left(N+1\right)],$ (3)

where $N$ is the total raw key (signals sent by Alice and detected by Bob),
$n$ is the number of detections in the key map (computational) basis,
$C^{\prime}$ ($Q^{\prime}$) is the observed $C$-parameter (error rate in the
key map basis) modified to account for finite statistics [26], and $I_{E}$ is
the upper bound on the information that Eve can gain. The other finite size
parameters are found in table 1. For scenario (i), we considered the entire
$60\text{\,}\mathrm{s}$ in determining the $C$-parameter and QBER. This gave a
total number of coincidences of $487936$, a computational basis QBER of
$4.2\%$ and a $C_{64}$-parameter of $0.8845$. With these parameters and
assuming perfect error correction, we calculate a secret key fraction of
$0.00571$, which amounts to $2785$ secret bits from the total raw key. For
(ii), considering the entire $60\text{\,}\mathrm{s}$ link would not produce
any key as the smearing of the $C$-parameter would be too great. However,
shorter integration times would not produce positive key due to low photon
statistics. For an error correction factor of $f=1.2$, no positive secret key
is achievable. However for an error correction factor of $f=1.0$, the only
positive key achieved for (ii) while using all other parameter values of table
1, is from one $13\text{\,}\mathrm{s}$ block that has a secret key fraction of
$0.00217$ or $229$ secret bits from the $120649$ coincidences of the
$13\text{\,}\mathrm{s}$ block. It is apparent that the coincidence count rates
need to be improved in order to produce a secret key that is similar with
state of the art systems reported by the literature [14, 15].

Table 1: Parameters of the finite size secret key analysis. Parameter | Description | Value
---|---|---
$\bar{\varepsilon}$ | Smooth min entropy estimation error | $2.5\times 10^{-9}$
$\varepsilon_{\text{EC}}$ | Error correction failure probability | $2.5\times 10^{-9}$
$\varepsilon_{\text{PA}}$ | Privacy amplification failure probability | $2.5\times 10^{-9}$
$\varepsilon_{\text{PE}}$ | Parameter estimation failure probability | $2.5\times 10^{-9}$
$n_{Q}$ | Fraction of computational basis used to estimate the QBER | $0.1$
$f$ | Error correction factor | 1.2

To our knowledge, our results are of the first fully passive time-bin RFI-QKD
implementation over a multi-mode fiber. We demonstrated time-bin QKD without
the need for any active phase stabilization of the sender and receiver
interferometers, and without actively correcting for spatial mode distortions.
Overall this system could therefore be considered a plug and play QKD system.
While using a hybrid entanglement-based 6-state 4-state RFI-QKD protocol, we
observed error rates of less than $7\%$ in the superposition basis, a
continuous asymptotic key rate of above
$0.06\text{\,}\mathrm{b}\mathrm{i}\mathrm{t}\mathrm{s}\mathrm{/}\text{coincidence}$,
and when accounting for finite size, achieved an estimated key length of
$2785$ secret bits over $60\text{\,}\mathrm{s}$. Our results are of great
practical importance since the ability to distribute entanglement correlations
with time-bin encoded photons can ease the deployment of quantum networks,
entanglement distribution experiments, and quantum sensing applications [28].
Furthermore, the purely passive scheme is particularly suitable for QKD with
moving platforms where the use of multi-mode interferometers are desirable
compared to single-mode interferometers due to better pointing requirements
and optical throughput. The next steps are to improve the performance of the
system in terms of photon throughput by working with a source that is
optimized for the apparatus wavelength and bandwidth, and increasing the EPS
heralding efficiency and visibility to near unity [29]. In addition, improving
the performance of the field widened interferometers with a better optical
design and custom optics to reduce the overall QBER of the system, for
instance using a compact optical layout [6]. Further theoretical studies on
the robustness of the protocol are to be conducted, particularly against rapid
phase fluctuations and finite size effects.

The authors would like to thank Dr. Norbert Lütkenhaus and Scott Johnstun for
discussions on the analysis. The authors are grateful for the support from the
Canada Foundation for Innovation, the Ontario Research Fund, the Canadian
Institute for Advanced Research, the Natural Sciences and Engineering Research
Council of Canada, Industry Canada, and the Research Centres of Excellence
program supported by the National Research Foundation (NRF) Singapore and the
Ministry of Education, Singapore. RT would like to thank the NSERC CGS-D for
personal funding.

## References

## References

  * [1] J. Brendel, N. Gisin, W. Tittel, and H. Zbinden.  Pulsed Energy-Time Entangled Twin-Photon Source for Quantum Communication.  Physical Review Letters, 82(12):2594–2597, March 1999.
  * [2] W. Tittel, J. Brendel, H. Zbinden, and N. Gisin.  Quantum Cryptography Using Entangled Photons in Energy-Time Bell States.  Physical Review Letters, 84(20):4737–4740, May 2000.
  * [3] R. T. Thew, S. Tanzilli, W. Tittel, H. Zbinden, and N. Gisin.  Experimental investigation of the robustness of partially entangled qubits over 11 km.  Physical Review A, 66(6):062304, December 2002.
  * [4] Jeongwan Jin, Jean-Philippe Bourgoin, Ramy Tannous, Sascha Agne, Christopher J Pugh, Katanya B Kuntz, Brendon L Higgins, and Thomas Jennewein.  Genuine time-bin-encoded quantum key distribution over a turbulent depolarizing free-space channel.  Optics express, 27(26):37214–37223, 2019.
  * [5] Jeongwan Jin, Sascha Agne, Jean-Philippe Bourgoin, Yanbao Zhang, Norbert Lütkenhaus, and Thomas Jennewein.  Demonstration of analyzers for multimode photonic time-bin qubits.  Phys. Rev. A, 97(4):043847, April 2018.  Publisher: American Physical Society.
  * [6] Clinton Cahall, Nurul T. Islam, Daniel J. Gauthier, and Jungsang Kim.  Multi-mode Time-delay Interferometer for Free-space Quantum Communication.  Physical Review Applied, 13(2):024047, February 2020.  arXiv: 1908.00852.
  * [7] C Gobby, aZL Yuan, and AJ Shields.  Quantum key distribution over 122 km of standard telecom fiber.  Applied Physics Letters, 84(19):3762–3764, 2004.
  * [8] ZL Yuan and AJ Shields.  Continuous operation of a one-way quantum key distribution system over installed telecom fibre.  Optics express, 13(2):660–665, 2005.
  * [9] Erik Fitzke, Lucas Bialowons, Till Dolejsky, Maximilian Tippmann, Oleg Nikiforov, Thomas Walther, Felix Wissel, and Matthias Gunkel.  Scalable network for simultaneous pairwise quantum key distribution via entanglement-based time-bin coding.  PRX Quantum, 3(2):020341, 2022.
  * [10] William S Rabinovich, Rita Mahon, Mike S Ferraro, Peter G Goetz, Mark Bashkansky, Rachel E Freeman, John Reintjes, and James L Murphy.  Free space quantum key distribution using modulating retro-reflectors.  Optics Express, 26(9):11331–11351, 2018.
  * [11] Xingyu Wang, Wei Liu, Tianyi Wu, Chang Guo, Yijun Zhang, Shanghong Zhao, and Chen Dong.  Free space measurement device independent quantum key distribution with modulating retro-reflectors under correlated turbulent channel.  Entropy, 23(10):1299, 2021.
  * [12] Anthony Laing, Valerio Scarani, John G Rarity, and Jeremy L O’Brien.  Reference-frame-independent quantum key distribution.  Physical Review A, 82(1):012304, 2010.
  * [13] Huan Chen, Jipeng Wang, Bangying Tang, Zhenhua Li, Bo Liu, and Shihai Sun.  Field demonstration of time-bin reference-frame-independent quantum key distribution via an intracity free-space link.  Opt. Lett., 45(11):3022, June 2020.
  * [14] Bao Feng, Ziyan Zhao, Shunyu Yang, Tianqi Dou, Zhenhua Li, Jipeng Wang, Zhongqi Sun, Fen Zhou, Yanxin Han, Yuqing Huang, and Haiqiang Ma.  Four-state reference-frame-independent quantum key distribution using heralded pair-coherent sources with source flaws.  Optoelectronics Letters, 17(10):636–640, October 2021.
  * [15] Bang-Ying Tang, Huan Chen, Ji-Peng Wang, Hui-Cun Yu, Lei Shi, Shi-Hai Sun, Wei Peng, Bo Liu, and Wan-Rong Yu.  Free-running long-distance reference-frame-independent quantum key distribution.  npj Quantum Information, 8(1):1–8, September 2022.  Number: 1 Publisher: Nature Publishing Group.
  * [16] Christopher J. Pugh, Sarah Kaiser, Jean-Philippe Bourgoin, Jeongwan Jin, Nigar Sultana, Sascha Agne, Elena Anisimova, Vadim Makarov, Eric Choi, Brendon L. Higgins, and Thomas Jennewein.  Airborne demonstration of a quantum key distribution receiver payload.  Quantum Sci. Technol., 2(2):024009, June 2017.  arXiv:1612.06396 [quant-ph].
  * [17] Sheng-Kai Liao, Wen-Qi Cai, Wei-Yue Liu, Liang Zhang, Yang Li, Ji-Gang Ren, Juan Yin, Qi Shen, Yuan Cao, Zheng-Ping Li, Feng-Zhi Li, Xia-Wei Chen, Li-Hua Sun, Jian-Jun Jia, Jin-Cai Wu, Xiao-Jun Jiang, Jian-Feng Wang, Yong-Mei Huang, Qiang Wang, Yi-Lin Zhou, Lei Deng, Tao Xi, Lu Ma, Tai Hu, Qiang Zhang, Yu-Ao Chen, Nai-Le Liu, Xiang-Bin Wang, Zhen-Cai Zhu, Chao-Yang Lu, Rong Shu, Cheng-Zhi Peng, Jian-Yu Wang, and Jian-Wei Pan.  Satellite-to-ground quantum key distribution.  Nature, 549(7670):43–47, September 2017.
  * [18] Hugh Podmore, I. D’Souza, J. Cain, T. Jennewein, B. L. Higgins, Y. S. Lee, A. Koujelev, D. Hudson, and A. McColgan.  QKD terminal for Canada’s Quantum Encryption and Science Satellite (QEYSSat).  In Zoran Sodnik, Bruno Cugny, and Nikos Karafolas, editors, International Conference on Space Optics — ICSO 2020, page 17, Online Only, France, June 2021. SPIE.
  * [19] Ramy Tannous, Zhangdong Ye, Jeongwan Jin, Katanya B Kuntz, Norbert Lütkenhaus, and Thomas Jennewein.  Demonstration of a 6 state-4 state reference frame independent channel for quantum key distribution.  Applied Physics Letters, 115(21):211103, 2019.
  * [20] Xiao-song Ma, Angie Qarry, Johannes Kofler, Thomas Jennewein, and Anton Zeilinger.  Experimental violation of a bell inequality with two different degrees of freedom of entangled particle pairs.  Physical Review A, 79(4):042101, 2009.
  * [21] Fumin Wang, Pei Zhang, Xiaoli Wang, and Fuli Li.  Valid conditions of the reference-frame-independent quantum key distribution.  Physical Review A, 94(6):062330, 2016.
  * [22] Juhwan Yoon, Tanumoy Pramanik, Byung-Kwon Park, Young-Wook Cho, Sang-Yun Lee, Sangin Kim, Sang-Wook Han, Sung Moon, and Yong-Su Kim.  Experimental comparison of various quantum key distribution protocols under reference frame rotation and fluctuation.  Optics Communications, 441:64–68, June 2019.
  * [23] J Robert Johansson, Paul D Nation, and Franco Nori.  Qutip: An open-source python framework for the dynamics of open quantum systems.  Computer Physics Communications, 183(8):1760–1772, 2012.
  * [24] J.R. Johansson, P.D. Nation, and Franco Nori.  Qutip 2: A python framework for the dynamics of open quantum systems.  Computer Physics Communications, 184(4):1234–1240, 2013.
  * [25] Alexander Lohrmann, Chithrabhanu Perumangatt, Aitor Villar, and Alexander Ling.  Broadband pumped polarization entangled photon-pair source in a linear beam displacement interferometer.  Applied Physics Letters, 116(2):021101, 2020.
  * [26] Lana Sheridan, Thinh Phuc Le, and Valerio Scarani.  Finite-key security against coherent attacks in quantum key distribution.  New Journal of Physics, 12(12):123019, 2010.
  * [27] Ian George, Jie Lin, and Norbert Lütkenhaus.  Numerical calculations of the finite key rate for general quantum key distribution protocols.  Physical Review Research, 3(1):013274, 2021.
  * [28] Shihan Sajeed and Thomas Jennewein.  Observing quantum coherence from photons scattered in free-space.  arXiv:2103.00298 [quant-ph], February 2021.  arXiv: 2103.00298.
  * [29] Ali Anwar, Chithrabhanu Perumangatt, Fabian Steinlechner, Thomas Jennewein, and Alexander Ling.  Entangled photon-pair sources based on three-wave mixing in bulk crystals.  Review of Scientific Instruments, 92(4):041101, 2021.