# Seeking Quality Diversity in Evolutionary Co-design of Morphology and
Control of Soft Tensegrity Modular Robots

Enrico Zardini  University of TrentoTrentoItaly 0000-0002-7475-7183
<EMAIL_ADDRESS>,  Davide Zappetti  École Polytechnique Fédérale
deLausanneSwitzerland<EMAIL_ADDRESS>,  Davide Zambrano  École
Polytechnique Fédérale deLausanneSwitzerland<EMAIL_ADDRESS>,
Giovanni Iacca  University of TrentoTrentoItaly 0000-0001-9723-1830
<EMAIL_ADDRESS>and  Dario Floreano  École Polytechnique Fédérale
deLausanneSwitzerland<EMAIL_ADDRESS>
(2021)

###### Abstract.

Designing optimal soft modular robots is difficult, due to non-trivial
interactions between morphology and controller. Evolutionary algorithms (EAs),
combined with physical simulators, represent a valid tool to overcome this
issue. In this work, we investigate algorithmic solutions to improve the
Quality Diversity of co-evolved designs of Tensegrity Soft Modular Robots
(TSMRs) for two robotic tasks, namely goal reaching and squeezing trough a
narrow passage. To this aim, we use three different EAs, i.e., MAP-Elites and
two custom algorithms: one based on Viability Evolution (ViE) and NEAT (ViE-
NEAT), the other named Double Map MAP-Elites (DM-ME) and devised to seek
diversity while co-evolving robot morphologies and neural network (NN)-based
controllers. In detail, DM-ME extends MAP-Elites in that it uses two distinct
feature maps, referring to morphologies and controllers respectively, and
integrates a mechanism to automatically define the NN-related feature
descriptor. Considering the fitness, in the goal-reaching task ViE-NEAT
outperforms MAP-Elites and results equivalent to DM-ME. Instead, when
considering diversity in terms of “illumination” of the feature space, DM-ME
outperforms the other two algorithms on both tasks, providing a richer pool of
possible robotic designs, whereas ViE-NEAT shows comparable performance to
MAP-Elites on goal reaching, although it does not exploit any map.

Soft Tensegrity Modular Robots, Co-evolution, NEAT, Viability Evolution, MAP-
Elites, Quality Diversity

††journalyear: 2021††copyright: acmcopyright††conference: 2021 Genetic and
Evolutionary Computation Conference; July 10–14, 2021; Lille,
France††booktitle: 2021 Genetic and Evolutionary Computation Conference (GECCO
’21), July 10–14, 2021, Lille, France††price: 15.00††doi:
10.1145/3449639.3459311††isbn: 978-1-4503-8350-9/21/07††ccs: Computing
methodologies Evolutionary robotics††ccs: Computing methodologies Genetic
algorithms

##  1\. Introduction

Soft robots might be one of the key technologies of the future. Indeed, their
robustness and adaptive morphology allow them to overcome situations and
perform tasks in which traditional hard robots are of limited applicability
(Mintchev and Floreano, 2016). For example, they can perform locomotion on
rough terrains with risking damage (Mintchev et al., 2018), or squeeze through
narrow passages (Cheney et al., 2015). Hence, they can be employed in the
exploration of hard-to-access environments (Hallawa et al., 2020), as well as
in the medical field (Zhang et al., 2019). Recently, researchers have applied
the soft robotic paradigm to modular robotics (Zhang et al., 2020), developing
novel soft modular robots able to explore a larger morpho-functional space
than their rigid counterparts while exhibiting robustness and mechanical
adaptability (Germann et al., 2014; Mintchev and Floreano, 2016).

Designing soft modular robots is notably difficult due to the hard-to-model
dynamics of soft materials and the non-trivial interactions between morphology
and controller (Laschi et al., 2016). These difficulties, together with the
lack of analytical methods, make evolutionary algorithms (EAs), coupled with
physics simulators, one of the most promising design tools (Howard et al.,
2019). Indeed, not only can EAs discover unconventional high-performing
solutions, but they can also provide a pool of possible designs by
”illuminating” the design space (Mouret and Clune, 2015). This latter aspect
is especially significant in soft modular robots, where the influence of the
morphology and controller parameters on the robot performance is hard to
predict, thus having a number of eligible candidates for the physical
realization is desirable. We should remark however that, although other
paradigms to design soft modular robots exist (Howison et al., 2020; Lee et
al., 2017; Vergara et al., 2017; Lee et al., 2016; Morin et al., 2014), only a
few involve the use of EAs: notably and extensively the voxel-based one (VSRs)
(Cheney et al., 2013; Shah et al., 2020) and the tensegrity soft modular
robots (TSMRs) (Zappetti et al., 2017), the latter in a preliminary work with
an open-loop controller (D. Zappetti, 2021) and in another work (Doney et al.,
2020) where the Multi-dimensional Archive of Phenotypic Elites (MAP-Elites)
(Mouret and Clune, 2015) was applied to diversify only behavioral properties
rather than morphological ones.

This work also focuses on TSMRs. More precisely, on the joint optimization of
TSMR morphologies and NN-based controllers, as the optimization of either
morphology only or controller only via EAs does not leverage the potential
synergy between them. Specifically, we perform the co-optimization using three
different approaches, namely MAP-Elites and two proposed custom EAs: 1)
Viability Evolution (ViE) (Maesani et al., 2014) coupled with Neuro-Evolution
of Augmenting Topologies (NEAT) (Stanley and Miikkulainen, 2002), denoted as
_ViE-NEAT_ , and 2) a method based on MAP-Elites, dubbed as _Double Map MAP-
Elites_ (DM-ME). The latter makes use of two maps, associated to an entity-
related and a controller-related feature descriptor (FD), respectively, and
includes a dimensionality reduction mechanism (Cully, 2019) for the automatic
definition of the controller-related FD.

We experimentally compare the three approaches on two robotic tasks: goal
reaching and squeezing through a narrow passage. Our results show that DM-ME
outperforms the other two approaches at illuminating the search space on both
tasks. In addition, although ViE-NEAT does not keep any map, it achieves very
similar performance to MAP-Elites on goal reaching (in terms of illumination).

The rest of this paper is structured as follows: Section 2 deals with the
background and related work; Section 3 introduces the proposed methods, i.e.,
ViE-NEAT and DM-ME; Section 4 describes the experimental setup; Section 5
presents the numerical results. Finally, Section 6 provides the main
conclusions of our study.

##  2\. Background

This section introduces the TSMRs employed in this work and the EAs related to
the algorithms we have considered. It also introduces the mechanisms that have
been integrated in DM-ME.

###  2.1. Tensegrity soft modular robots

The term _tensegrity_ has been coined by the architect R. Buckminster Fuller
(Fuller, 1961) and denotes a structure that maintains its mechanical integrity
through the pre-stretching of some elements constantly in tension (known as
strings or cables) connected in a network with other elements constantly under
compression (called rods or struts). Several kinds of tensegrity structures
exist: this work focuses on the icosahedron tensegrity, which consists of 6
rods connected by 24 pre-stretched cables. This particular structure has been
chosen due to its properties: the ability to deform in all directions; the
symmetry; the small number of rods and cables; the possibility to carry a
payload in the inner volume (Zappetti et al., 2017). The last point is
particularly important since it allows to equip the structure with a
servomotor, which enables structural contractions, and different kinds of
sensors.

The term _modular_ denotes robots composed of a variable number of building
blocks. In this work, the icosahedron tensegrity has been used as a module. In
detail, the nodes of the icosahedron draft eight triangular faces, which are
exploited to link the modules together by connecting the triangles’ vertices.
In the experiments, we take in consideration only linear chains of modules,
both for simplicity and feasibility of a future hardware implementation. The
appearance of a TSMR in the simulation framework that we have used (see
Section 4.1) is shown in Figure 1.

The term _soft_ refers to the fact that these robots are not rigid. Indeed,
they are able to deform their structure by contracting and expanding their
modules.

###  2.2. Neuro-evolution of augmenting topologies

NEAT (Stanley and Miikkulainen, 2002) is a well-known neuro-evolutionary
algorithm capable of evolving both topology and weights of artificial neural
networks (NNs). It exploits an advanced genetic encoding based on the
_historical marking_ mechanism. In practice, each topological feature is
marked with a number that allows to match homologous genes while performing
crossover; in this way different genotypes can be aligned and crossed over.
Moreover, the initialization of the population is addressed through a
_complexification_ process: the initial population consists of minimal
networks, which are gradually made more complex as the search progresses. In
this way the evolution will focus on the most promising parts of a very high-
dimensional search space. Diversity and innovations are preserved through
speciation: the population is split into species by means of a clustering
method that groups individuals by similarity and the competition is limited to
individuals within each species. It is also worth mentioning the mechanism
known as species-elitism, which preserves a certain number of species (the
best ones) from extinction due to stagnation, i.e., no improvement for long
time. In this work, we used the NEAT implementation provided by NEAT-
Python111https://github.com/CodeReclaimers/neat-python.

###  2.3. Diversity-driven evolutionary algorithms

Among the various EAs proposed in the literature for seeking diversity (either
explicitly or implicitly) the following two have been taken into account here:
Viability Evolution (“ViE”) (Maesani et al., 2014, 2016), which is
characterized by a good balance between exploitation and exploration, and MAP-
Elites (Mouret and Clune, 2015), which is focused on “illuminating” the search
space.

ViE is an evolutionary paradigm based on the elimination of unviable
individuals. In detail, ViE makes use of viability boundaries that
discriminate between viable individuals, which satisfy specific requirements,
and unviable ones, not satisfying them. The boundaries in question are
dynamic: at the beginning they are set so as to include all the individuals
belonging to the population; then, they are progressively shrunk either
towards the optimal values or towards target values that correspond to some
constraints. In this work we used ViE to evolve a population of TSMR
morphologies; since there are no specific constraints in this case, the
viability boundary represents a limit on the worst fitness.

MAP-Elites is a Quality Diversity (QD) algorithm, i.e., an algorithm that aims
at generating a large collection (archive) of diverse and high-performing
solutions (see (Pugh et al., 2015, 2016) for more details on QD). In
particular, it aims at discovering the best-performing solution of each cell
of a user-defined feature space, which is a lower dimensional space with
respect to the original search space. In order to be added to the archive, a
solution is first projected onto the feature space, producing a feature
descriptor; if the corresponding cell is empty or contains an individual with
worse performance, the solution is added to the archive (replacing the
occupant, if present), otherwise it is discarded. In this work we used MAP-
Elites to evolve morphology-controller pairs.

###  2.4. MAP-Elites with multiple maps

Multiple feature maps have already been employed in previous research,
although for the evolution of _a single entity_. For example, in (Pugh et al.,
2016), two MAP-Elites archives are employed for the evolution of the
controller of a wheeled robot in a maze navigation task. The two FDs used
there are the endpoint of the simulation and the main orientation of the robot
for each fifth of the simulation; thus, they are two behavioral properties of
the evolved controller. Multiple archives are used also in (Doncieux and
Coninx, 2018) (although with a different QD algorithm), but even there the FDs
represent only behavioral features. Instead, the DM-ME algorithm introduced
here is devised for the evolution of _two entities_ : a generic entity (in our
case, a TSMR morphology) and its related controller (a NN). Hence, the FDs of
the two archives refer to distinct objects.

###  2.5. Automatic feature descriptor definition

The selection of the FD is a critical factor for MAP-Elites. A mechanism to
define it in an automatic way has been introduced in (Cully, 2019). In
practice, the FD is obtained by applying a Dimensionality Reduction (DR)
algorithm, i.e., the Principal Component Analysis (PCA), on the sensory data
extracted during the evaluation of the individual; the DR’s output vector
represents the FD that is used for the insertion in the archive. Nevertheless,
the FD in question consists of continuous values, which is acceptable for an
unstructured archive, such as the one used in (Cully, 2019), but not for MAP-
Elites. Hence, an additional processing step to discretize the feature values
is required. Basically, a first discretization to determine the cells
boundaries is performed at the beginning of the evolution using the values
retrieved by the PCA for the initial population. Then, the boundaries are
recomputed every time the PCA is fitted to data or a value outside the current
boundaries is returned. It is worth also mentioning that the PCA is scale-
sensitive; hence, the sensory data are not directly provided as input to the
PCA, but is standardized first by removing the mean and dividing by the
standard deviation. As regards the standardizer, it is fitted to data every
time the fitting is done for the PCA. In conclusion, the only information that
is required is the size of the FD (i.e., the number of cells per dimension).

##  3\. Proposed Methods

We present now the algorithmic details of the proposed methods, i.e., ViE-NEAT
and Double Map MAP-Elites (DM-ME).

###  3.1. ViE-NEAT

ViE-NEAT consists in the parallel evolution of two populations, a population
of entities (morphologies) and one of NN-based controllers. The former is
evolved using ViE, the latter using NEAT. The only interaction between the two
populations happens at the evaluation time: in detail, at each generation the
individuals belonging to the two populations are randomly paired for
evaluation. In principle pairing should be one-to-one. However, since in the
NEAT implementation used here the population is not fixed in size, it may
happen that an individual from one population is paired with two individuals
from the other, due to the bigger size of the other population at play. The
pseudocode is shown in Algorithm 1.

###  3.2. Double Map MAP-Elites (DM-ME)

DM-ME deals with individuals represented by an ¡entity-NN¿ pair: the entity
encoding depends on the domain of application, whereas NEAT’s genetic encoding
is used for NNs.

Algorithm 1 ViE-NEAT

1:procedure main():

2: $P_{ViE}$ $\leftarrow$ initViEPopulation() $\triangleright$ entity pop.

3: $P_{NEAT}$ $\leftarrow$ initNEATPopulation() $\triangleright$ NN pop.

4: for $t$ = 1 $\rightarrow$ $T$ do

5: evaluate($P_{ViE}$, $P_{NEAT}$) $\triangleright$ random pairing used

6: $P_{ViE}$ $\leftarrow$ runOneGeneration($P_{ViE}$)

7: $P_{NEAT}$ $\leftarrow$ runOneGeneration($P_{NEAT}$)

Algorithm 2 Double Map Map-Elites (DM-ME)

1:procedure main($numInitialSolutions$, $archiveSize$):

2: $A_{E}$ $\leftarrow$ $\emptyset$ $\triangleright$ entity archive

3: $A_{NN}$ $\leftarrow$ $\emptyset$ $\triangleright$ NN archive

4: $PCA$ $\leftarrow$ initPCA($archiveSize$)

5: $X$ $\leftarrow$ generateRandomSolutions($numInitialSolutions$)

6: evaluate($X$) $\triangleright$ get fitness and sensory data

7: addAllToMap($A_{E}$, $X$)

8: addToNNMap($A_{NN}$, $PCA$, $X$)

9: for $t$ = 1 $\rightarrow$ $T$ do

10: $X$ $\leftarrow$ randomSelection($A_{E}$, $A_{NN}$)

11: $X$ $\leftarrow$ randomVariation($X$)

12: evaluate($X$)

13: addAllToMap($A_{E}$, $X$)

14: addToNNMap($A_{NN}$, $PCA$, $X$)

15:procedure addAllToMap($A$, $X$):

16: for $x$ $\in$ $X$ do

17: addToMap($A$, $x$)

18:procedure addToNNMap($A_{NN}$, $PCA$, $X$):

19: $SD_{X}$ $\leftarrow$ getSensoryData($X$)

20: if $PCA$.notFitted() then

21: $PCA$.fit($SD_{X}$)

22: else if $isFittingTime$ then

23: $SD_{A}$ $\leftarrow$ getSensoryData($A_{NN}$)

24: $PCA$.fit($SD_{A}$)

25: reInsert($A_{NN}$, $PCA$)

26: $PCA$.transform($X$, $SD_{X}$) $\triangleright$ compute NN FD

27: if outOfBounds($X$) then

28: $PCA$.reComputeBounds($A_{NN}$, $X$)

29: reInsert($A_{NN}$, $PCA$)

30: $PCA$.transform($X$, $SD_{X}$)

31: addAllToMap($A_{NN}$, $X$)

The functioning of DM-ME is illustrated in Algorithm 2. At the beginning, the
archives are initialized with a large number of randomly generated ¡entity-NN¿
pairs, exploiting the corresponding sensory data to fit the PCA. Then, until
the desired number of generations has been reached, a batch of elites is
sampled from the two maps, half from each of them. Next, the individuals are
mutated as follows: with equal probability, either only the entity is mutated,
only the NN is mutated, or both of them are mutated. The mutation of the
entity is domain-dependent, whereas the NNs are mutated according to NEAT’s
mutation operator (crossover is not used). After mutation, the individuals are
evaluated and added to the archives based on the respective FDs. In this
regard, each of the two archives is an independent MAP-Elites map, and the
addToMap() archive insertion procedure is handled accordingly: an individual
(an ¡entity-NN¿ pair) is added to an archive if the corresponding cell is
empty or occupied by a worse individual, hence, it may be added to one map,
but not to the other, depending on the fitness of the currently stored
individuals. The peculiarity lies in the fact that a FD related to the entity
is used for the first map, whereas a FD related to the NN is employed for the
second one. As shown in (Pugh et al., 2016), the use of multiple adequate FDs
(in the form of multiple maps for MAP-Elites) leads to better results,
especially in difficult domains, since different perspectives are taken into
account. As for our DM-ME, the two FDs employed are obtained in different
ways. In particular, the entity is projected as usual on a user-defined
feature space, in order to obtain the corresponding FD. Instead, for the NN,
the mechanism described in Section 2.5 is used. More specifically, the PCA is
fitted using the sensory data of the individuals contained in the NN archive
with a frequency that decreases exponentially, as in (Cully, 2019); obviously,
at every fitting, it is necessary to re-compute the FDs of the individuals
that are already in the archive and re-insert them, which is done also when
the space is discretized between two PCA updates due to the presence of values
outside the current cells boundaries.

##  4\. Experimental setup

This section presents the details of our experimental setup.

###  4.1. Implementation details

Apart from NEAT, for which we used the existing NEAT-Python library, the other
EAs involved in the experimentation were implemented from scratch in Python,
and coupled with a custom TSMR simulation framework developed in C++, named
Tensoft. The latter is based on the NASA Tensegrity Robotics
Toolkit222https://github.com/NASA-Tensegrity-Robotics-Toolkit/NTRTsim, a
collection of tools for modeling and simulating tensegrity robots, built on
the top of the Bullet
Physics333https://github.com/bulletphysics/bullet3/archive/2.88.tar.gz engine.
Our source code is available at: https://github.com/lis-epfl/Tensoft-G21.

###  4.2. Tasks

TSMRs have only recently become a subject of study, as demonstrated by (Lee et
al., 2020); hence, apart for locomotion (D. Zappetti, 2021), their
capabilities are still largely unknown. Here, we investigate two tasks that
are fairly more complex than plain locomotion: _goal reaching_ and
_squeezing_.

Goal reaching consists in having the robot moving towards a target placed
somewhere in the environment. While this is a typical task in robotics
research, modular soft structures have been barely taken into account in this
regard. In practice, the objective consists in finding a controller able to
drive the associated morphology towards different targets, based on the
sensory information provided, i.e., distance and bearing to the target. In
detail, distance and bearing are computed with respect to the position and the
direction of the front face of the head module, respectively. An example of
final situation of the task in Tensoft is shown in Figure 1 (left).

Table 1. Main parameters for all the experiments. | Goal Reaching | Squeezing
---|---|---
# Runs | 10 | 10
Run’s budget | 45000 evaluations | 45000 evaluations
# Targets | 4 | 2
Distance | $45\text{\,}\mathrm{c}\mathrm{m}$ | $60\text{\,}\mathrm{c}\mathrm{m}$
Bearing | 90° l, 45° l, 45° r, 90° r | 5° l, 5° r
# Simulation seeds | 2 | 2
Simulation time | $40\text{\,}\mathrm{s}$ | $40\text{\,}\mathrm{s}$
Table 2. Parameter setting employed for ViE-NEAT.

ViE
---
Population size (init.) | 48
# Mutants | 48

NEAT
---
Population size | 54
Individual elitism | 3
Species elitism | 2
Compatibility threshold | 2.85

Table 3. Parameter setting for MAP-Elites and DM-ME. | MAP-Elites | DM-ME
---|---|---
# Initial solutions | 1080 | 1080
Batch size | 24 | 24
Archive(s) size | [9, 10, 9, 10] | [9, 10] [9,10]
PCA updates | [0, 50, 150, 350, 750, 1550]
Trajectory sampling | $1\text{\,}\mathrm{s}$ | $1\text{\,}\mathrm{s}$

Squeezing consists in having the robot shrinking through a restricted aperture
that is narrower than the robot size in terms of width and/or height. Due to
the high degree of complexity, this task has not been addressed extensively in
previous research. In this work, squeezing has been dealt with as an advanced
form of goal reaching, by positioning the target beyond an aperture narrower
than the maximum robot width. At the beginning of the simulation, the robot is
enclosed by four walls, and it has to pass through the aperture to reach the
target. In this case the sensory information provided includes distance and
bearing to the target, distance and bearing with respect to the center of the
entrance of the aperture (computed analogously to the ones referred to the
target), and the presence of obstacles in front of the head module within a
certain range (set to $10\text{\,}\mathrm{c}\mathrm{m}$ in our experiments).
An example of final situation of the task in Tensoft is shown in Figure 1
(right).

Figure 1. Example of final situation for goal reaching (left) and squeezing
(right). The cyan cylinder depicts the target.

###  4.3. Encoding

The TSMR morphology encoding includes global properties, affecting the entire
robot, and local ones, specific to each module. In detail, the global genes
are the number of modules and the stiffness, which determines the degree of
deformability of the modules. Instead, the local ones include only the number
$\rho$, which identifies the triangular face that is connected to the next
module in the chain. This influences not only the weight balance, but also how
the robot moves since the servomotor acts always on the same pair of faces.
The allowed mutations include: the addition of a new module in a random
position; the deletion of a random module in the chain; the change of the
robot stiffness; the change of the connection face to the next module (local
mutation).

As concerns the controllers, the encoding and the allowed mutations are
defined by NEAT. In particular, the TSMR controller has been implemented as a
feed-forward NN that at each timestep of the simulation takes as input the
sensory information provided for the considered task and produces as output
the actuation parameters, i.e., the frequency $f$ and phase $\phi$ of a
sinusoidal signal that controls contraction/expansion for all modules.

###  4.4. Algorithm and task configurations

Three different approaches have been applied to both the goal reaching and
squeezing tasks:

  * •

co-evolution of morphology, evolved through ViE, and controller, evolved
through NEAT (i.e., using ViE-NEAT);

  * •

evolution of morphology-controller pairs through MAP-Elites;

  * •

evolution of morphology-controller pairs using DM-ME.

The main parameters for all the experiments, common to the three approaches,
are reported for each task in Table 1. For a fair comparison, the budget,
limited by the computational resources available, has been defined for all
algorithms in terms of total number of solution evaluations.

Each solution is evaluated against different targets placed at different
orientations and, since the simulations are noisy, for each target multiple
simulations are run, with different seeds for the actuation noise. All the
simulations contribute to the fitness of the individual, to be minimized,
which is computed for both tasks as:

$\displaystyle\overline{d_{t}}=\frac{1}{N_{s}}\sum_{s=1}^{N_{s}}d_{ts}\quad
t\in\\{1,...\ ,N_{t}\\}$ (2)
$\displaystyle\begin{multlined}f=\frac{1}{N_{t}}\sum_{t=1}^{N_{t}}\overline{d_{t}}+\frac{1}{2}\times(\max_{t=1\dots
N_{t}}\overline{d_{t}}-\frac{1}{N_{t}}\sum_{t=1}^{N_{t}}\overline{d_{t}})\end{multlined}f=\frac{1}{N_{t}}\sum_{t=1}^{N_{t}}\overline{d_{t}}+\frac{1}{2}\times(\max_{t=1\dots
N_{t}}\overline{d_{t}}-\frac{1}{N_{t}}\sum_{t=1}^{N_{t}}\overline{d_{t}})$

where $N_{s}$ is the number of simulation seeds, $N_{t}$ is the number of
targets, $d_{ts}$ is the final distance of the individual from the $t$-th
target using the $s$-th seed, and $f$ is the resulting fitness.

In practice, first the average final distance across simulation seeds is
calculated for each target. Then, the fitness of the individual is computed as
the mean value of the average distances plus a penalty equal to half of the
difference between the worst (maximum) and the mean average distance across
targets. Actually, in the squeezing experiments $d_{ts}$ represents the final
distance with any potential bonus already deducted: if the robot succeeds in
crossing the entrance of the aperture by at least
$4\text{\,}\mathrm{c}\mathrm{m}$, a fixed bonus of
$4\text{\,}\mathrm{c}\mathrm{m}$ is deducted from the final distance. In this
way, the individuals capable of entering the aperture are favored by
evolution.

As concerns the squeezing experiments, we set the wall’s width to
$16.67\text{\,}\mathrm{c}\mathrm{m}$, its height to
$12\text{\,}\mathrm{c}\mathrm{m}$, and the aperture width to
$8\text{\,}\mathrm{c}\mathrm{m}$. In particular, the wall’s width, which
corresponds to the length of the aperture, is set to 1.5 times the length of a
TSMR module. Instead, the width of the aperture is set to 72% of a module’s
width (note that the width of a module is equal to its length, i.e., about
$11\text{\,}\mathrm{c}\mathrm{m}$).

A brief description of the approach-specific configurations, which are common
to the two tasks, is provided in the following paragraphs. Indeed, only the
number of controller input nodes is task-dependent (2 for goal reaching, 5 for
squeezing, see Section 4.2).

#### ViE-NEAT

The parameter configuration of the ViE-NEAT approach is shown in Table 2. This
configuration is characterized by elitism at both individual and species
level, a one-to-one (on average) morphology-controller pairing for evaluation
(a slightly larger NEAT population is required to compensate for the presence
of elite individuals), and a relatively low threshold for speciation.

If an individual from one population is paired with two individuals from the
other (due to the bigger size of the other population at play), the fitness of
that individual is computed as:

(3) $f=\frac{1}{N_{p}}\sum_{p=1}^{N_{p}}f_{p}+\frac{1}{2}\times(\max_{p=1\dots
N_{p}}f_{p}-\frac{1}{N_{p}}\sum_{p=1}^{N_{p}}f_{p})$

where $N_{p}$ is the number of paired individuals belonging to the other
population, $f_{p}$ is the fitness obtained with the $p$-th paired individual,
according to Eq. (2), and $f$ is the resulting fitness. Since the controllers
are not associated to a fixed morphology, the number of output nodes, which is
morphology-dependent, is set to two times the maximum allowed number of
modules, i.e., $20$ $(2\times 10)$, and only the required number of outputs is
considered at the evaluation time (the same holds for the MAP-Elites and DM-ME
experiments).

#### MAP-Elites and DM-ME

The configuration used for each of the two algorithms is reported in Table 3;
the only difference lies in the archive(s) size. The mutation of individuals
(morphology-controller pairs) in MAP-Elites is performed as in DM-ME (see
Section 3).

As regards the feature descriptors, the FD used for MAP-Elites corresponds to
the concatenation of the FDs employed in the two archives of DM-ME, i.e., the
morphology-related and the controller-related FDs (see the archives size in
Table 3). As a consequence, the archive filled by MAP-Elites has a higher
number of cells than the ones explored by DM-ME (8100 cells vs 90 cells
$\times$ 2 archives). Nevertheless, a comparison with equal number of cells
would require using different FDs for the two algorithms, which would make the
comparison unfair. Moreover, DM-ME has been specifically designed to reduce
the total number of cells and thus to increase the selective pressure (with
equal FDs) w.r.t. the original MAP-Elites.

In detail, the morphology-related FD is represented by the number of modules
and the module stiffness; instead, the controller-related FD is provided by
the PCA, exploiting the trajectory of the head module as sensory data. Since
each morphology-controller pair is simulated with multiple targets, the
sensory data vector are defined as the concatenation of the various
trajectories: 320 elements in the case of goal reaching (4 targets $\times$
$40\text{\,}\mathrm{s}$ $\times$ 2 coordinates), and 160 elements in the case
of squeezing (2 targets instead of 4). In addition, since a solution is
simulated multiple times with different seeds, the average trajectory across
seeds is taken as representative for each target. Actually, the size of the
controller-related FD has been chosen so that the number of cells is the same
for the morphology-related and the controller-related feature space.

##  5\. Results

In the following, we present the results achieved in the goal reaching and
squeezing experiments. We must remark that the fitness trends, the heatmaps
and the statistics shown here have been computed after applying the following
transformation to the fitness values obtained throughout the evolution:

(4) $f_{s}=d_{init}-f$

where $d_{init}$ is the initial distance from the target
($45\text{\,}\mathrm{c}\mathrm{m}$ for goal reaching and
$60\text{\,}\mathrm{c}\mathrm{m}$ for squeezing, see Table 1), $f$ is the
fitness of the individual measured according to Eq. (2) or (3) during the
evolutionary process, and $f_{s}$ is the fitness used for plots and statistics
computation. In this way, the minimization problem is presented as a
maximization one, to facilitate the analysis in terms of QD.

Figure 2. Best fitness (mean $\pm$ std. dev. across $10$ runs) over
evaluations for the goal reaching task. In particular, the trends shown are
related to: the morphology population for ViE-NEAT; the one and only archive
for MAP-Elites (ME); the morphology archive for DM-ME. The fitness values (to
be maximized) shown here have been obtained through Eq. (4). Figure 3. Best
fitness (mean $\pm$ std. dev. across $10$ runs) over evaluations for the
squeezing task. In particular, the trends shown are related to: the morphology
population for VIE-NEAT; the one and only archive for MAP-Elites (ME); the
morphology archive for DM-ME. The fitness values (to be maximized) shown here
have been obtained through Eq. (4).

All the experiments have been executed using a High Performance Computing
(HPC) infrastructure. In detail, for each run, one node with 2 Intel Broadwell
processors running at 2.6 GHz, with 14 cores each (28 cores in total), and
16/32 GB of RAM has been used.

Table 4. Descriptive statistics of the best fitness across $10$ runs for the
goal reaching task (up); Wilcoxon rank-sum test ($\alpha\,{=}\,0.05$) applied
to the distributions of best fitness (bottom).

Algorithm | Mean | Std | Min | Max | Median
---|---|---|---|---|---
ViE-NEAT | 26.21 | 3.22 | 21.76 | 30.79 | 26.08
MAP-Elites | 22.16 | 4.13 | 15.40 | 29.14 | 21.50
DM-ME | 25.76 | 4.93 | 20.66 | 39.36 | 24.22

Algorithms | p-value
---|---
ViE-NEAT | MAP-Elites | 0.034293
ViE-NEAT | DM-ME | 0.545350
MAP-Elites | DM-ME | 0.082099

Table 5. Descriptive statistics of the best fitness across $10$ runs for the
squeezing task (up); Wilcoxon rank-sum test ($\alpha\,{=}\,0.05$) applied to
the distributions of best fitness (bottom).

Algorithm | Mean | Std | Min | Max | Median
---|---|---|---|---|---
ViE-NEAT | 24.88 | 5.63 | 18.81 | 35.94 | 22.75
MAP-Elites | 22.63 | 2.74 | 18.60 | 26.05 | 22.63
DM-ME | 24.16 | 4.07 | 18.99 | 31.00 | 23.49

Algorithms | p-value
---|---
ViE-NEAT | MAP-Elites | 0.496292
ViE-NEAT | DM-ME | 0.820596
MAP-Elites | DM-ME | 0.496292

Table 6. Descriptive statistics of the QD score of the morphology archives
generated in the goal reaching (up) and squeezing (bottom) experiments.

Algorithm | Mean | Std | Min | Max | Median
---|---|---|---|---|---
ViE-NEAT | 535.61 | 59.21 | 461.81 | 641.67 | 524.60
MAP-Elites | 552.27 | 35.41 | 519.24 | 642.73 | 544.26
DM-ME | 656.08 | 35.05 | 589.36 | 707.55 | 655.60

Algorithm | Mean | Std | Min | Max | Median
---|---|---|---|---|---
ViE-NEAT | 1196.26 | 60.05 | 1099.57 | 1291.28 | 1203.39
MAP-Elites | 1293.25 | 15.30 | 1264.05 | 1314.11 | 1296.67
DM-ME | 1322.52 | 14.60 | 1296.39 | 1341.12 | 1324.52

Table 7. Wilcoxon rank-sum test ($\alpha\,{=}\,0.05$) applied to the QD score of the morphology archives generated in the goal reaching (GR) and squeezing (SQ) experiments. Algorithms | p-value (GR) | p-value (SQ)
---|---|---
ViE-NEAT | MAP-Elites | 0.325751 | 0.000881
ViE-NEAT | DM-ME | 0.000881 | 0.000157
MAP-Elites | DM-ME | 0.000507 | 0.002497

#### Goal reaching

Figure 2 shows the best fitness trend (mean $\pm$ std. dev. across $10$ runs)
for ViE-NEAT, MAP-Elites and DM-ME. In detail, ViE-NEAT and DM-ME exhibit an
almost equivalent performance on average, whereas MAP-Elites turns out to be
worse. The descriptive statistics of the best fitness across runs are reported
in Table 4 (up) for each algorithm. It can be seen that both ViE-NEAT and DM-
ME obtain a higher mean best fitness value than MAP-Elites. It is also worth
highlighting that the best individual across all the goal reaching experiments
has been discovered by DM-ME444The videos of the best individuals are
available at https://tinyurl.com/ynaav7ek.. Both the equivalence of ViE-NEAT
and DM-ME and the superiority of ViE-NEAT with respect to MAP-Elites are
statistically significant, see the results of the Wilcoxon rank-sum test in
Table 4 (bottom). The table also shows that in this case DM-ME does not
statistically outperform MAP-Elites.

The other evaluation metric we are interested in evaluating here is the
illumination of the feature space(s). With reference to Figure 4 (left column,
first three rows), which provides a heatmap representation of the archives, or
archive projections, generated in five runs of the goal reaching experiments,
it can be noted that DM-ME appears to perform better in terms of morphology
archives produced (although the difference is not so pronounced), whereas ViE-
NEAT shows almost equivalent performance to MAP-Elites. As regards the
controller archives, it is worth remembering that each cell corresponds to a
different trajectory according to the PCA and, since not all trajectories can
lead to good results, it is reasonable to find also bad performing individuals
inside them. In addition, the coverage of the controller archives is
influenced by the re-discretization of the features space. Hence, it is
reasonable to focus on the morphology archives. The difference in terms of
illumination capability between the three algorithms is confirmed by the
statistics computed on the QD score (i.e, the sum of the fitness of the
individuals in the archive (Pugh et al., 2015)) of the morphology archives.
All the ten runs per algorithm have been taken into account for the
computation of these statistics, see Table 6 (up). Both the superiority of DM-
ME on the other two methods and the equivalence between ViE-NEAT and MAP-
Elites in terms of QD score are statistically significant, see the results of
the Wilcoxon rank-sum test in Table 7.

(a) ViE-NEAT, morphology archive (projection), goal reaching

(b) ViE-NEAT, morphology archive (projection), squeezing

(c) MAP-Elites, morphology archive (projection), goal reaching

(d) MAP-Elites, morphology archive (projection), squeezing

(e) DM-ME, morphology archive, goal reaching

(f) DM-ME, morphology archive, squeezing

(g) MAP-Elites, controller archive (projection), goal reaching

(h) MAP-Elites, controller archive (projection), squeezing

(i) DM-ME, controller archive, goal reaching

(j) DM-ME, controller archive, squeezing

Figure 4. Archives generated by five runs of each experiment configuration.
The first row shows the archives generated by ViE-NEAT, obtained by projecting
on the morphology feature space the individuals discovered throughout the
evolution (the reconstruction of the controller archives would require
modifying the algorithm since ViE-NEAT does not include the PCA). The second
and third row show the morphology-related archives generated by MAP-Elites
(i.e., the projection from the resulting single archive) and DM-ME,
respectively. The fourth and fifth row display the corresponding controller-
related archives. In all cases, the fitness values (to be maximized) shown
here have been obtained through Eq. (4).

An inverse analysis has been also performed on MAP-Elites and DM-ME. The
double archives generated by DM-ME have been projected onto single ones, and
then compared with the single archives produced by MAP-Elites. This requires
to re-simulate the individuals contained in the morphology archives of DM-ME
in order to collect their trajectories and compute the controller-related FDs.
Specifically, the PCA and the cell boundaries used for the projection (i.e.,
the final ones) are potentially different from the ones used for the
evaluation of the considered individuals during the evolution. Hence, if the
reconstructed single archives were projected back onto separate morphology-
controller ones, the controller archives would be potentially different from
the original ones. The results of the double-to-single projection show that
most of the individuals contained in the controller archives generated by DM-
ME tend to be concentrated in a small number of cells in the morphological
feature space, typically with a medium-high number of modules and very low
stiffness. This is reasonable, since those morphological properties allow a
greater variety of mid-performing behaviors. Instead, as expected, the single
archives generated by MAP-Elites show a higher coverage of the 4-dimensional
feature space, especially in the region of medium-low stiffness. The heatmap
representation of the single archives can be found in the Supplementary
Material.

Concerning the execution time, MAP-Elites ($154.3\pm 9.22$ hours) and DM-ME
($145.4\pm 10.76$ hours) turn out to be considerably more expensive than ViE-
NEAT ($94.3\pm 13.15$ hours). It is worth mentioning however that both MAP-
Elites and DM-ME present an overhead due to the need for collecting the
sensory data (trajectories), which is written by the simulator into temporary
files.

#### Squeezing

Figure 3 shows the best fitness trend (mean $\pm$ std. dev. across $10$ runs)
for the three algorithms. The black dashed line represents the fitness limit
corresponding to a robot entering the aperture in a single simulation. By
looking at the plot, it turns out that the three algorithms discover well-
performing solutions. In this case, the best individual has been discovered by
ViE-NEAT, and again both ViE-NEAT and DM-ME obtain a higher mean best fitness
value than MAP-Elites, see the descriptive statistics reported in Table 5
(up). On the other hand, the three algorithms result statistically equivalent,
see the results of the Wilcoxon rank-sum test provided in Table 5 (bottom).
Nevertheless, by looking at the behavior of the best individuals555The videos
of the best individuals are available at https://tinyurl.com/y3e9neej., it
turns out that the best individual evolved by DM-ME is the only one, among the
best individuals obtained by all algorithms, that is able to pass through the
aperture and reach the target if the simulation is allowed to continue beyond
the $40\text{\,}\mathrm{s}$ used in the evolutionary process. This also
demonstrates that the task is achievable.

As concerns the illumination of the search space, Figure 4 (right column)
provides a heatmap representation of the archives, or archive projections,
generated in five runs of the squeezing experiments. As in the case of goal
reaching, we focus the analysis on the morphology archives. In this case, some
runs of ViE-NEAT are unable to fill all the cells of the grid. Except for
that, the archives generated appear to be very similar, especially the ones
produced by MAP-Elites and DM-ME. Nevertheless, by looking at the statistics
of the QD score of the archives generated in the ten runs, reported in Table 6
(bottom), it turns out that DM-ME outperforms MAP-Elites, which in turn
outperforms ViE-NEAT. Table 7 provides the statistical evidence based on the
Wilcoxon rank-sum test.

The heatmap representation of the single archives, the ones generated by MAP-
Elites and the ones resulting from the projection of the double archives
produced by DM-ME, can be found in the Supplementary Material. The
observations made for goal reaching, as concerns this inverse analysis, hold
also for squeezing.

As in goal reaching, the execution time of MAP-Elites ($103.3\pm 5.87$ hours)
and DM-ME ($94.9\pm 12.23$ hours) turns out to be higher than that of ViE-NEAT
($70.9\pm 6.07$ hours). The lower execution time with respect to goal reaching
is mainly due to the number of targets considered for each individual, which
in this case is 2 instead of 4 (see Table 1). However, a single simulation of
the squeezing task tends to require a little more time than one of the goal
reaching task, due to the interactions between the robot and the walls.

##  6\. Conclusions

In this work we have addressed the joint optimization of morphology and
controller of TSMRs. In detail, we have considered three different
evolutionary approaches, i.e., MAP-Elites, ViE-NEAT, and DM-ME, with the last
two being algorithms proposed here for co-evolving morphologies and
controllers. In order to compare the three algorithms, we have conducted an
experimental campaign on two robotic tasks: goal reaching and squeezing. As
concerns goal reaching, ViE-NEAT outperforms MAP-Elites and results equivalent
to DM-ME in terms of best fitness. As regards squeezing, the three algorithms
achieve similar results in terms of quality of discovered solutions. Moreover,
DM-ME outperforms the other two approaches in terms of
illumination/exploration of the feature space (measured with QD score) in both
tasks. The higher total fitness of the resulting archive implies a higher
number of well-performing candidates for the physical realization. Notably,
ViE-NEAT achieves very similar performance to MAP-Elites in terms of
illumination of the search space in the goal reaching task, although it does
not exploit any map. Future work includes testing other configurations of DM-
ME and co-evolution, other robotic tasks, and more complex TSMR structures.

###### Acknowledgements.

We thank Daniele Bissoli for implementing the original framework on which this
work was built.

## References

  * (1)
  * Cheney et al. (2015) Nick Cheney, Josh Bongard, and Hod Lipson. 2015\.  Evolving Soft Robots in Tight Spaces. In _Proceedings of the Genetic and Evolutionary Computation Conference_. ACM, New York, NY, USA, 935–942.
  * Cheney et al. (2013) Nick Cheney, Robert MacCurdy, Jeff Clune, and Hod Lipson. 2013\.  Unshackling evolution: evolving soft robots with multiple materials and a powerful generative encoding. In _Proceedings of the Genetic and Evolutionary Computation Conference_. ACM, New York, NY, USA, 167–174.
  * Cully (2019) Antoine Cully. 2019\.  Autonomous skill discovery with quality-diversity and unsupervised descriptors. In _Proceedings of the Genetic and Evolutionary Computation Conference_. ACM, New York, NY, USA, 81–89.
  * D. Zappetti (2021) D. Floreano D. Zappetti, J. M. Bejjani. 2021\.  Evolutionary co-design of morphology and control of soft tensegrity modular robots with programmable stiffness.  arXiv:2101.11772 [cs.RO]
  * Doncieux and Coninx (2018) Stephane Doncieux and Alexandre Coninx. 2018.  Open-Ended Evolution with Multi-Containers QD. In _Proceedings of the Genetic and Evolutionary Computation Conference Companion_ (Kyoto, Japan) _(GECCO ’18)_. Association for Computing Machinery, New York, NY, USA, 107–108.  https://doi.org/10.1145/3205651.3205705
  * Doney et al. (2020) K. Doney, A. Petridou, J. Karaul, A. Khan, G. Liu, and J. Rieffel. 2020\.  Behavioral Repertoires for Soft Tensegrity Robots. In _2020 IEEE Symposium Series on Computational Intelligence (SSCI)_. IEEE, Canberra, Australia, 2265–2271.  https://doi.org/10.1109/SSCI47803.2020.9308218
  * Fuller (1961) Buckminster Fuller. 1961\.  Tensegrity.  _Portfolio Artnews Annual_ 4 (1961), 112–127.
  * Germann et al. (2014) Jürg Germann, Andrea Maesani, Ramon Pericet-Camara, and Dario Floreano. 2014.  Soft cells for programmable self-assembly of robotic modules.  _Soft Robotics_ 1, 4 (2014), 239–245.
  * Hallawa et al. (2020) Ahmed Hallawa, Giovanni Iacca, Cagatay Sariman, Touhidur Rahman, Michael Cochez, and Gerd Ascheid. 2020\.  Morphological evolution for pipe inspection using Robot Operating System (ROS).  _Materials and Manufacturing Processes_ 35, 6 (2020), 714–724.
  * Howard et al. (2019) David Howard, Agoston E Eiben, Danielle Frances Kennedy, Jean-Baptiste Mouret, Philip Valencia, and Dave Winkler. 2019.  Evolving embodied intelligence from materials to machines.  _Nature Machine Intelligence_ 1, 1 (2019), 12–19.
  * Howison et al. (2020) Toby Howison, Simon Hauser, Josie Hughes, and Fumiya Iida. 2020.  Reality-assisted evolution of soft robots through large-scale physical experimentation: a review.  arXiv:2009.13960 [cs.RO]
  * Laschi et al. (2016) Cecilia Laschi, Barbara Mazzolai, and Matteo Cianchetti. 2016\.  Soft robotics: Technologies and systems pushing the boundaries of robot abilities.  _Science Robotics_ 1, 1 (2016), 11 pages.
  * Lee et al. (2017) Chiwon Lee, Myungjoon Kim, Yoon Jae Kim, Nhayoung Hong, Seungwan Ryu, H Jin Kim, and Sungwan Kim. 2017.  Soft robot review.  _International Journal of Control, Automation and Systems_ 15, 1 (2017), 3–15.
  * Lee et al. (2020) Hajun Lee, Yeonwoo Jang, Jun Kyu Choe, Suwoo Lee, Hyeonseo Song, Jin Pyo Lee, Nasreena Lone, and Jiyun Kim. 2020\.  3D-printed programmable tensegrity for soft robotics.  _Science Robotics_ 5, 45 (2020), 11 pages.
  * Lee et al. (2016) Jun-Young Lee, Woong-Bae Kim, Woo-Young Choi, and Kyu-Jin Cho. 2016\.  Soft robotic blocks: Introducing SoBL, a fast-build modularized design block.  _IEEE Robotics & Automation Magazine_ 23, 3 (2016), 30–41.
  * Maesani et al. (2014) Andrea Maesani, P. N. R. Fernando, and Dario Floreano. 2014\.  Artificial Evolution by Viability Rather than Competition.  _PLoS ONE_ 9 (2014), 1–12.
  * Maesani et al. (2016) Andrea Maesani, Giovanni Iacca, and Dario Floreano. 2016\.  Memetic Viability Evolution for Constrained Optimization.  _IEEE Transactions on Evolutionary Computation_ 20, 1 (2016), 125–144.
  * Mintchev and Floreano (2016) Stefano Mintchev and Dario Floreano. 2016.  Adaptive morphology: A design principle for multimodal and multifunctional robots.  _IEEE Robotics & Automation Magazine_ 23, 3 (2016), 42–54.
  * Mintchev et al. (2018) Stefano Mintchev, Davide Zappetti, Jerome Willemin, and Dario Floreano. 2018.  A soft robot for random exploration of terrestrial environments. In _2018 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, Piscataway, NJ, USA, 7492–7497.
  * Morin et al. (2014) Stephen A Morin, Sen Wai Kwok, Joshua Lessing, Jason Ting, Robert F Shepherd, Adam A Stokes, and George M Whitesides. 2014.  Elastomeric tiles for the fabrication of inflatable structures.  _Advanced Functional Materials_ 24, 35 (2014), 5541–5549.
  * Mouret and Clune (2015) Jean-Baptiste Mouret and Jeff Clune. 2015.  Illuminating search spaces by mapping elites.  arXiv:1504.04909 [cs.AI]
  * Pugh et al. (2015) Justin Pugh, L. Soros, Paul Szerlip, and Kenneth Stanley. 2015\.  Confronting the Challenge of Quality Diversity. In _Proceedings of the Genetic and Evolutionary Computation Conference_. ACM, New York, NY, USA, 967–974.
  * Pugh et al. (2016) Justin K. Pugh, Lisa B. Soros, and Kenneth O. Stanley. 2016\.  Quality Diversity: A New Frontier for Evolutionary Computation.  _Frontiers in Robotics and AI_ 3 (2016), 17 pages.
  * Shah et al. (2020) Dylan Shah, Bilige Yang, Sam Kriegman, Michael Levin, Josh Bongard, and Rebecca Kramer-Bottiglio. 2020.  Shape Changing Robots: Bioinspiration, Simulation, and Physical Realization.  _Advanced Materials_ 2002882 (2020), 12 pages.
  * Stanley and Miikkulainen (2002) Kenneth O. Stanley and Risto Miikkulainen. 2002.  Evolving Neural Networks through Augmenting Topologies.  _Evolutionary Computation_ 10 (2002), 99–127.
  * Vergara et al. (2017) Andrea Vergara, Yi-sheng Lau, Ricardo-Franco Mendoza-Garcia, and Juan Cristóbal Zagal. 2017.  Soft modular robotic cubes: toward replicating morphogenetic movements of the embryo.  _PloS one_ 12, 1 (2017), e0169179.
  * Zappetti et al. (2017) Davide Zappetti, Stefano Mintchev, Jun Shintake, and Dario Floreano. 2017.  Bio-inspired Tensegrity Soft Modular Robots. In _Proceedings of Conference Living Machines: Biomimetic and Biohybrid Systems_. Springer, Cham, Switzerland, 497–508.
  * Zhang et al. (2019) Boyu Zhang, Yingwei Fan, Penghui Yang, Tianle Cao, and Hongen Liao. 2019.  Worm-like soft robot for complicated tubular environments.  _Soft robotics_ 6, 3 (2019), 399–413.
  * Zhang et al. (2020) Chao Zhang, Pingan Zhu, Yangqiao Lin, Zhongdong Jiao, and Jun Zou. 2020.  Modular Soft Robotics: Modular Units, Connection Mechanisms, and Applications.  _Advanced Intelligent Systems_ 2, 6 (2020), 1900166\.

## Supplementary Material

MAP-Elites

Double Map MAP-Elites (DM-ME)

Figure 5. Single archives for the goal reaching task. In detail, the first two
rows show the single archives generated by the 10 runs of MAP-Elites, whereas
the other two show the projection of the double archives produced by Double
Map MAP-Elites (DM-ME) onto a single one.

MAP-Elites

Double Map MAP-Elites (DM-ME)

Figure 6. Single archives for the squeezing task. In detail, the first two
rows show the single archives generated by the 10 runs of MAP-Elites, whereas
the other two show the projection of the double archives produced by Double
Map MAP-Elites (DM-ME) onto a single one.


# Charge-$4e$ superconductivity from nematic superconductors in 2D and 3D

Shao-Kai Jian  Condensed Matter Theory Center, Department of Physics,
University of Maryland, College Park, Maryland 20742, USA  Yingyi Huang
Institute for Advanced Study, Tsinghua University, Beijing 100084, China  Hong
Yao<EMAIL_ADDRESS>Institute for Advanced Study, Tsinghua
University, Beijing 100084, China  State Key Laboratory of Low Dimensional
Quantum Physics, Tsinghua University, Beijing 100084, China

###### Abstract

Charge-$4e$ superconductivity as a novel phase of matter remains elusive so
far. Here we show that charge-4e phase can arise as a vestigial order above
the nematic superconducting transition temperature in time-reversal-invariant
nematic superconductors. On the one hand, the nontrivial topological
defect—nematic vortex—is energetically favored over the superconducting phase
vortex when the nematic stiffness is less than the superfluid stiffness;
consequently the charge-$4e$ phase emerges by proliferation of nematic
vortices upon increasing temperatures. On the other hand, the Ginzburg-Landau
theory of the nematic superconductors has two distinct decoupling channels to
either charge-$4e$ orders or nematic orders; by analyzing the competition
between the effective mass of the charge-$4e$ order and the cubic potential of
the nematic order, we find a sizable regime where the charge-$4e$ order is
favored. These two analysis consistently show that nematic superconductors can
provide a promising route to realize charge-$4e$ phases, which may apply to
candidate nematic superconductors such as PbTaSe2 and twisted bilayer
graphene.

Introduction.—Featuring the condensation of quartets with four times the
fundamental electron charges, charge-$4e$ superconductivity Kivelson  _et al._
(1990); Wu (2005); Berg  _et al._ (2009a); Jiang  _et al._ (2017) is
intrinsically distinct from the conventional charge-$2e$ superconductivity
discovered more than a century ago. One hallmark of charge-4e superconductors
is the magnetic flux quantization with period $hc/4e$, which is half that of
the usual superconducting flux quantum. Unlike the conventional
superconductivity that is well described by the seminal
Bardeen–Cooper–Schrieffer (BCS) theory, many properties of charge-$4e$
superconductivity remain less well understood. The mystery of this phase is
not only because it defies any BCS-type analysis, but also due to the lack of
experimental realization so far. Previous studies in Ref. Berg  _et al._
(2009a) suggested that the charge-$4e$ superconductivity can arise as a
vestigial long-range order above the transition temperature of a pair-density-
wave (PDW) superconductor Fulde and Ferrell (1964); Larkin and Ovchinnikov
(1965); Radzihovsky and Vishwanath (2009); Agterberg and Tsunetsugu (2008);
Wang  _et al._ (2010); Cho  _et al._ (2012); Lee (2014); Maciejko and
Nandkishore (2014); Jian  _et al._ (2015, 2017, 2020); Han  _et al._ (2020);
Agterberg  _et al._ (2020), whose order parameter varies periodically in the
real space. The additional breaking of translational symmetry in the PDW state
is essential in realizing the charge-$4e$ superconductivity as it provides
nontrivial topological defects—dislocations. As a consequence, if the
dislocation proliferates as the temperature is raised, the equilibrium state
will restore the translational symmetry by doubling the original charge-$2e$
condensate, and therefore lead to the charge-$4e$ phase.

The crucial ingredient of the underlying condensate having more than one
quantum number Chung  _et al._ (2012); Moon (2012); Xu and Balents (2018)
provides a general guide to search for charge-$4e$ superconductors. It was
proposed that in certain two-dimensional high-temperature superconductors
(e.g. La2-xBaxCuO4 and La1.6-xNd0.4SrxCuO4 Fujita  _et al._ (2004); Tranquada
_et al._ (2008); Hücker  _et al._ (2011)) exhibiting stripe superconducting
orders Berg  _et al._ (2007, 2009b), the charge-$4e$ superconductivity may
occur above the transition temperature Berg  _et al._ (2009a). Nevertheless,
the experimental signature of the charge-$4e$ superconductivity in these high-
temperature superconductors has not been observed so far. Apart from PDW
superconductors, the nematic superconductor that breaks both charge
conservation and lattice rotational symmetry also hosts an order parameter
that carries multiple quantum numbers. The additional topological defect in
nematic superconductors is created by the intersection of domain walls
separating different ground states. More importantly, many experimental
progresses have been made in achieving the nematic superconductivity in
various systems. For instance, recent experiments found that the spin
susceptibility below the superconducting temperature breaks the three-fold
lattice rotational symmetry in doped topological insulator CuxBi2Se3 Matano
_et al._ (2016); Yonezawa  _et al._ (2017); Tao  _et al._ (2018); Yonezawa
(2019), where the breaking of rotational symmetry suggests that the order
parameter in the superconducting phase is a two-dimensional $E_{u}$
representation of point group $D_{3d}$ Fu (2014); Venderbos  _et al._ (2016);
Hecker and Schmalian (2018). Providing the growing experimental evidence of
the nematic superconductivity in various systems, including doped topological
insulators CuxBi2Se3 Matano  _et al._ (2016); Yonezawa  _et al._ (2017); Tao
_et al._ (2018); Yonezawa (2019), SrxBi2Se3 Wang  _et al._ (2019); Kostylev
_et al._ (2020) and NbxBi2Se3 Shen  _et al._ (2017); Asaba  _et al._ (2017),
superconducting topological semimetals PbTaSe2 Le  _et al._ (2020), and more
recently twisted bilayer graphene Kerelsky  _et al._ (2019); Cao  _et al._
(2020), as well as increasing interest in novel charge-$4e$ ordering, it is
thus interesting to ask whether the multiple-component superconducting order
parameters in nematic superconductors is able to realize charge-$4e$ phases
above their charge-$2e$ transition temperature.

In this paper, we answer this question in the affirmative; namely, we provide
two analysis, suitable for 2D and 3D, respectively, to support the possibility
of a vestigial charge-$4e$ superconducting phase from nematic superconductors.
In 2D, because a vortex is a point-like object, we analyze the fate of various
topological defects to determine the phase diagram. After identifying three
distinct topological defects that are responsible for three different orders
out of the nematic superconducting phase, we find that the competition between
the superfluid stiffness and the elastic constant leads to a rich phase
diagram, as shown in Fig. 1. In particular, if the elastic constant of the
material is less than one third of the superfluid stiffness, the energetically
favored nematic vortices are proliferated, resulting in the novel charge-$4e$
phase, when the temperature gets raised above the transition temperature. In
3D, we employ a Ginzburg-Landau theory near the transition point where the
nematic superconducting order parameter is well described by a simple field
theory up to quartic terms. We are able to show analytically that the
effective mass of the charge-$4e$ order is less than that of the nematic order
by treating the three-fold anisotropy perturbatively, which indicates the
charge-$4e$ order is favored. We further confirm our results by numerically
solving the saddle-point equation. Our result in 3D suggest that materials
whose dispersion along the third direction is weaker than the in-plane
dispersion such as in SrxBi2Se3 are promising in realizing the charge-$4e$
phase.

Charge-4$e$ phase from proliferating topological defect.—We start by analyzing
the symmetries in the nematic superconducting phase providing an appropriate
language for clarifying the topological defects, and then turn to
renormalization group analysis of the defect theory. As mentioned above, we
consider the order parameter of the nematic superconducting phase being a two-
component complex boson $\Delta=(\Delta_{x},\Delta_{y})^{T}$ carrying $E_{u}$
representation in $D_{3d}$ group, and each complex component additionally
hosts the $U(1)$ quantum number corresponding to the charge conservation,
i.e., it is the conventional charge-$2e$ condensate. In the phase respecting
time-reversal symmetry, the relative phase between two components
$\Delta_{x},\Delta_{y}$ is pined at $0$ or $\pi$. As we focus on the time-
reversal-invariant nematic superconductor, we hereafter assume the relative
phase between $\Delta_{x}$ and $\Delta_{y}$ is pined at $0$ or $\pi$ in the
analysis of low-energy physics.

In terms of the basis $\Delta_{\pm}=\Delta_{x}\pm i\Delta_{y}$, one can bring
the the four-dimensional field configuration into three phase modes and one
amplitude mode. Moreover, the time reversal invariance makes one of the phase
modes heavier, thus we have two phase modes and one amplitude mode in the low-
energy sector,

$\displaystyle\Delta_{+}=|\Delta|e^{i(\theta+\phi)},\quad\Delta_{-}=|\Delta|e^{i(\theta-\phi)}.$
(1)

where $\theta$ and $\phi$ denote the two phase modes. The $\phi$ field
describes the $U(1)$ rotation between amplitudes of two components (i.e. the
spatial rotation), whereas the $\theta$ field is the $U(1)$ phase conjugated
to the global charge. The two $U(1)$ phases can also be understood by
secondary orders, i.e., the nematic order
$Q\sim\Delta_{-}^{\dagger}\Delta_{+}$ and the charge-$4e$ order
$\Delta_{4e}\sim\Delta_{+}\Delta_{-}$, that transform as

$\displaystyle Q\rightarrow
Qe^{2i\phi},\quad\Delta_{4e}\rightarrow\Delta_{4e}e^{2i\theta}.$ (2)

The single-valueness of nematic superconducting order parameter $\Delta_{\pm}$
uniquely determines that the topological defects are given by
$(\delta\theta,\delta\phi)=(2\pi,0),(0,2\pi),(\pi,\pi)$, where $\delta\theta$,
$\delta\phi$ denote the winding of the phase around a defect. Physically, they
correspond to the superconducting vortex, the nematic double vortex and the
superconducting half-vortex binding with a single nematic vortex,
respectively. In the nematic superconducting phase, the effective Hamiltonian
characterizing the phase modes is

$\displaystyle
H=\frac{\rho}{2}(\partial\theta)^{2}+\frac{\kappa}{2}(\partial\phi)^{2}-g_{6}\cos
6\pi\phi,$ (3)

where $\rho$ is the superfluid stiffness (superfluid density) and $\kappa$ is
the nematic stiffness (elastic constant). Notice that in writing this
Hamiltonian, we have changed the compactification of the phase mode from
$2\pi$ to $1$ for convenience. Note that $\cos 6\pi\phi$ is allowed owing to
the three-fold anisotropy, i.e., the action is invariant under
$\phi\\!\rightarrow\\!\phi\\!+\\!1/3$. Thermal proliferation of topological
defects can be then described by the dual bosons $\tilde{\phi}$ and
$\tilde{\theta}$ Berg  _et al._ (2009a),

$\displaystyle H$ $\displaystyle=$
$\displaystyle\frac{T}{2\rho}(\partial\tilde{\theta})^{2}+\frac{T}{2\kappa}(\partial\tilde{\phi})^{2}-g_{6}\cos
6\pi\phi-g_{2,0}\cos 2\pi\tilde{\theta}$ (4) $\displaystyle-g_{0,2}\cos
2\pi\tilde{\phi}-g_{1,1}\cos\pi\tilde{\theta}\cos\pi\tilde{\phi},$

where $g_{2,0}$, $g_{0,2}$ and $g_{1,1}$ are couplings characterizing the
strength of creating (annihilating) each kind of topological defects,
respectively, and $T$ is the temperature.

Figure 1: Schematic phase diagram of the defect theory. $\kappa$, $\rho$ and
$T$ denote the nematic stiffness, the superfluid stiffness and temperature,
respectively. The solid lines refer to phase boundaries. Blow the dashed line,
the three-fold anisotropy is relevant from the lowest order calculation.
Nevertheless, the true transition to the charge-$4e$ phase belongs to the
three-state Potts model universality class Wu (1982) where the three-fold
anisotropy is relevant at the transition point.

The standard renormalization group flow for Eq. (4) to the lowest-order reads

$\displaystyle\frac{dg_{2,0}}{dl}$ $\displaystyle=$
$\displaystyle\left(2-\frac{\pi\rho}{T}\right)g_{2,0},$ (5)
$\displaystyle\frac{dg_{0,2}}{dl}$ $\displaystyle=$
$\displaystyle\left(2-\frac{\pi\kappa}{T}\right)g_{0,2},$ (6)
$\displaystyle\frac{dg_{1,1}}{dl}$ $\displaystyle=$
$\displaystyle\left[2-\frac{\pi}{4T}(\rho+\kappa)\right]g_{0,2},$ (7)
$\displaystyle\frac{dg_{6}}{dl}$ $\displaystyle=$
$\displaystyle\left(2-\frac{9\pi T}{\kappa}\right)g_{6}.$ (8)

The coefficient in front of the coupling on the right-hand side of the
renormalization group equation determines whether the corresponding process is
relevant or not. For instance, if $\frac{T}{\rho}>\frac{\pi}{2}$, the creation
and annihilation process of the superconducting vortex is relevant, leading to
the proliferation of superconducting phase defects and, consequently,
destroying the superconducting phase coherence driving the system out of
superconducting orders. In this way, the phase diagram as shown in Fig. 1 is
mapped out by the renormalization group equations. Clearly, when the nematic
stiffness is less than the superfluid stiffness, more specifically
$\kappa<\frac{\rho}{3}$, it is easier to create nematic vortices than
superconducting vortices such that raising temperature can more efficiently
proliferate the nematic vortices, restoring the lattice rotational symmetry
but not the $U(1)$ charge symmetry. As a result, the charge-$4e$ order emerges
since any charge-$2e$ condensate from the nematic superconductivity breaks
lattice rotational symmetry whereas the charge-$4e$ order is blind to the
nematic vortex as shown by the symmetry transformation law in Eq. (2).

Besides the interesting charge-$4e$ phase, the competition between the nematic
stiffness and the superfluid stiffness results in a rich phase diagram in Fig.
1. Namely, when $\frac{1}{3}<\frac{\kappa}{\rho}<3$, raising the temperature
causes a direct transition from the nematic superconductivity to the normal
phase since proliferating a superconducting half-vortex bounded with a nematic
vortex is favored, whereas, when $\frac{\kappa}{\rho}>3$, a vestigial of
nematic phase emerges since proliferating the normal superconducting phase
vortex is favored Cho  _et al._ (2020).

Charge-$4e$ phase from superconducting fluctuations.— In 3D, the Ginzburg-
Landau theory works better since the quantum fluctuation is generally
suppressed as dimension increases, so we expect a mean-field analysis to 3D
nematic superconductivity near the phase boundary can describe the essential
physics. Such a theory was also used in Ref. Hecker and Schmalian (2018) to
investigate the vestigial nematic order, but the authors did not analyze the
possibility of charge-$4e$ orders. The Ginzburg-Landau theory of the nematic
superconducting order parameter $\Delta=(\Delta_{x},\Delta_{y})^{T}$ near the
phase boundary reads Hecker and Schmalian (2018)

$\displaystyle
S=\int_{p}\Delta^{\dagger}G_{p}^{-1}\Delta+\int_{x}(u(\Delta^{\dagger}\Delta)^{2}+v(\Delta^{\dagger}\tau^{y}\Delta)^{2}),$
(9)

where $\int_{p}=\int\frac{dp^{3}}{(2\pi)^{3}}$, $\int_{x}=\int d^{3}x$, and

$\displaystyle
G_{p}^{-1}=m_{0}(\bm{p})\tau^{0}+m_{1}(\bm{p})\tau^{z}+m_{2}(\bm{p})\tau^{x},$
(10)

is the inverse propagator of the nematic superconducting order parameter with
$m_{0}(\bm{p})=d_{\parallel}(p_{x}^{2}+p_{y}^{2})+d_{z}p_{z}^{2}+r_{0}$,
$m_{1}(\bm{p})=d^{\prime}(p_{x}^{2}-p_{y}^{2})+\bar{d}p_{y}p_{z}$, and
$m_{2}(\bm{p})=d^{\prime}2p_{x}p_{y}+\bar{d}p_{x}p_{z}$. Notice that while
$m_{0}(\bm{p})$ enjoys a continuous rotational symmetry, $m_{1}(\bm{p})$ and
$m_{2}(\bm{p})$ lower it down to three-fold rotation. $r_{0}$, $u$, $v$ are
real parameters allowed by the symmetry and $d_{\parallel}$, $d_{z}$,
$d^{\prime}$, $\bar{d}$ characterize the kinetic energy. In the following, we
consider $v>0$, $u>0$, as the case for $v<0$ favors the chiral
superconductivity that breaks time reversal symmetry. It is crucial to realize
that there are two ways given by the following two different Fierz identities
to decouple the last quartic term in Eq. (9)

$\displaystyle\tau^{y}_{\alpha\beta}\tau^{y}_{\gamma\delta}$ $\displaystyle=$
$\displaystyle
2\delta_{\alpha\delta}\delta_{\beta\gamma}-\delta_{\alpha\beta}\delta_{\gamma\delta}-\tau^{x}_{\alpha\beta}\tau^{x}_{\gamma\delta}-\tau^{z}_{\alpha\beta}\tau^{z}_{\gamma\delta},$
(11) $\displaystyle\tau^{y}_{\alpha\beta}\tau^{y}_{\gamma\delta}$
$\displaystyle=$
$\displaystyle\delta_{\alpha\beta}\delta_{\gamma\delta}-\delta_{\alpha\gamma}\delta_{\beta\delta}.$
(12)

This leads to either nematic channel for the first one Eq. (11) or charge-$4e$
channel for the second one Eq. (12). After the decoupling, one is able to
integrate out the quadratic nematic superconducting order parameter. We leave
the details to the Supplemental Materials SM , and present the main results
here. The Ginzburg-Landau theory for the nematic order and the charge-$4e$
order (assuming to be homogeneous in real space) are given by

$\displaystyle S_{nem}$ $\displaystyle=$
$\displaystyle\text{Tr}\log\chi_{p}^{-1}+\int_{x}\left(\frac{1}{4v}(Q_{x}^{2}+Q_{y}^{2})-\frac{1}{4u^{\prime}}R^{2}\right),$
(13) $\displaystyle S_{4e}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\text{Tr}\log
D_{p}^{-1}+\int_{x}\left(\frac{1}{4v}|\Delta_{4e}|^{2}-\frac{1}{4u^{\prime}}R^{2}\right),$
(14)

where $u^{\prime}=u+v$, $Q_{x},Q_{y}$ are the nematic orders, $\Delta_{4e}$ is
the charge-$4e$ order, and $R\sim\Delta^{\dagger}\Delta$ is a scalar that
decouples the first quartic interaction in Eq. (9). The $\text{Tr}\log$ term
comes from integrating out the nematic superconducting order parameter, and
$\chi_{p}^{-1}=G_{p}^{-1}+R+Q_{x}\tau^{z}+Q_{y}\tau^{x}$, and
$D_{p}^{-1}=G^{-1}_{p}+R+\Delta_{4e}\rho^{+}+\Delta_{4e}^{\ast}\rho^{-}$ are
the modified inverse propagators in the presence of the nematic order and the
charge-$4e$ order, respectively. $\rho^{\pm}=\frac{1}{2}(\rho^{x}\pm
i\rho^{y})$ where $\rho^{x}$ and $\rho^{y}$ are Pauli matrices acting on the
“Nambu” space, i.e., $\Delta_{4e}$ and $\Delta_{4e}^{\ast}$, and the extra
factor $\frac{1}{2}$ in front of the $\text{Tr}\log$ term in Eq. (14) is due
to the redundancy of the enlarged Nambu space.

The Fierz identity and the bosonic nature of the nematic superconducting order
parameter $\Delta$ cause the similar Ginzburg-Landau theory for the nematic
order Eq. (13) and the charge-$4e$ orders Eq. (14). The only difference comes
from the $\text{Tr}\log$ term as the nematic order carries a real space
(nematic) quantum number while the charge-$4e$ order carries a Nambu space
quantum number. Unlike the fermionic field, the bosonic nematic
superconducting field is blind to the Nambu space as one can easily observe
that the inverse propagator $G_{p}^{-1}$ (10) is an identity in the Nambu
space. Thus if one artificially turns off $m_{1}$ and $m_{2}$ such that the
nematic superconducting field has an accidental continuous real space
rotational symmetry at the leading order, the nematic order and the
charge-$4e$ order will be degenerate since now the nematic superconducting
field is blind to both nematic space and Nambu space (the factor of
$\frac{1}{2}$ in Eq. (14) is perfectly compensated by the enlarged Nambu
space). So even without any calculation, one can identify a specific albeit
artificial point where the charge-$4e$ order is degenerate with the nematic
order. Now the question is how the nonvanishing $m_{1}$ and $m_{2}$,
equivelantly $d^{\prime}$ and $\bar{d}$, affect the two instabilities. In the
following, we provide an analytical result showing the charge-$4e$ order is
favored over the nematic order when $d^{\prime}\gg\bar{d}$, which makes
quasi-2D materials whose out-of-plane dispersion is much weaker than the in-
plane dispersion promising candidates for realizing charge-$4e$ phases. In
addition, we also numerically solve the saddle-point equation for both nematic
orders and charge-$4e$ orders, which confirms the analytical result.

Figure 2: The Feynman diagrams that contribute to the effective mass of the
charge-$4e$ order and the nematic order, respectively. The difference of (a)
and (b) arises from the different vertices denoted by $\bm{\rho}$ and
$\bm{\tau}$.

Near the transition point, the effective Ginzburg-Landau theory can be
obtained by assuming order parameters small and expanding the action order by
order. The effective mass that is crucial in determining the phase boundary is
given by evaluating the Feynman diagram in Fig. 2, and as we discussed above
the difference lies in the vertices of the two orders as shown in Fig. 2 and
Fig. 2. Since nonvanishing $m_{1}$ and $m_{2}$ are diagonal in the Nambu space
the contributions from the $\tau^{z}$ and $\tau^{x}$ kinetic energy are
additive to the charge-$4e$ order’s mass; however, since they are not diagonal
in the nematic space the contributions tend to cancel each other. We calculate
the effective mass perturbatively in $\bar{d}$, and the lowest-order result is

$\displaystyle m_{4e}$ $\displaystyle=$ $\displaystyle m_{nem}-\delta m_{2},$
(15) $\displaystyle\delta m_{2}$ $\displaystyle=$
$\displaystyle\frac{1}{d_{\parallel}\sqrt{d_{z}(R+r_{0})}}\frac{2\gamma-(1-\gamma^{2})\log\left(\frac{1+\gamma}{1-\gamma}\right)}{32\pi\gamma(1-\gamma^{2})},$
(16)

where $\gamma=\frac{d^{\prime}}{d_{\parallel}}$. It is not hard to see that
$\delta m_{2}$ is positive and monotonic in $0<\gamma<1$, indicating the
anisotropy induced by nonvanishing $d^{\prime}$ favors the charge-$4e$ order.

However, this enhancement of the charge-$4e$ instability needs to compete with
the cubic potential of nematic orders arised from the three-fold anisotropy.
The cubic term appears at least at the quadratic order in $\bar{d}$ and the
linear order in $d^{\prime}$ (actually, the dependence of $d^{\prime}$ at
quadratic order of $\bar{d}$ can be obtained explicitly), and through explicit
calculation, the cubic term is

$\displaystyle S_{nem}$ $\displaystyle\ni$ $\displaystyle
w\int(Q_{x}^{3}-3Q_{x}Q_{y}^{2}),$ (17) $\displaystyle w$
$\displaystyle\approx$
$\displaystyle\frac{1}{960\pi}\frac{1}{[d_{z}(R+r_{0})]^{3/2}}\frac{\bar{d}^{2}d^{\prime}}{d_{\parallel}^{3}}.$
(18)

It is a simple matter of fact that the cubic term enhances the instability at
least in the quadratic orders, i.e., $\delta m_{3}\propto
w^{2}\propto\bar{d}^{4}$, so the lowest order is the quartic order in
$\bar{d}$. On the other hand, the enhancement of the charge-$4e$
superconductivity $\delta m_{2}$ appears already in the zeroth order of
$\bar{d}$. For the systems described by Eq. (9), there exists a parameter
regime, where $\bar{d}\ll d^{\prime}$, such that the charge-$4e$ order is
favored over the nematic order. Notice that we have neglected the sixth order
terms like $\left(\Delta^{\dagger}(\tau^{z}+i\tau^{x})\Delta\right)^{3}+h.c.$
in Eq. (9) that are allowed by the three-fold anisotropy, and these terms can
also render cubic contributions to the nematic order effective action Eq.
(17). These higher order terms are in general irrelevant for the long
wavelength physics near and slightly above the nematic superconducting
transition temperature at which the Ginzburg-Landau theory Eq. (9) can apply.
(Whereas if the temperature is lower than the transition temperature, this
becomes dangerously irrelevant Jian  _et al._ (2020) and selects one of the
degenerate ground states). Unlike the universal physics that can be obtained
with the help of the Ginzburg-Landau action in Eq. (9), it is not clear away
from the applicability of such an action what the fate of the competition
between the charge-$4e$ order and the nematic order is. Nevertheless, the 2D
calculation shown in Fig. 1 provides a complementary qualitative
understanding: the large anisotropy from higher-order potentials can increase
the nematic stiffness and prohibit the charge-$4e$ phase Cho  _et al._ (2020).

Figure 3: The mean field phase diagram of the charge-$4e$ order (a) and the
nematic order (b). $|\Delta_{4e}|$, and $Q=\sqrt{Q_{x}^{2}+Q_{y}^{2}}$ are the
amplitude of the charge-$4e$ and the nematic order respectively. $\delta
r=r_{0}-r_{0}^{\ast}\propto T-T^{\ast}$, where $r^{\ast}$ ($T^{\ast}$) denotes
the transition point (transition temperature) without vestigial orders. The
parameters are $d_{\parallel}=1,d^{\prime}=0.5,d_{z}=0.1,\bar{d}=0.1,v=1,u=4$.

To further support our result, we carry out a numerical calculation of the
saddle-point equations in both the charge-$4e$ channel and the nematic
channel. We leave the saddle-point equation in the Supplemental Materials SM .
The phase diagram is shown in Fig. 3, where $|\Delta_{4e}|$, and
$Q=\sqrt{Q_{x}^{2}+Q_{y}^{2}}$ are the amplitude of the charge-$4e$ and the
nematic order respectively. $\delta r=r_{0}-r_{0}^{\ast}\propto T-T^{\ast}$,
where $r^{\ast}$ ($T^{\ast}$) denotes the transition point (transition
temperature) without vestigial orders. That the system is in favor of the
charge-$4e$ phases when $\bar{d}\ll d^{\prime}$ is demonstrated explicitly by
$T_{nem}<T_{4e}$ in Fig. 3. One should also notice the first-order jump of the
nematic transition due to the cubic anisotropy, on the contrary, the phase
transition to the nematic phase at two dimensions is continuous due to the
strong quantum fluctuations Wu (1982).

Conclusions.—We have revealed the promosing possibility of realizing a novel
charge-$4e$ superconducting phase above the transition temperature of time-
reversal-invariant nematic superconductors in 2D and quasi-2D. Since
experimental evidences of rotational symmetry breaking superconducting phases
in various systems, including doped topological insulators CuxBi2Se3 Matano
_et al._ (2016); Yonezawa  _et al._ (2017); Tao  _et al._ (2018); Yonezawa
(2019), SrxBi2Se3 Wang  _et al._ (2019); Kostylev  _et al._ (2020), and
NbxBi2Se3 Shen  _et al._ (2017); Asaba  _et al._ (2017), superconducting
topological semimetal PbTaSe2 Le  _et al._ (2020), and more recently twisted
bilayer graphene Kerelsky  _et al._ (2019); Cao  _et al._ (2020), have
accumulated, we believe that our results pave an important step toward
experimentally realizing charge-$4e$ phases in quantum materials.

Acknowledgement.— We thank Wen Huang for helpful discussions. This work is
supported in part by the NSFC under Grant No. 11825404 (HY), the MOSTC under
Grant Nos. 2016YFA0301001 and 2018YFA0305604 (HY), and the Strategic Priority
Research Program of Chinese Academy of Sciences under Grant No. XDB28000000
(HY). SKJ is supported by the Simons Foundation via the It From Qubit
Collaboration.

Note added: While the present paper was close to be completed, we notice an
interesting work on a similar topic in Ref. Fernandes and Fu (2021). Our
results of Ginzburg-Landau theory analysis in three dimensions are
qualitatively consistent with ones in Ref. Fernandes and Fu (2021). We
additionally analyzed the topological defect theory in two dimensions which
supports the possibility of realizing charge-4e phase.

## References

  * Kivelson  _et al._ (1990) S. Kivelson, V. Emery, and H. Lin, Phys. Rev. B 42, 6523 (1990).
  * Wu (2005) C. Wu, Phys. Rev. Lett. 95, 266404 (2005).
  * Berg  _et al._ (2009a) E. Berg, E. Fradkin, and S. A. Kivelson, Nature Physics 5, 830 (2009a).
  * Jiang  _et al._ (2017) Y.-F. Jiang, Z.-X. Li, S. A. Kivelson, and H. Yao, Phys. Rev. B 95, 241103 (2017).
  * Fulde and Ferrell (1964) P. Fulde and R. A. Ferrell, Phys. Rev. 135, A550 (1964).
  * Larkin and Ovchinnikov (1965) A. Larkin and Y. N. Ovchinnikov, JETP 20, 762 (1965).
  * Radzihovsky and Vishwanath (2009) L. Radzihovsky and A. Vishwanath, Phys. Rev. Lett. 103, 010404 (2009).
  * Agterberg and Tsunetsugu (2008) D. Agterberg and H. Tsunetsugu, Nature Physics 4, 639 (2008).
  * Wang  _et al._ (2010) C. Wang, C. Gao, C.-M. Jian, and H. Zhai, Phys. Rev. Lett. 105, 160403 (2010).
  * Cho  _et al._ (2012) G. Y. Cho, J. H. Bardarson, Y.-M. Lu, and J. E. Moore, Phys. Rev. B 86, 214514 (2012).
  * Lee (2014) P. A. Lee, Phys. Rev. X 4, 031017 (2014).
  * Maciejko and Nandkishore (2014) J. Maciejko and R. Nandkishore, Phys. Rev. B 90, 035126 (2014).
  * Jian  _et al._ (2015) S.-K. Jian, Y.-F. Jiang, and H. Yao, Phys. Rev. Lett. 114, 237001 (2015).
  * Jian  _et al._ (2017) S.-K. Jian, C.-H. Lin, J. Maciejko, and H. Yao, Phys. Rev. Lett. 118, 166802 (2017).
  * Jian  _et al._ (2020) S.-K. Jian, M. M. Scherer, and H. Yao, Phys. Rev. Research 2, 013034 (2020).
  * Han  _et al._ (2020) Z. Han, S. A. Kivelson, and H. Yao, Phys. Rev. Lett. 125, 167001 (2020).
  * Agterberg  _et al._ (2020) D. F. Agterberg, J. S. Davis, S. D. Edkins, E. Fradkin, D. J. Van Harlingen, S. A. Kivelson, P. A. Lee, L. Radzihovsky, J. M. Tranquada, and Y. Wang, Annual Review of Condensed Matter Physics 11, 231 (2020).
  * Chung  _et al._ (2012) S. B. Chung, S. Raghu, A. Kapitulnik, and S. A. Kivelson, Physical Review B 86, 064525 (2012).
  * Moon (2012) E.-G. Moon, Phys. Rev. B 85, 245123 (2012).
  * Xu and Balents (2018) C. Xu and L. Balents, Physical review letters 121, 087001 (2018).
  * Fujita  _et al._ (2004) M. Fujita, H. Goka, K. Yamada, J. Tranquada, and L. Regnault, Phys. Rev. B 70, 104517 (2004).
  * Tranquada  _et al._ (2008) J. M. Tranquada, G. D. Gu, M. Hücker, Q. Jie, H.-J. Kang, R. Klingeler, Q. Li, N. Tristan, J. S. Wen, G. Y. Xu, Z. J. Xu, J. Zhou, and M. v. Zimmermann, Phys. Rev. B 78, 174529 (2008).
  * Hücker  _et al._ (2011) M. Hücker, M. v. Zimmermann, G. Gu, Z. Xu, J. Wen, G. Xu, H. Kang, A. Zheludev, and J. M. Tranquada, Phys. Rev. B 83, 104506 (2011).
  * Berg  _et al._ (2007) E. Berg, E. Fradkin, E.-A. Kim, S. A. Kivelson, V. Oganesyan, J. M. Tranquada, and S.-C. Zhang, Phys. Rev. Lett. 99, 127003 (2007).
  * Berg  _et al._ (2009b) E. Berg, E. Fradkin, and S. A. Kivelson, Phys. Rev. B 79, 064515 (2009b).
  * Matano  _et al._ (2016) K. Matano, M. Kriener, K. Segawa, Y. Ando, and G.-q. Zheng, Nature Physics 12, 852 (2016).
  * Yonezawa  _et al._ (2017) S. Yonezawa, K. Tajiri, S. Nakata, Y. Nagai, Z. Wang, K. Segawa, Y. Ando, and Y. Maeno, Nature Physics 13, 123 (2017).
  * Tao  _et al._ (2018) R. Tao, Y.-J. Yan, X. Liu, Z.-W. Wang, Y. Ando, Q.-H. Wang, T. Zhang, and D.-L. Feng, Phys. Rev. X 8, 041024 (2018).
  * Yonezawa (2019) S. Yonezawa, Condensed Matter 4, 2 (2019).
  * Fu (2014) L. Fu, Phys. Rev. B 90, 100509 (2014).
  * Venderbos  _et al._ (2016) J. W. Venderbos, V. Kozii, and L. Fu, Phys. Rev. B 94, 180504 (2016).
  * Hecker and Schmalian (2018) M. Hecker and J. Schmalian, npj Quantum Materials 3, 1 (2018).
  * Wang  _et al._ (2019) J. Wang, K. Ran, S. Li, Z. Ma, S. Bao, Z. Cai, Y. Zhang, K. Nakajima, S. Ohira-Kawamura, P. Čermák, A. Schneidewind, S. Y. Savrasov, X. Wan, and J. Wen, Nature communications 10, 2802 (2019).
  * Kostylev  _et al._ (2020) I. Kostylev, S. Yonezawa, Z. Wang, Y. Ando, and Y. Maeno, Nature communications 11, 1 (2020).
  * Shen  _et al._ (2017) J. Shen, W.-Y. He, N. F. Q. Yuan, Z. Huang, C.-w. Cho, S. H. Lee, Y. San Hor, K. T. Law, and R. Lortz, npj Quantum Materials 2, 1 (2017).
  * Asaba  _et al._ (2017) T. Asaba, B. Lawson, C. Tinsman, L. Chen, P. Corbae, G. Li, Y. Qiu, Y. S. Hor, L. Fu, and L. Li, Physical Review X 7, 011009 (2017).
  * Le  _et al._ (2020) T. Le, Y. Sun, H.-K. Jin, L. Che, L. Yin, J. Li, G. Pang, C. Xu, L. Zhao, S. Kittaka, T. Sakakibara, K. Machida, R. Sankar, H. Yuan, G. Chen, X. Xu, S. Li, Y. Zhou, and X. Lu, Science Bulletin 65, 1349 (2020).
  * Kerelsky  _et al._ (2019) A. Kerelsky, L. J. McGilly, D. M. Kennes, L. Xian, M. Yankowitz, S. Chen, K. Watanabe, T. Taniguchi, J. Hone, C. Dean, A. Rubio, and A. N. Pasupathy, Nature 572, 95 (2019).
  * Cao  _et al._ (2020) Y. Cao, D. Rodan-Legrain, J. M. Park, F. N. Yuan, K. Watanabe, T. Taniguchi, R. M. Fernandes, L. Fu, and P. Jarillo-Herrero, arXiv:2004.04148 (2020).
  * Wu (1982) F.-Y. Wu, Rev. Mod. Phys. 54, 235 (1982).
  * Cho  _et al._ (2020) C.-w. Cho, J. Shen, J. Lyu, O. Atanov, Q. Chen, S. H. Lee, Y. San Hor, D. J. Gawryluk, E. Pomjakushina, M. Bartkowiak, _et al._ , Nature communications 11, 1 (2020).
  * (42) See Supplemental Materials for the technical details of derivation the saddle-point equation.
  * Fernandes and Fu (2021) R. M. Fernandes and L. Fu, arXiv:2101.07943 (2021).

Supplemental Materials

There are two ways to decouple the quartic term in Eq. (9), given by the
following two different Fierz identities

$\displaystyle\tau^{y}_{\alpha\beta}\tau^{y}_{\gamma\delta}$ $\displaystyle=$
$\displaystyle
2\delta_{\alpha\delta}\delta_{\beta\gamma}-\delta_{\alpha\beta}\delta_{\gamma\delta}-\tau^{x}_{\alpha\beta}\tau^{x}_{\gamma\delta}-\tau^{z}_{\alpha\beta}\tau^{z}_{\gamma\delta},$
(S1) $\displaystyle\tau^{y}_{\alpha\beta}\tau^{y}_{\gamma\delta}$
$\displaystyle=$
$\displaystyle\delta_{\alpha\beta}\delta_{\gamma\delta}-\delta_{\alpha\gamma}\delta_{\beta\delta},$
(S2)

which lead to two exact rewritings of Eq. (9)

$\displaystyle S$ $\displaystyle=$
$\displaystyle\int_{p}\Delta^{\dagger}G^{-1}_{p}\Delta+\int_{x}(u^{\prime}(\Delta^{\dagger}\Delta)^{2}-v[(\Delta^{\dagger}\tau^{x}\Delta)^{2}+(\Delta^{\dagger}\tau^{z}\Delta)^{2}]),$
(S3) $\displaystyle S$ $\displaystyle=$
$\displaystyle\int_{p}\Delta^{\dagger}G^{-1}_{p}\Delta+\int_{x}(u^{\prime}(\Delta^{\dagger}\Delta)^{2}-v|\Delta^{T}\Delta|^{2}),$
(S4)

where $u^{\prime}=u+v$.

The first equation (S3) is readily decoupled by nematic orders,
$(Q_{x},Q_{y})$,

$\displaystyle S$ $\displaystyle=$
$\displaystyle\int_{x}\left(\Delta^{\dagger}(G^{-1}+R+Q_{x}\tau^{z}+Q_{y}\tau^{x})\Delta+\frac{1}{4v}(Q_{x}^{2}+Q_{y}^{2})-\frac{1}{4u^{\prime}}R^{2})\right),$
(S5)

where $R$ is another boson field that decouple $u^{\prime}$ term, and it does
not break any symmetry. On the other hand, the second equation (S4) is readily
decoupled by charge-$4e$ SC orders, $(\Delta_{4e},\Delta_{4e}^{\ast})$,

$\displaystyle S$ $\displaystyle=$
$\displaystyle\int_{x}\left(\frac{1}{2}\hat{\Delta}^{\dagger}\left(\begin{array}[]{cccc}G^{-1}+R&\Delta_{4e}\\\
\Delta_{4e}^{\ast}&G^{-1}+R\end{array}\right)\hat{\Delta}+\frac{1}{4v}|\Delta_{4e}|^{2}-\frac{1}{4u^{\prime}}R^{2}\right)$
(S8)

where $\hat{\Delta}=(\Delta^{\dagger},\Delta^{T})$.

Since now the nematic SC order is quadratic in the action, it is
straightforward to integrate it over.

$\displaystyle S$ $\displaystyle=$
$\displaystyle\text{Tr}\log(G^{-1}+R+Q_{x}\tau^{z}+Q_{y}\tau^{x})+\int_{x}\left(\frac{1}{4v}(Q_{x}^{2}+Q_{y}^{2})-\frac{1}{4u^{\prime}}R^{2}\right),$
(S9) $\displaystyle S$ $\displaystyle=$
$\displaystyle\frac{1}{2}\text{Tr}\log\left(\begin{array}[]{cccc}G^{-1}+R&\Delta_{4e}\\\
\Delta_{4e}^{\ast}&G^{-1}+R\end{array}\right)+\int_{x}\left(\frac{1}{4v}|\Delta_{4e}|^{2}-\frac{1}{4u^{\prime}}R^{2}\right).$
(S12)

From the effective action (13), the saddle-point equation for homogeneous
nematic orders, $Q_{i}(\bm{x})=Q_{i}$, $i=x,y$ reads

$\displaystyle R$ $\displaystyle=$ $\displaystyle
2u^{\prime}\Delta_{0}^{\dagger}\Delta_{0}+2u^{\prime}\int_{p}\text{Tr}(G_{p}^{-1}+R+Q_{x}\tau^{z}+Q_{y}\tau^{x})^{-1},$
(S13) $\displaystyle Q_{x}$ $\displaystyle=$
$\displaystyle-2v\Delta_{0}^{\dagger}\tau^{z}\Delta_{0}-2v\int_{p}\text{Tr}(G_{p}^{-1}+R+Q_{x}\tau^{z}+Q_{y}\tau^{x})^{-1}\tau^{z},$
(S14) $\displaystyle Q_{y}$ $\displaystyle=$
$\displaystyle-2v\Delta_{0}^{\dagger}\tau^{x}\Delta_{0}-2v\int_{p}\text{Tr}(G_{p}^{-1}+R+Q_{x}\tau^{z}+Q_{y}\tau^{x})^{-1}\tau^{x},$
(S15) $\displaystyle 0$ $\displaystyle=$
$\displaystyle\chi_{0}^{-1}\Delta_{0}.$ (S16)

On the other hand, the saddle-point equation for homogeneous charge-$4e$ SC
orders $\Delta_{4e}(\bm{x})=\Delta_{4e}^{\prime}+i\Delta_{4e}^{\prime\prime}$
where $\Delta_{4e}^{\prime},\Delta_{4e}^{\prime\prime}\in\mathbb{R}$ from Eq.
(14) reads

$\displaystyle R$ $\displaystyle=$ $\displaystyle
2u^{\prime}\Delta_{0}^{\dagger}\Delta_{0}+u^{\prime}\int_{p}\text{Tr}\left(\begin{array}[]{cccc}G^{-1}_{p}+R&\Delta_{4e}^{\prime}+i\Delta_{4e}^{\prime\prime}\\\
\Delta_{4e}^{\prime}-i\Delta_{4e}^{\prime\prime}&G^{-1}_{p}+R\end{array}\right)^{-1}$
(S19) $\displaystyle\Delta_{4e}^{\prime}$ $\displaystyle=$
$\displaystyle-2v\text{Re}[\Delta_{0}^{T}\Delta_{0}]-v\int_{p}\text{Tr}\left(\begin{array}[]{cccc}G^{-1}_{p}+R&\Delta_{4e}^{\prime}+i\Delta_{4e}^{\prime\prime}\\\
\Delta_{4e}^{\prime}-i\Delta_{4e}^{\prime\prime}&G^{-1}_{p}+R\end{array}\right)^{-1}\rho^{x},$
(S22) $\displaystyle\Delta_{4e}^{\prime\prime}$ $\displaystyle=$
$\displaystyle-2v\text{Im}[\Delta_{0}^{T}\Delta_{0}]-v\int_{p}\text{Tr}\left(\begin{array}[]{cccc}G^{-1}_{p}+R&\Delta_{4e}^{\prime}+i\Delta_{4e}^{\prime\prime}\\\
\Delta_{4e}^{\prime}-i\Delta_{4e}^{\prime\prime}&G^{-1}_{p}+R\end{array}\right)^{-1}(-\rho^{y}),$
(S25) $\displaystyle 0$ $\displaystyle=$
$\displaystyle\chi_{0}^{-1}\Delta_{0}.$ (S26)

In the saddle-point equations, we also assume $R(\bm{x})=R$ is homogeneous.
The saddle-point equations are solved numerically, and the results are shown
in Fig. 3.