# Exploring the structure-property relations of thin-walled, 2D extruded
lattices using neural networks

Junyan He1 Shashank Kushwaha1 Diab Abueidda2 Iwona Jasiuk1 1 Department of
Mechanical Science and Engineering, University of Illinois at Urbana-
Champaign, Champaign, IL, USA
2 National Center for Supercomputing Applications, University of Illinois at
Urbana-Champaign, Champaign, IL, USA<EMAIL_ADDRESS>
###### Abstract

This paper investigates the structure-property relations of thin-walled
lattices under dynamic longitudinal compression, characterized by their cross-
sections and heights. These relations elucidate the interactions of different
geometric features of a design on mechanical response, including energy
absorption. We proposed a combinatorial, key-based design system to generate
different lattice designs and used the finite element method to simulate their
response with the Johnson-Cook material model. Using an autoencoder, we
encoded the cross-sectional images of the lattices into latent design feature
vectors, which were supplied to the neural network model to generate
predictions. The trained models can accurately predict lattice energy
absorption curves in the key-based design system and can be extended to new
designs outside of the system via transfer learning.

## Graphical abstract

###### keywords:

Thin-walled lattices , Structure-property relations , Johnson-Cook model ,
Neural networks

††journal: Computers & Structures

*[inlinelist,1]label=(),

##  1 Introduction

Lattice-filled sandwich panels see increasingly widespread use as energy-
absorbing structures in different engineering applications such as energy
absorbers [1, 2], sacrificial cladding [3, 4], and armor plates [5, 6].
Previous studies have established that the core design plays a significant
role in the energy absorption capabilities of the sandwich panel [7]. The
honeycomb [8] lattice structure, along with many of its variants, such as the
square honeycomb [9], bio-inspired honeycombs [10], and hierarchical
honeycombs [11], have been studied extensively in the literature as cores of
sandwich panels. They are examples of thin-walled 2D extruded lattices, whose
wall thickness is significantly smaller than their in-plane dimensions, and
the unit cells can be described by their cross-sectional designs and heights.
To generate optimized designs that maximize specific energy absorption,
parameter optimizations have been performed on design variables such as
thickness [12], unit cell size [13], wave periodicity, and amplitude (for some
bio-inspired honeycombs) [10].

The aforementioned sizing optimization efforts and the construction of
response surfaces for the cores [13, 14] provide insights into the structure-
property relations relating the lattice structure to key performance metrics
like stress-strain curve and energy absorption during deformation for a
limited set of honeycomb-like lattice core designs. It is of great research
and industrial interests to further extend such structure-property relations
to a diverse set of lattice core designs with different cross-sectional
geometries, and in particular, to elucidate how different geometric features
such as the addition of wavy unit-cell walls, hierarchical designs, auxetic
designs, and fractal designs combine and interact to affect the specific
energy absorption of the lattice core design. The exploration of structure-
property relations inherently involves surveying many different lattice core
designs. Most current thin-walled lattice core designs in the literature were
generated based on engineering intuition, experiments, and/or bio-inspiration
[15, 16, 17]. Other studies leveraged optimization-based methods like the
hybrid cellular automata [18, 19], ground structure approach [20], and the ant
colony optimization method [21] to systematically generate new, optimized
cross sections for lattices under longitudinal compression. Nonetheless, there
is a lack of a systematic compilation of the mechanical response and energy
absorption characteristics for these new designs and an attempt to reveal
potential structure-property relations for various lattice cross-sectional
designs.

Our current work aims to develop a systematic framework to generate distinct
lattice cross-section designs for longitudinal compression. To explore the
structure-property relations, the considered structures do not need to have
optimized energy absorption characteristics. Thus, we chose to generate new
lattice cross-section designs via a simple combinatorial framework.
Combinatorial methods have been applied to generate trussed-based lattices
[22, 23, 24], triply periodic minimal surface lattices [25, 26], and thin-
walled structures [27], which create new designs by combining different pre-
selected geometric descriptors. It is a novel attempt to leverage a
combinatorial framework to generate 2D extruded lattice core designs for
longitudinal compression systematically. The second objective of this work is
to approximate the underlying structure-property relations via a neural
network (NN) model. NN models have seen increasing use in solid mechanics,
such as to predict stress-strain curve and toughness of composites [28, 29,
30, 31, 32] and to predict properties of lattices and rapidly evaluate
performance [23, 33, 34, 35, 36]. However, NN models have not been widely
applied to approximate the structure-property relations of thin-walled lattice
cores.

This paper is organized as follows: Section 2 presents an overview of the
combinatorial generation framework, numerical simulation, preprocessing of the
design images, and the architecture of the NN model. Section 3 presents the
results, and Section 4 discusses the quality of the image preprocessing and NN
predictions. Section 5 summarizes the outcomes and highlights possible future
works.

##  2 Methods

###  2.1 A framework to create lattices using geometric descriptors

Diverse design space of lattice cross-sections is needed such that the
subsequent NN can extract trends from various lattice architectures. Drawing
inspiration from the truss descriptors used to construct truss-based lattices
in Bastek et al. [23] and Zok et al. [37], we define a set of geometric
descriptors for generating lattice cross-sections. In this work, we focus on
2D extruded lattices with constant cross-sections. Such lattices can be
uniquely defined by their cross-sections and heights. We further require that
the unit cell of the periodic lattice fits within a square bounding box. In
this case, a 2D tessellation can be generated by repeating the unit cell along
the X and Y directions. With this simplification, one can formulate three
types of geometric descriptors:

  * •

Vertex. The vertices of the unit cell bounding box can take on different
styles: it can be a vertex or the vertex can be replaced by a straight edge or
an arc, see Fig. 1.

(a)

(b)

(c)

Figure 1:  LABEL:sub@fig:vs1 The vertex of the unit cell remains as-is.
LABEL:sub@fig:vs2 The vertex is replaced by a straight edge. LABEL:sub@fig:vs3
The vertex is replaced by an arc.

  * •

Edge segment. Each edge of the unit cell bounding box can be divided into
different edge segments, and each can take on a different style: it can be a
straight edge or two arcs (inspired by the bio-inspired honeycomb studied in
the work of Ha et al. [10]), see Fig. 2.

(a)

(b)

Figure 2:  LABEL:sub@fig:es1 The edge segment of the unit cell remains as-is.
LABEL:sub@fig:es2 The edge segment is replaced by two arcs.

  * •

Interior support. Supports can be added to the interior of the unit cell: it
can be no support, a $+$-shape, or an $X$-shape support. A circle can be added
to the center of the unit cell, see Fig. 3.

(a)

(b)

(c)

(d)

(e)

Figure 3:  LABEL:sub@fig:is1 No interior support is added. LABEL:sub@fig:is2 A
$+$-shape support is added. LABEL:sub@fig:is3 A $X$-shape support is added.
LABEL:sub@fig:is4 A $+$-shape support is added along with a circle.
LABEL:sub@fig:is5 A $X$-shape support is added along with a circle.

It is important to point out that the lattice geometric descriptors are not
unique; other choices are possible to generate a different design space. The
descriptors and their sub-options chosen here aim to generate lattices that
bear similar design cues to some traditional or bio-inspired lattices seen in
the literature [10, 38, 39].

A unique lattice design can then be generated in a combinatorial manner by
forming different combinations of the geometric descriptors. Some geometric
descriptors also allow for sub-design options. In the straight edge vertex
case (Fig. 1b), the edge can take on either edge style depicted in Fig. 2. In
the arc vertex case (Fig. 1c), the arc can point inward or outward to the unit
cell. Vertical and horizontal edges of the unit cell bounding box can be
divided into a different number of edge segments. Still, we limit to a minimum
of 2 and a maximum of 4 divisions for all edges. The lattice design created by
a combination of the design options and their sub-options can be conveniently
denoted by an 8-digit key. Each digit represents a choice on each option/sub-
option. A list of all possible options is presented in Table 1. This key-based
design system contains 660 unique lattice design keys. Besides the discrete
design variables, the framework also features a continuous design variable,
namely the thickness of the lattice walls. In the analysis, we fixed the
number of unit cells to be 4 in the cross-section, forming a 2-by-2
arrangement. Note that different numbers of unit cells and their arrangements
can be used to construct the lattice structure and will affect the mechanical
response of the lattice. Studying the effects of the number of unit cells and
unit cell arrangement is a subject of our future works. A survey of some
periodic lattices generated from the system is presented in Fig. 4.

Table 1: Explanation of the 8-digit design key system Digit | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8
---|---|---|---|---|---|---|---|---
Option |  | Vertex
---
| Arc
---
sub-opt
| # segments
---
horizontal
| # segments
---
vertical
| Edge
---
sub-opt
horizontal
| Edge
---
sub-opt
vertical
Interior |  | Interior
---
sub-opt
| Possible
---
values
|  0 (1a)
---
1 (1b)
2 (1c)
| 0 (in)
---
1 (out)
| 2
---
3
4
| 2
---
3
4
|  0 (2a)
---
1 (2b)
|  0 (2a)
---
1 (2b)
|  0 (3a)
---
1 (3b)
2 (3c)
|  0 (3d)
---
1 (3e)
(a) Design 00231121 |  (b) Design 10331011 |  (c) Design 11220020
---|---|---
(d) Design 11331011 |  (e) Design 20220000 |  (f) Design 21221110
Figure 4: A survey of some lattice designs generated by the key-based system.

###  2.2 Finite element simulations

A Python script was developed to generate cross-sectional sketches in the
finite element (FE) analysis package Abaqus [40] based on given design keys.
The sketches were scaled so that the 2-by-2 lattice touches a square bounding
box with a side length of 20 mm. Three-dimensional shell parts were created in
Abaqus via extrusion to a fixed height of 10 mm. All lattice designs were
discretized with 4-node shell elements with reduced integration (S4R) and
uniform mesh sizes of 0.25 mm (in the X-Y plane) and 0.8 mm (along the height
of the lattice structure) were used.

To accurately capture the response of lattices during high strain-rate impact
loading, which is common in lattice-filled sandwich panels applications, the
Johnson-Cook plasticity and damage models developed by Johnson and Cook [41]
were used in the Abaqus/Explicit dynamic simulations. The material is
Ti-6Al-4V, whose properties follow from those used in the work of Wang and Shi
[42] (see Table 2 therein). As a simple approximation to high strain-rate
impact loading, we sandwiched the lattice between two rigid plates, and the
lattices were subjected to dynamic longitudinal compression. The bottom plate,
where the reaction force was measured, was held fixed, and the top plate
traveled downward with a constant velocity determined by the user-defined
strain rate. All sidewalls were traction-free and were free to deform. All
simulations shared a constant final displacement of 2 mm, corresponding to 20%
nominal compressive strain along the Z-axis. Mass scaling was applied to the
lattice to shorten simulation run time so that a minimum stable time increment
of $2\times 10^{-8}$ s was achieved. The reaction force at the bottom plate,
displacement of the top plate, plastic dissipation, damage dissipation, and
elastic strain energy of the lattice structure were outputs of the FE
simulations. Fig. 5 depicts the FE model assembly and a typical deformed
lattice at the end of dynamic compression.

(a)

(b)

Figure 5: FE model setup and results: LABEL:sub@fig:cae Lattice structure and
two rigid plates set up for dynamic longitudinal compression.
LABEL:sub@fig:peeq A typical deformed lattice at the end of the simulation,
colored by the equivalent plastic strain.

A total of 15000 simulations were conducted on an Intel i7-11800H processor
with 8 cores. Each simulation took about one to two minutes to complete. The
design key, nominal strain rate, and thickness of the shell elements were
sampled using a random number generator. The nominal strain rates were
uniformly sampled on the log scale from the range [$10^{2},10^{5}$] s-1,
representating strain rates during blast loading [43]. Wall thicknesses of the
lattice designs were uniformly sampled from the range [0.25,0.75] mm.

###  2.3 Encoding lattice cross-sectional images

As a general way to contain the geometry information in the lattice cross-
section, we converted the 2D cross-sectional sketch into a 128-by-128 binary
image revealing the skeleton of the cross-section. Note that the image did not
include the shell thickness, so the resulting skeleton is always a single
pixel wide. However, this means that the image matrices are highly sparse, and
a naïve attempt to directly feed the image matrix to the NN as input is likely
wasteful on computer memory.

In this work, we employed a simple autoencoder to extract information in the
lattice cross-sectional image. The encoder consists of two dense layers, with
each having 100 neurons with the rectified linear unit as an activation
function. It transforms a 128-by-128 input image into a 100-by-1 latent
feature vector. The decoder consists of a single dense layer with 1282 neurons
and uses the sigmoid function as activation. The autoencoder was developed and
tested in Keras [44] with a TensorFlow [45] backend. To train and test the
autoencoder, 600 design images were randomly chosen from a total of 660
designs, and an 80-20 split was adopted for training and testing. We used the
mean squared error (MSE) as the loss function and the Adam optimizer [46] to
minimize it. The autoencoder was trained for 80 epochs with a batch size of
50. To gauge the ability of the autoencoder to encode lattice structures never
seen in a training set, an additional testing set consisting only of the 60
lattice designs not included in the training data set, was formed.

###  2.4 Neural network for sequence prediction

####  2.4.1 Input data, data augmentation, and loss function

A sampling of the input space was described in Section 2.2. The output arrays
were extracted from Abaqus and were downsampled from 100 time steps down to
50-time steps to reduce size of the training data set. Two groups of inputs
were generated for the NN model. The first group contains information on the
lattice cross-sections in the form of 100-by-1 encoded latent feature vectors,
which remain constant in time. The second group of inputs has six physics-
informed temporal information arrays, which are:

  1. 1.

Lattice shell thickness111All scalars are expressed in the form of a constant
array with a length of 50, same as the number of time steps..

  2. 2.

Final nominal compression strain1.

  3. 3.

$\rm{log}_{\rm 10}\dot{\bar{\epsilon}}\;^{1}$. The logarithm of strain rate is
used instead of the strain rate since the strain rate dependence in Johnson-
Cook model is logarithmic.

  4. 4.

Nominal compression strain value at each output time point.

  5. 5.

Current time value at each output time point.

  6. 6.

A binary elastic wave indicator is defined as:

$I_{e}=\begin{cases}0,&\quad{}t\leq t_{e}\\\ 1,&\quad{}t>t_{e}\\\
\end{cases},$ (1)

where $t_{e}=\frac{H}{\sqrt{E/\rho}}$ is an estimated time for the elastic
stress wave to travel through the height of the lattice $H$. This indicates
that the reaction force at the rigid bottom plate should remain 0 before the
impact stress wave arrives.

Both groups of inputs were normalized by a standard scaler in Scikit-Learn
[47] prior to training. The scaler was fitted to the training data points to
avoid information leakage [28].

To increase the amount of training data, data augmentation was applied. For
each FE simulation (compressed to a constant 20% strain), twelve final nominal
strains in the range [5%,20%] were randomly sampled, and all inputs and
outputs were linearly interpolated up to the selected final strain level. This
effectively generated data points corresponding to the same strain rate but
different final strain level and increased the total number of input data
points to 180000. These data points were divided into training (600 designs
keys, approximately 61.8%), validation (10.9%), testing on lattice designs
seen in the training set (denoted Test1, 18.2%), and testing on lattice
designs unseen by the training set (denoted Test2, 60 design keys,
approximately 9.1%).

We employed the mean absolute error (MAE) [48] as the loss function, defined
as:

${\rm{MAE}}=\frac{\sum^{N}_{i=1}|\bf{Y}_{i}-\hat{\bf{Y}}_{i}|}{N},$ (2)

where $N,\bf{Y}_{i},\hat{\bf{Y}}_{i}$ denote the number of training data
points, ground-truth outputs, and the NN predictions, respectively. The MSE
was chosen as a metric.

####  2.4.2 Neural network model

In solid mechanics, a recurrent NN model known as the gated recurrent unit
(GRU) model has been widely used to predict sequences [49, 50]. As
investigated by Abueidda et al. [49], although the GRU model is more
computationally expensive than the long short-term memory (LSTM) model and the
temporal convolutional network (TCN) model, it gives most accurate predictions
on energy absorption curves of an elastoplastic material undergoing complex
deformations. Due to the capability of the GRU model to learn from complex
deformation histories, we utilize it to predict the energy absorption of
lattices under dynamic compression.

The developed GRU model consists of three stacked layers, each of 300 GRU
units with hyperbolic tangent (tanh) activation, leading to a model with 1.45
million trainable parameters. The loss function was minimized using an Adam
optimizer [46] with an inverse time decay learning rate schedule and an
initial learning rate of $1\times 10^{-3}$. The model was trained for 150
epochs with a batch size of 600, and training was repeated 10 times to obtain
average training time and model accuracy. The complete data set was reshuffled
and partitioned in each training repetition, as described in Section 2.4.1.
The GRU model was implemented and tested in Keras [44] with a TensorFlow [45]
backend. All training was conducted on an Intel i7-11800H processor with GPU
acceleration on a Nvidia GeForce RTX 3050 GPU.

##  3 Results

###  3.1 Lattice design encoding and reconstruction

Using a latent space dimension of 100, we compare the original lattice designs
to the reconstructed ones generated by the trained autoencoder for six
randomly selected designs, see Fig. 6. To quantitatively describe the
reconstruction quality, the dice similarity coefficient (DSC) [17] is used:

${\rm{DSC}}=\frac{2|\bm{I}\cap\bm{I}_{r}|}{|\bm{I}|+|\bm{I}_{r}|},$ (3)

where $\bm{I}$ and $\bm{I}_{r}$ are the ground-truth and reconstructed binary
skeleton images (1 denotes presence of material and 0 indicates void),
respectively.

Figure 6: Comparison of original and reconstructed lattice design skeletons
for six randomly chosen designs.

###  3.2 Predicting energy outputs

To access the number of input data points required in training to obtain
accurate results, different percentages of the total input data were used to
train the model independently and were tested on two identical testing sets;
the results are shown in Fig. 7a. To obtain a measure of the average
performance of the GRU model, training was repeated ten times using the
partition described in Section 2.4.1. A typical convergence curve showing the
loss and metric during training is shown in Fig. 7b.

(a)

(b)

Figure 7: Convergence plot for GRU model training process:
LABEL:sub@fig:pct_plot Scaled mean squared error when different percentage of
the total data is used in training. LABEL:sub@fig:loss Scaled mean absolute
error evolution during training. Note that the MAE shown here is the MAE
computed on the variables scaled by the standard scaler.

As described in Section 2.4.1, we test our model on two testing sets: (1)
Test1, which was chosen from the 600 lattice designs that were seen by
training, and (2) Test2, which contains 60 lattice designs that were
completely unseen by training. Mean and standard deviation of training time,
relative MAE for each component in each test set, are reported in Table 2.
Relative MAE is presented here since the four output arrays considered have
vastly different scales; thus, we normalize the MAE by the range over time
steps of the respective ground truth 222To avoid erroneous relative MAE
calculations when the range of the ground truth is very close to 0 (i.e., for
cases where the stress wave has not arrived yet), a minimum range of 0.25 N is
enforced for reaction force, and $1\times 10^{-2}$ J for energies. .

Table 2: Mean and standard deviation of training time and relative MAE |  | Train time
---
$[$s$]$
| rMAE,
---
RF1, $[$%$]$
| rMAE,
---
PD2, $[$%$]$
| rMAE,
---
DMD3, $[$%$]$
| rMAE,
---
ELSE4, $[$%$]$
Test1 | 3597.44 (23.15) | 6.26 (0.19) | 0.30 (0.01) | 4.15 (0.16) | 0.95 (0.01)
Test2 | 9.32 (1.08) | 1.04 (0.18) | 6.93 (0.61) | 2.38 (0.26)
| Percent
---
difference
\ | 48.75% | 242.58% | 67.09% | 151.83%
1: Reaction force; 2: Plastic dissipation; 3: Damage dissipation; 4: Elastic
strain energy.

A comparison between ground truths and NN predictions, ranked by the
percentile of MAE for each output array, is presented in Figures 8 and 9. The
median model (one that gives the median overall MAE among the 10 training
repetitions) was used to generate the plots.

Figure 8: Comparison of ground truths and GRU predictions on set Test1
(contains designs seen in training), ranked by percentile of MAE to provide a
representative sampling. Note that the MAE is ranked independently for each of
the four output arrays. Figure 9: Comparison of ground truths and GRU
predictions on set Test2 (contains designs unseen by training), ranked by
percentile of MAE to provide a representative sampling. Note that the MAE is
ranked independently for each of the four output arrays.

###  3.3 Extending to new lattice designs through transfer learning

Comparison of model performance between sets Test1 and Test2 shows how well
the trained GRU model can generalize to previously unseen but similar designs
(all were generated by the key-based system detailed in Section 2.1). It is of
interest to further investigate the performance of the as-trained model on
lattice designs that are generated outside of the key-based design system.

Three new geometries were generated manually, and the comparison of the
lattice design with its autoencoder reconstruction is shown in Fig. 10. The
mean and standard deviation of relative MAE for each output array in each test
set using the as-trained model are listed in Table 3. For easier comparison,
results from Test2 in Table 2 are repeated here.

Figure 10: Comparison of original and reconstructed lattice design skeletons for new lattice geometries. Table 3: Mean and standard deviation of relative MAE on new geometries |  | rMAE,
---
RF1, $[$%$]$
| rMAE,
---
PD2, $[$%$]$
| rMAE,
---
DMD3, $[$%$]$
| rMAE,
---
ELSE4, $[$%$]$
Test2 | 9.32 (1.08) | 1.04 (0.18) | 6.93 (0.61) | 2.38 (0.26)
New geometry 1 | 132.45 (70.19) | 7.05 (0.86) | 27.59 (5.38) | 13.47 (1.03)
New geometry 2 | 338.09 (85.60) | 26.35 (2.05) | 38.95 (7.86) | 35.82 (1.98)
New geometry 3 | 186.81 (63.80) | 17.88 (1.62) | 28.52 (4.99) | 26.00 (1.30)
1: Reaction force; 2: Plastic dissipation; 3: Damage dissipation; 4: Elastic
strain energy.

Besides the as-trained model, it is worth investigating how the model can be
extended to unseen geometries by exposing it to a small amount of new input
data via transfer learning. Transfer learning can be used to extend the
trained model to a new field by transferring its knowledge learned previously
from a similar field [51]. The trained GRU model was retrained for 20 epochs
on a new training set of data points corresponding to new geometries. To
prevent catastrophic forgetting during transfer learning, we add 5000 original
training data points corresponding to geometries in the key-based design
system to the training set, although other more involved methods exist [52,
53]. To gauge the presence and amount of catastrophic forgetting, we test our
data on two testing sets: one consists of 60 lattice design keys unseen during
the training of both the original GRU model and the transferred model, while
the other one contains data points for the three new geometries (that are
unseen by the transfer training). Fifty strain rates and thicknesses were
randomly sampled from the same ranges as described in Section 2.2 for each of
the three new geometries, leading to 150 FE simulations. A similar data
augmentation technique detailed in Section 2.4.1 was used to generate 7500
input data points. Plots comparing the relative MAE in both testing sets are
shown in Fig. 11.

Figure 11: Comparison of relative MAE on different output arrays using
different percentages of the total data points corresponding to the new
geometries. The relative MAE of the as-trained GRU model is plotted in red as
a reference.

A comparison between ground truths and transferred GRU model predictions,
ranked by the percentile of MAE for each output array, is presented in Fig.
12. The model retrained with 40% of the new data points was used to generate
the plot.

Figure 12: Comparison of ground truths and GRU predictions on new geometries,
ranked by percentile of MAE to provide a representative sampling. Note that
the MAE is ranked independently for each of the four output arrays.

##  4 Discussion

Fig. 6 demonstrates the ground truth images and their autoencoder
reconstructions. For six randomly chosen cases, we observe a close resemblance
of the two, which indicates that a 100-by-1 latent vector is sufficient to
encode the input images from within the key-based design system. This confirms
autoencoder as an effective tool to encode image inputs [54, 55]. However, as
Fig. 10 shows, the reconstruction quality drops significantly on images from
outside of the key-based design system. Especially for the third new
geometries shown in Fig. 10, its autoencoder reconstruction resembles the
first geometry, although the ground truths are vastly different. This
confusion in the autoencoder indicates that it is specialized in encoding
images within the key-based system, which all bear some design feature
resemblance due to the combinatoric nature of the framework. One possible way
to enhance the encoding capability of the autoencoder is to expose it to a
more diverse design space. As training for the autoencoder is separated from
the training of the GRU model, this additional training can be done
efficiently without the need for any expensive FE simulations. It is also
hypothesized that the prediction accuracy of the subsequent GRU model will
increase if the autoencoder effectively and accurately extracts important
latent features from the given design.

Two validation sets, Test1 and Test2, were always employed for examining the
model’s ability to extend the structure-property relations learned from a
limited input space to new design geometries. Fig. 7a shows the accuracy of
the two sets when a different number of data points were used in training.
Although loss keeps decreasing on set Test1, model performance plateaus on set
Test2 after about 50% of points are used. For best performance in both sets,
all subsequent models in this work were trained using 61.8% of total data
points, as described in Section 2.4.1. It can be inferred from Fig. 7b that no
severe overfitting occurred during training. Comparing the results in Table 2,
we see that the prediction accuracy on set Test2 has deteriorated by about 50
to 240 percent compared to set Test1, but we note that all relative MAEs in
Test2 remain smaller than 10% of the ground truth range. The statistical
distributions of MAEs are shown in Fig. 8 and Fig. 9. Comparing the first
three columns (up to 75th percentile), we observe that the GRU predictions
closely follow the FE simulation results, except when the ground truth values
are very small (e.g., row 3 column 1 in Fig. 8). Deviations from ground truth
grow in later stages of loading, after about 30 time steps. However, even on
the highest MAE case, the GRU model can predict response curves that bear
similar trends as the ground truth in both testing sets. It is worth noting
that a positivity constraint was not enforced on the energy outputs; thus,
some non-physical negative energy values arise when the ground truths are
small. In most cases, the model can accurately predict the reaction curve and
energy absorption history for various lattice designs, with an observed
accuracy similar to that reported by Yang et al. [28], which employs a
convolutional NN architecture as opposed to the recurrent NN architecture in
this work. As is apparent in Table 3, the trained GRU model has a certain
generalization ability on set Test2, but falls short on new geometries outside
of the key-based design system. It is hypothesized that this limited
generalization ability is due to both geometries in sets Test1 and Test2
coming from the same key-based design system and can be remedied by having a
more diverse training set that covers a wider range of input designs. Previous
studies also suggest alternative ways to improve the generalization ability of
trained models, such as altering the input format [56] or using a neural
network ensemble [57].

It is unlikely that the model needs to predict mechanical response and energy
absorption of a new lattice geometry given only its cross-section. It is often
possible to obtain additional input data from limited testing. The model would
then extend predictions to the new lattice design with different wall
thicknesses and/or loading strain rates. Fig. 11 demonstrates the results of
transfer learning with different percentages of new input data. We see that
transfer learning based on the GRU model trained in Section 3.2 is highly
effective. With as few as 20 training epochs using only 30% of the 150 FE
simulations on the new lattice designs, the prediction accuracy on the energy
absorption curves is within 5% of the model performance on the set Test2.
Relative MAE of the reaction force prediction remains relatively high at 30.5%
with inclusion of 30% new data points, but it is still a dramatic decrease
from over 300% in the as-trained model. Also, note the transferred model
retains similar accuracy on set Test2 to the original GRU model, indicating
that catastrophic forgetting has not occurred. Fig. 12 depicts the MAE
distribution of the transferred model. Predictions up to the 75th percentile
show relatively close agreement with FE simulation results, and large
deviations only emerge after about halfway through the loading. This shows
that transfer learning effectively extends a trained GRU model to new lattice
designs outside of the original, finite training design space with limited new
data.

##  5 Conclusions and future work

In this work, a GRU model was trained to predict the reaction force curve and
energy absorption curves of thin-walled lattices during dynamic longitudinal
compression. A novel combinatorial, key-based system was developed to generate
lattice geometries that share some of the design cues with bio-inspired
honeycomb lattices, which the authors believe is the first time that a
combinatorial framework is used to systematically generate thin-walled lattice
designs. An autoencoder was used to effectively encode the lattice cross-
sectional image to reduce input dimensions. The GRU model is able to
accurately predict the energy responses in unseen geometries from the key-
based design system but falls short on directly extending the predicting
capability to lattice designs outside of the key-based design system.
Nonetheless, when enhanced with data points on the new geometries via transfer
learning, the transferred GRU model reclaims much of the accuracy on new
designs, while maintaining high accuracy over the designs in the key-based
system.

We conclude that the trained GRU model can accurately approximate the
relations between the lattice structure and key performance measurements such
as stress-strain curve and energy absorption as a function of compression
strain for a wide variety of thin-walled lattices under dynamic longitudinal
compression. The ability of the trained model to quickly generate performance
predictions for new lattice designs even on low-end laptop platforms renders
itself a suitable guide in preliminary design stages to quickly sift out
potentially performant designs for more detailed analyses.

In future work, we will focus on leveraging the gradients of the trained GRU
model to inversely generate new designs with optimized specific energy
absorption during high strain rate loading, similar to the work by Chen et al.
[58] and Zheng et al. [59]. Different ways to parameterize the lattice design
space, such as using the adjacency matrix [60], will be explored to enhance
the generalizability and generative ability of future NN models. The effect of
a number of unit cells and unit cell arrangements, as well as the role of
periodic boundary conditions, will also be explored.

## Data availability

The data and source code that support the findings of this study can be found
at: https://github.com/Jasiuk-Research-Group/LatticeResponse_NN_Prediction.

## Conflict of interest

The authors declare that they have no conflict of interest.

## Acknowledgements

We (I. J.) acknowledges the support of the Army Research Office contract (No.
W 911NF-18-2-0067) and the National Science Foundation grant (MOMS-1926353).

## CRediT author contributions

Junyan He: Conceptualization, Methodology, Software, Formal analysis,
Investigation, Data Curation, Writing - Original Draft. Shashank Kushwaha:
Methodology, Software, Formal analysis, Investigation, Writing - Original
Draft. Diab Abueidda: Supervision, Writing - Review & Editing. Iwona Jasiuk:
Supervision, Resources, Writing - Review & Editing, Funding Acquisition.

## References

  * Shin et al. [2008] Kwang Bok Shin, Jae Youl Lee, and Se Hyun Cho.  An experimental study of low-velocity impact responses of sandwich panels for korean low floor bus.  _Composite Structures_ , 84(3):228–240, 2008\.
  * Xie and Zhou [2014] Suchao Xie and Hui Zhou.  Impact characteristics of a composite energy absorbing bearing structure for railway vehicles.  _Composites Part B: Engineering_ , 67:455–463, 2014.
  * Van Paepegem et al. [2014] Wim Van Paepegem, Sivakumar Palanivelu, Joris Degrieck, John Vantomme, Bruno Reymen, Dimitrios Kakogiannis, Danny Van Hemelrijck, and Jan Wastiels.  Blast performance of a sacrificial cladding with composite tubes for protection of civil engineering structures.  _Composites Part B: Engineering_ , 65:131–146, 2014.
  * Codina et al. [2017] Ramón Codina, Daniel Ambrosini, and Fernanda de Borbón.  New sacrificial cladding system for the reduction of blast damage in reinforced concrete structures.  _International Journal of Protective Structures_ , 8(2):221–236, 2017.
  * Qi et al. [2013] Chang Qi, Shu Yang, Li-Jun Yang, Zhi-Yong Wei, and Zhen-Hua Lu.  Blast resistance and multi-objective optimization of aluminum foam-cored sandwich panels.  _Composite Structures_ , 105:45–57, 2013.
  * Gama et al. [2001] Bazle A Gama, Travis A Bogetti, Bruce K Fink, Chin-Jye Yu, T Dennis Claar, Harald H Eifert, and John W Gillespie Jr.  Aluminum foam integral armor: a new dimension in armor design.  _Composite structures_ , 52(3-4):381–395, 2001\.
  * Tarlochan et al. [2007] F Tarlochan, AMS Hamouda, E Mahdi, and BB Sahari.  Composite sandwich structures for crashworthiness applications.  _Proceedings of the Institution of Mechanical Engineers, Part L: Journal of Materials: Design and Applications_ , 221(2):121–130, 2007.
  * Tarlochan [2021] Faris Tarlochan.  Sandwich structures for energy absorption applications: A review.  _Materials_ , 14(16):4731, 2021.
  * Xue and Hutchinson [2006] Zhenyu Xue and John W Hutchinson.  Crush dynamics of square honeycomb sandwich cores.  _International Journal for Numerical Methods in Engineering_ , 65(13):2221–2245, 2006.
  * Ha et al. [2019] Ngoc San Ha, Guoxing Lu, and Xinmei Xiang.  Energy absorption of a bio-inspired honeycomb sandwich panel.  _Journal of Materials Science_ , 54(8):6286–6300, 2019.
  * Qiao and Chen [2016] Jinxiu Qiao and Changqing Chen.  In-plane crushing of a hierarchical honeycomb.  _International Journal of Solids and Structures_ , 85:57–66, 2016.
  * Paz et al. [2014] J Paz, J Díaz, L Romera, and M Costas.  Crushing analysis and multi-objective crashworthiness optimization of gfrp honeycomb-filled energy absorption devices.  _Finite Elements in Analysis and Design_ , 91:30–39, 2014\.
  * Sun et al. [2010] Guangyong Sun, Guangyao Li, Michael Stone, and Qing Li.  A two-stage multi-fidelity optimization procedure for honeycomb-type cellular materials.  _Computational Materials Science_ , 49(3):500–511, 2010.
  * Panda et al. [2018] Biranchi Panda, Marco Leite, Bibhuti Bhusan Biswal, Xiaodong Niu, and Akhil Garg.  Experimental and numerical modelling of mechanical properties of 3d printed honeycomb structures.  _Measurement_ , 116:495–506, 2018.
  * Christensen [2000] R.M. Christensen.  Mechanics of cellular and other low-density materials.  _International Journal of Solids and Structures_ , 37(1):93–104, 2000.  ISSN 0020-7683.  doi: https://doi.org/10.1016/S0020-7683(99)00080-3.  URL https://www.sciencedirect.com/science/article/pii/S0020768399000803.
  * Gao et al. [2018] Jie Gao, Hao Li, Liang Gao, and Mi Xiao.  Topological shape optimization of 3d micro-structured materials using energy-based homogenization method.  _Advances in Engineering Software_ , 116:89–102, 2018.
  * Abueidda et al. [2020] Diab W Abueidda, Seid Koric, and Nahil A Sobh.  Topology optimization of 2d structures with nonlinearities using deep learning.  _Computers & Structures_, 237:106283, 2020.
  * Duddeck et al. [2016] Fabian Duddeck, Stephan Hunkeler, Pablo Lozano, Erich Wehrle, and Duo Zeng.  Topology optimization for crashworthiness of thin-walled structures under axial impact using hybrid cellular automata.  _Structural and Multidisciplinary Optimization_ , 54(3):415–428, 2016.
  * Zeng and Duddeck [2017] Duo Zeng and Fabian Duddeck.  Improved hybrid cellular automata for crashworthiness optimization of thin-walled structures.  _Structural and Multidisciplinary Optimization_ , 56(1):101–115, 2017.
  * Guo et al. [2021] Guikai Guo, Yanfang Zhao, Weihe Su, and Wenjie Zuo.  Topology optimization of thin-walled cross section using moving morphable components approach.  _Structural and Multidisciplinary Optimization_ , 63(5):2159–2176, 2021.
  * Sharafi et al. [2014] Pezhman Sharafi, Lip H Teh, and Muhammad NS Hadi.  Shape optimization of thin-walled steel sections using graph theory and aco algorithm.  _Journal of Constructional Steel Research_ , 101:331–341, 2014.
  * Verma et al. [2020] Chaman Singh Verma, Behzad Rankouhi, and Krishnan Suresh.  A combinatorial approach for constructing lattice structures.  _Journal of Mechanical Design_ , 142(4), 2020.
  * Bastek et al. [2022] Jan-Hendrik Bastek, Siddhant Kumar, Bastian Telgen, Raphaël N Glaesener, and Dennis M Kochmann.  Inverting the structure–property map of truss metamaterials by deep learning.  _Proceedings of the National Academy of Sciences_ , 119(1), 2022.
  * Baykasoğlu et al. [2020] Adil Baykasoğlu, Cengiz Baykasoğlu, and Erhan Cetin.  Multi-objective crashworthiness optimization of lattice structure filled thin-walled tubes.  _Thin-Walled Structures_ , 149:106630, 2020.
  * Wang et al. [2021] Jun Wang, Jesse Callanan, Oladapo Ogunbodede, and Rahul Rai.  Hierarchical combinatorial design and optimization of non-periodic metamaterial structures.  _Additive Manufacturing_ , 37:101710, 2021.
  * Callanan et al. [2018] Jesse Callanan, Oladapo Ogunbodede, Maulikkumar Dhameliya, Jun Wang, and Rahul Rai.  Hierarchical combinatorial design and optimization of quasi-periodic metamaterial structures.  In _International Design Engineering Technical Conferences and Computers and Information in Engineering Conference_ , volume 51760, page V02BT03A011. American Society of Mechanical Engineers, 2018.
  * Wang et al. [2018] ZhiJin Wang, HuaZhi Zhou, VI Khaliulin, and AV Shabalov.  Six-ray folded configurations as the geometric basis of thin-walled elements in engineering structures.  _Thin-Walled Structures_ , 130:435–448, 2018.
  * Yang et al. [2020] Charles Yang, Youngsoo Kim, Seunghwa Ryu, and Grace X Gu.  Prediction of composite microstructure stress-strain curves using convolutional neural networks.  _Materials & Design_, 189:108509, 2020.
  * Gu et al. [2018] Grace X Gu, Chun-Teh Chen, and Markus J Buehler.  De novo composite design based on machine learning algorithm.  _Extreme Mechanics Letters_ , 18:19–28, 2018.
  * Chen and Gu [2019] Chun-Teh Chen and Grace X Gu.  Machine learning for composite materials.  _MRS Communications_ , 9(2):556–566, 2019.
  * Yang et al. [2019] Charles Yang, Youngsoo Kim, Seunghwa Ryu, and Grace X Gu.  Using convolutional neural networks to predict composite properties beyond the elastic limit.  _MRS Communications_ , 9(2):609–617, 2019.
  * Abueidda et al. [2019] Diab W. Abueidda, Mohammad Almasri, Rami Ammourah, Umberto Ravaioli, Iwona M. Jasiuk, and Nahil A. Sobh.  Prediction and optimization of mechanical properties of composites using convolutional neural networks.  _Composite Structures_ , 227:111264, 2019.  ISSN 0263-8223.  doi: https://doi.org/10.1016/j.compstruct.2019.111264.  URL https://www.sciencedirect.com/science/article/pii/S0263822319312383.
  * Laban et al. [2020] Othman Laban, Samer Gowid, and Elsadig Mahdi.  Experimental investigation and uncertainty prediction of the load-carrying capacity of composite double hat for lattice core sandwich panels using artificial neural network.  In _2020 IEEE International Conference on Informatics, IoT, and Enabling Technologies (ICIoT)_ , pages 67–72. IEEE, 2020.
  * Messner [2019] Mark C. Messner.  Convolutional Neural Network Surrogate Models for the Mechanical Properties of Periodic Structures.  _Journal of Mechanical Design_ , 142(2), 10 2019.  ISSN 1050-0472.  doi: 10.1115/1.4045040.  URL https://doi.org/10.1115/1.4045040.  024503\.
  * Garland et al. [2020] Anthony P Garland, Benjamin C White, Bradley H Jared, Michael Heiden, Emily Donahue, and Brad L Boyce.  Deep convolutional neural networks as a rapid screening tool for complex additively manufactured structures.  _Additive Manufacturing_ , 35:101217, 2020.
  * Hassanin et al. [2020] Hany Hassanin, Yusra Alkendi, Mahmoud Elsayed, Khamis Essa, and Yahya Zweiri.  Controlling the properties of additively manufactured cellular structures using machine learning approaches.  _Advanced Engineering Materials_ , 22(3):1901338, 2020.
  * Zok et al. [2016] Frank W Zok, Ryan M Latture, and Matthew R Begley.  Periodic truss structures.  _Journal of the Mechanics and Physics of Solids_ , 96:184–203, 2016.
  * Yang et al. [2018] Xianfeng Yang, Yuxin Sun, Jialing Yang, and Qifan Pan.  Out-of-plane crashworthiness analysis of bio-inspired aluminum honeycomb patterned with horseshoe mesostructure.  _Thin-Walled Structures_ , 125:1–11, 2018.
  * San Ha and Lu [2020] Ngoc San Ha and Guoxing Lu.  A review of recent research on bio-inspired structures and materials for energy absorption applications.  _Composites Part B: Engineering_ , 181:107496, 2020.
  * SIMULIA [2021] SIMULIA.  Abaqus, 2021.
  * Johnson and Cook [1985] Gordon R Johnson and William H Cook.  Fracture characteristics of three metals subjected to various strains, strain rates, temperatures and pressures.  _Engineering Fracture Mechanics_ , 21(1):31–48, 1985.
  * Wang and Shi [2013] Xuemei Wang and Jun Shi.  Validation of johnson-cook plasticity and damage model using impact experiment.  _International Journal of Impact Engineering_ , 60:67–75, 2013.
  * Othman and Marzouk [2016] Hesham Othman and H Marzouk.  Strain rate sensitivity of fiber-reinforced cementitious composites.  _American Concrete Institute Materials Journal_ , 113(2):143–150, 2016.
  * Chollet et al. [2015] Francois Chollet et al.  Keras, 2015.  URL https://github.com/fchollet/keras.
  * Abadi et al. [2015] Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S. Corrado, Andy Davis, Jeffrey Dean, Matthieu Devin, Sanjay Ghemawat, Ian Goodfellow, Andrew Harp, Geoffrey Irving, Michael Isard, Yangqing Jia, Rafal Jozefowicz, Lukasz Kaiser, Manjunath Kudlur, Josh Levenberg, Dandelion Mané, Rajat Monga, Sherry Moore, Derek Murray, Chris Olah, Mike Schuster, Jonathon Shlens, Benoit Steiner, Ilya Sutskever, Kunal Talwar, Paul Tucker, Vincent Vanhoucke, Vijay Vasudevan, Fernanda Viégas, Oriol Vinyals, Pete Warden, Martin Wattenberg, Martin Wicke, Yuan Yu, and Xiaoqiang Zheng.  TensorFlow: Large-scale machine learning on heterogeneous systems, 2015\.  URL https://www.tensorflow.org/.  Software available from tensorflow.org.
  * Kingma and Ba [2014] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Pedregosa et al. [2011] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.  Scikit-learn: Machine learning in Python.  _Journal of Machine Learning Research_ , 12:2825–2830, 2011\.
  * Willmott and Matsuura [2005] Cort J Willmott and Kenji Matsuura.  Advantages of the mean absolute error (mae) over the root mean square error (rmse) in assessing average model performance.  _Climate research_ , 30(1):79–82, 2005.
  * Abueidda et al. [2021] Diab W Abueidda, Seid Koric, Nahil A Sobh, and Huseyin Sehitoglu.  Deep learning for plasticity and thermo-viscoplasticity.  _International Journal of Plasticity_ , 136:102852, 2021\.
  * Frankel et al. [2019] Ari L Frankel, Reese E Jones, Coleman Alleman, and Jeremy A Templeton.  Predicting the mechanical response of oligocrystals with deep learning.  _Computational Materials Science_ , 169:109099, 2019.
  * Weiss et al. [2016] Karl Weiss, Taghi M Khoshgoftaar, and DingDing Wang.  A survey of transfer learning.  _Journal of Big Data_ , 3(1):1–40, 2016.
  * Lee et al. [2017] Sang-Woo Lee, Jin-Hwa Kim, Jaehyun Jun, Jung-Woo Ha, and Byoung-Tak Zhang.  Overcoming catastrophic forgetting by incremental moment matching.  _Advances in neural information processing systems_ , 30, 2017.
  * Chen et al. [2019] Xinyang Chen, Sinan Wang, Bo Fu, Mingsheng Long, and Jianmin Wang.  Catastrophic forgetting meets negative transfer: Batch spectral shrinkage for safe transfer learning.  _Advances in Neural Information Processing Systems_ , 32, 2019.
  * Kim et al. [2021] Yongju Kim, Hyung Keun Park, Jaimyun Jung, Peyman Asghari-Rad, Seungchul Lee, Jin You Kim, Hwan Gyo Jung, and Hyoung Seop Kim.  Exploration of optimal microstructure and mechanical properties in continuous microstructure space using a variational autoencoder.  _Materials & Design_, 202:109544, 2021.  ISSN 0264-1275.  doi: https://doi.org/10.1016/j.matdes.2021.109544.  URL https://www.sciencedirect.com/science/article/pii/S0264127521000976.
  * Jung et al. [2020] Jaimyun Jung, Jae Ik Yoon, Hyung Keun Park, Hyeontae Jo, and Hyoung Seop Kim.  Microstructure design using machine learning generated low dimensional and continuous design space.  _Materialia_ , 11:100690, 2020.  ISSN 2589-1529.  doi: https://doi.org/10.1016/j.mtla.2020.100690.  URL https://www.sciencedirect.com/science/article/pii/S2589152920301071.
  * Zhang et al. [2019] Yiquan Zhang, Bo Peng, Xiaoyi Zhou, Cheng Xiang, and Dalei Wang.  A deep Convolutional Neural Network for topology optimization with strong generalization ability.  _arXiv e-prints_ , art. arXiv:1901.07761, January 2019.
  * Yang et al. [2013] Jing Yang, Xiaoqin Zeng, Shuiming Zhong, and Shengli Wu.  Effective neural network ensemble approach for improving generalization performance.  _IEEE Transactions on Neural Networks and Learning Systems_ , 24(6):878–887, 2013.  doi: 10.1109/TNNLS.2013.2246578.
  * Chen and Gu [2020] Chun-Teh Chen and Grace X Gu.  Generative deep neural networks for inverse materials design using backpropagation and active learning.  _Advanced Science_ , 7(5):1902607, 2020.
  * Zheng et al. [2021] Li Zheng, Siddhant Kumar, and Dennis M. Kochmann.  Data-driven topology optimization of spinodoid metamaterials with seamlessly tunable anisotropy.  _Computer Methods in Applied Mechanics and Engineering_ , 383:113894, 2021.  ISSN 0045-7825.  doi: https://doi.org/10.1016/j.cma.2021.113894.  URL https://www.sciencedirect.com/science/article/pii/S0045782521002310.
  * Li et al. [2022] Shunning Li, Yuanji Liu, Dong Chen, Yi Jiang, Zhiwei Nie, and Feng Pan.  Encoding the atomic structure for machine learning in materials science.  _Wiley Interdisciplinary Reviews: Computational Molecular Science_ , 12(1):e1558, 2022.


2em0em

# Percolation of ‘Civilisation’ in a Homogeneous Isotropic Universe

Allan L. Alinea, Cedrix Jake C. Jadrin

###### Abstract

In this work, we consider the spread of a ‘civilisation’ in an idealised
homogeneous isotropic universe where all the planets of interest are
habitable. Following a framework that goes beyond the usual idea of
percolation, we investigate the behaviour of the number of colonised planets
with time, and the total colonisation time for three types of universes. These
include static, dark energy-dominated, and matter-dominated universes. For all
these types of universes, we find a remarkable fit with the Logistic Growth
Function for the number of colonised planets with time. This is in spite of
the fact that for the matter- and dark-energy dominated universes, the space
itself is expanding. For the total colonisation time, $T$, the case for a dark
energy-dominated universe is marked with divergence beyond the linear regime
characterised by small values of the Hubble parameter, $H$. Not all planets in
a spherical section of this universe can be ‘colonised’ due to the presence of
a shrinking Hubble sphere. In other words, the recession speeds of other
planets go beyond the speed of light making them impossible to reach. On the
other hand, for a matter-dominated universe, while there is an apparent
horizon, the Hubble sphere is growing instead of shrinking. This leads to a
finite total colonisation time that depends on the Hubble parameter
characterising the universe; in particular, we find $T\sim H$ for small $H$
and $T\sim H^{2}$ for large $H$.

keywords: percolation, homogeneous isotropic universe, dark energy-dominated,
matter-dominatedddominated, Hubble horizon, FLRW metric

Institute of Mathematical Sciences and Physics,
University of the Philippines Los Baños,
College, Los Baños, Laguna 4031 Philippines<EMAIL_ADDRESS>

##  1 Introduction

The question of whether we are alone or not in the Universe dates back to
antiquity; perhaps, as early as the times when modern humans started wondering
about the night sky. With the realisation that the sun is simply one of the
seemingly countless stars out there, many of which hosting their own solar
systems [1, 2, 3], we gained a dazzling insight that the probability of
existence of life beyond the Sun’s influence, cannot be zero. In fact,
counting the number of solar systems and planets in the Goldilocks zone, there
could be a multitude of planets in our galaxy and beyond, harnessing life [4,
5]. Pushing this idea further, with about a 14-billion year old Universe based
on the $\Lambda$CDM model [6, 7, 8], life in other planets may not be limited
to plants and wild animals alone, but could include more technologically
advanced civilisations compared to ours.

In spite of this, we have yet to find a definitive proof for the existence of
alien life, much less made a contact with intelligent extraterrestrial beings.
It is possible that advanced alien civilisations, should they exist, prefer
not to interfere with human affairs or those of other intelligent beings, in
general; letting them take their own course of development [9, 10]. Or,
perhaps more reasonably, the density of ‘living’ planets in the Universe is
extremely low [11]. The extremely large distances between these ‘living’
planets and an even more sparse distribution of putative intelligent
civilisations, effectively hinder discovery of alien life and communication
amongst intelligent beings in different solar systems. Possibly, it is also
for this main reason that the colonisation of other planets and the
propagation of civilisations in the network of stars constituting the Universe
is extremely slow. Although it is about 14 billion years old, it may not be
old enough for life to significantly propagate to habitable or terraformable
planets.111Readers interested in the so called Fermi-paradox and colonisation,
their analysis and/or proposed solutions within the framework of an
effectively static section(s) of the Universe and percolation/cellular
automata, may refer to Refs. [12, 13, 14, 15, 16] and the references therein.

In line with this, one may wonder from a computational physics perspective,
given a spherical section of the Universe that could be bounded by a Hubble
horizon222For cosmological perturbations in the early Universe and its
relation to the horizon, readers may refer to Refs. [17, 18, 19, 20].[21], how
long would it take an advanced civilisation to colonise all the habitable or
terraformable planets in it? Coupled with this question, given long enough
time, can an advanced civilisation possibly colonise all habitable or
terraformable planets in a given spherical section of the Universe?
Considering the dynamics of colonisation, we may also ask, with the
progressive spreading of a civilisation, how would the number of occupied
planets behave with time? Given the complexity of the Universe, taking into
account the variation in planet characteristics, conditions of host stars, and
variations in shapes and age of galaxies, amongst others, these and other
related questions can be challenging to fully answer. Nonetheless, if we are
to look for the distant future of humanity, noting that we are most likely
bound to leave the solar system, they are sensible questions to ponder rooted
in our impulse to survive and inclination to know the cosmos.

This study aims to answer the three questions mentioned above, albeit within a
more modest and doable framework333Perhaps a good realistic step beyond the
idealisation in this work is by using self-replicating probes [23, 24] for
space exploration.. The idea is to take small but significant steps by
investigating a simplified system to pave the way for more in-depth and
complex researches in the future. As part of our simplifying assumptions, we
consider a homogeneous and isotropic system, that is, an FLRW universe [6,
21], composed solely of cells representing habitable planets. Then we take a
three-dimensional spherical section of this system and let a ‘civilisation’
start at the centremost cell. From here, it spreads to other cells with a
constant propagation speed. The nearest cells are occupied first, and each
newly occupied cell becomes a new undying source for the propagation of
‘civilisation’ to neighbouring cells; please refer to the following sections
for more details. Given this scenario, we investigate three cases, namely, (a)
static, (b) dark energy-dominated, and (c) matter-dominated universes [25]. In
contrast to the first one, the last two are both dynamically expanding
universes with different expansion rates dictated by their matter-energy
content.

At a glance, this may look like the propagation of life or civilisation
modeled based on percolation theory [26, 27, 28, 29, 30]. Indeed, it is
similar to those familiar cases involving the spread of forest fire [31] or
disease [32], implemented on a square or cubic lattice. However, the current
study is interesting in that it differs from the usual ‘textbook’ percolation
in at least two ways. First, the cells representing planets are not arranged
following a regular lattice444For percolation on irregular but static lattice,
curious readers may refer to Ref. [33] and references therein.. Rather, they
are located randomly inside a sphere. However, because the locations are
random with each point inside the sphere equally likely to have a cell, the
system is effectively homogeneous and isotropic. Second, apart from the first
case of a static universe, the space itself is expanding555This is one of the
main reasons why it is not straightforward to implement the propagation of
‘civilisation’, as we describe it in this paper, as the usual percolation on a
square or cubic lattice found in textbooks on computational physics [26, 27,
28].. The expansion follows that of an FLRW universe (like our own Universe)
described by the metric given by [6, 21]

$\displaystyle{ds^{2}=-dt^{2}+a^{2}(t)\delta_{ij}dx^{i}dx^{j}},$ (1)

where $t$ and $x^{i}$ are the time and spatial comoving coordinates,
respectively, $\delta_{ij}$ is the Kronecker delta, and $a(t)$ is the scale
factor; here, the Einstein summation convention is followed for the spatial
indices $i,j=1,2,3$. With these two modifications, particularly the second
one, our work may be seen as an effective extension of percolation theory to
the regime involving an irregular and dynamical lattice.

With our sight to fulfill the objectives of this work, we divide the rest of
the paper into four sections. Section 2 involves the case for a static
universe. As we shall see, starting from the centremost cell, ‘civilisation’
propagates to other cells with time following the Logistic Growth Model [34].
In Sec. 3, we deal with an expanding universe driven by a constant dark energy
in Einstein’s general theory of relativity. The way of propagation of
‘civilisation’ is the same as that for the static universe but the space where
the cells are embedded is expanding. We shall find that the unabated expansion
due to dark energy leads to the existence of a horizon preventing the
occupation of all cells in our sphere. Following this is another case of
expanding universe, but this time, it is driven by matter instead of dark
energy. In Sec. 4, we cover a matter-dominated universe with a ‘tamer’
expansion rate. We shall see that although there appears to be a ‘horizon’ as
in Sec. 3, it differs in that the corresponding Hubble sphere grows with time
faster than the propagation of ‘civilisation’. This ultimately results in the
occupation of all cells despite the expansion of the universe. Lastly, we
state our concluding remarks and future prospects in Sec. 5.

##  2 Static Universe

Following a short description about the propagation of civilisation in the
preceding section, we consider an ideal static universe whose only
constituents are habitable planets represented by cells in an empty three
dimensional space. Considering a spherical section of this universe, we have
the following general algorithm for our simulation:

  * •

Randomly distribute $N$ number of cells inside the spherical section following
a uniform probability distribution corresponding to a homogenenous isotropic
universe. In our program, this is done using Permuted Congruential Generator
[35] for the cell coordinates $(x,y,z)$ with the constraint
$x^{2}+y^{2}+z^{2}<r^{2}$, where $r$ is the radius of the spherical section.

  * •

Initially, the centremost cell represent the only ‘living’ planet; all other
cells are yet to be occupied. The time counter for the spread of civilisation
is set to zero and the counter for the number of occupied cells is set to one.
The civilisation from this cell then propagates following a simplified 2D
pattern shown in Fig. 1; the actual simulation involves a 3D space.

  * •

From the centremost cell (1), the civilisation travels at constant speed, $v$,
to the nearest cell (2), where a new civilisation is established; ideally,
within effectively zero amount of time upon arrival as in the usual
percolation of forest fire in a rectangular lattice.

  * •

With two cells (1 and 2) now occupied in the figure, these two civilisations
propagate to their corresponding nearest target cells (3 and 4, respectively).
The main idea is, whenever a cell is occupied, its civilisation spreads to the
nearest cell at the same speed $v$.

  * •

Programmatically, each cell has several flags and state variables indicating
whether they are (a) inhabited or not, (b) targeted or not by a neighbouring
‘civilised’ cell, and (c) measuring the time of travel to a neighbouring cell.
Each time a new ‘uncivilised’ cell is ‘civilised’, the time counter and number
of occupied cells are incremented based on (c) and one, respectively. The
cells are stored in an array of structures with these flags.

  * •

Because neigbouring cells have different distances from the cells targeting
them, they are not reached at the same time. In the figure, cell 2 reaches its
target cell 3 first before cell 1 reaches its target cell 4 because the
distance between cells 2 and 3 is smaller than that for cells 1 and 4. The
propagation of civilisation continues from here following the same logic until
all cells are occupied or ‘civilised’.

Figure 1: Simplified illustration for the occupation of cells or ‘planets’ in
a static universe. The centremost cell is occupied first and ‘life’ or
‘civilisation’ propagates towards neighbouring cells following the order (1)
$\rightarrow$ (2) $\rightarrow$ (3) $\rightarrow$ (4).

Note that in our model, all occupied cells are undying. Once a cell is
occupied, it remains alive and a source of spreading civilisation throughout
the duration of the simulation666We leave it for future work to deal with
death of a civilisation, variation in planet habitability, limited
propagation, etc.. Furthermore, each occupied cell only targets one cell at a
time for the propagation of its civilisation—the nearest unoccupied cell.
Given the highly unlikely scenario where two neighbouring nearest cells are of
(nearly) equal distance from an occupied cell, only one of these two becomes a
target; the target is chosen randomly. For the time of arrival from one
occupied cell located at a position represented by a vector, $\vec{r}_{i}$,
measured from the centremost cell, to the nearest unoccupied (target) cell at
$\vec{r}_{j}$, we have $|\vec{r}_{j}-\vec{r}_{i}|/v$. Certainly, because $v$
is constant, two or more neighbouring occupied cells will, in general, not
arrive simultaneously at their target cells as illustrated in Fig. 1. This
serves to remind us that unlike that of the usual percolation on 2D square
lattice (modeling, for instance, the spread of forest fire,) our cells are not
on a regular lattice. This leads to non-uniform time steps between occupations
of neighbouring cells.

Having established our algorithm above, we perform a simulation involving
$N=5000$ cells located in a sphere with a radius of $L=5.0$ units. For
simplicity, we set the scale factor in the metric given by (1) and the
constant speed of propagation to unity $(v=1)$. Not to be confused with the
speed of light, $c=1$, usually employed in General Relativity and Quantum
Field Theory, the choice of units, $v=1$, simply implies our unit of distance
is the same as the unit of time; e.g., $1\,\text{km}=1\,\text{s}$. Figure 2
shows our simulation result for the number of cells occupied with respect to
time. This is an averaged result over 500 trials. For each trial, 5000 cells
are randomly distributed inside the sphere using the Permuted Congruential
Generator [35] and civilisation is allowed to spread from the centremost cell
with the intention to occupy all cells. With 500 trials, the standard errors
of the mean times to reach cells 2 to 5000, divided by the corresponding mean
times, are all below 0.5%; e.g., the average time to reach cell 5000 is
$\bar{t}=15.58$ units with a standard error of $\sigma_{\bar{t}}=0.019$
corresponding to $\sigma_{\bar{t}}/\bar{t}\approx 0.12\%$. In other words, for
the purposes of this study, 500 trials seem to be more than good enough for
the ‘smooth’ plot shown in Fig. 2.

Figure 2: Behaviour of the number of cells occupied with time in a static
universe, for $N=5000$ cells and 500 trials.

Focusing on the behaviour of the spread of civilisation shown in the figure,
we find a slow start. With a few sources of civilisations, the spread is
correspondingly slow. Then it picks up speed as more and more cells are
occupied corresponding to many sources of civilisations. As is evident in the
figure, the maximum rate of propagation of civilisation happens when around
half of all cells are occupied; about $n=2500$ corresponding to the
approximate time $t\approx 7$ units. Beyond this, while there are more sources
of civilisations, there are fewer cells for them to spread to in our sphere of
interest. The graph of the number of occupied cells, $n$, tapers off with $t$
near $n=N=5000$. Note that increasing the speed, $v$, for the propagation of
civilisation can only increase the number of occupied cells with time. But the
behaviour of $n$ with $t$ remains the same. Informally, the ‘S’-shape curve in
Fig. 2 is simply ‘compressed’ to the left with higher $v$.

Knowing this, our ideal model for the spread of civilisation in a static
universe, appears to follow the Logistic Growth Model [34]. It seems to stand
side-by-side with the other applications of this model such as population
growth [36], spread of communicable disease [37], and chemical reactions in a
closed system [38], amongst others. The prominent behaviour of these dynamical
systems is a relatively slow start in the propagation due to limited sources
(e.g., disease carriers). Then, as the number of sources grows, with many
others yet to be reached by these sources, the propagation speeds up. This
happens until there are significantly more sources than those that can be
infected or occupied in cases of disease spread and colonisation,
respectively. Eventually, the propagation slows down with less number of cells
to occupy (colonisation) or less reactants to react (chemical reactions) until
the elements of the system are exhausted.

Indeed, the behaviour of $n$ with respect to $t$ most likely follows that of
the Logistic Growth Model. Curve fitting our simulation result with the
logistic function given by

$\displaystyle{n(t)=\frac{N_{0}}{1+e^{-k(t-t_{0})}}},$ (2)

where $N_{0},\,k,$ and $t_{0}$ are the adjustable parameters, yields the
equation for the best fit curve indicated in Fig. 2. The corresponding
coefficient of determination, $R^{2}=0.9996$, is extremely close to unity,
indicating an excellent fit.

In the next two sections, we shall see if the Logistic Growth Model is still
followed when the space itself is expanding; that is, homogeneously and
isotropically for dark energy-dominated and matter-dominated universes.

##  3 Dark Energy-Dominated Universe

Let us consider a dark energy-dominated universe characterised by a constant
Hubble parameter, $H\equiv\dot{a}/a$, where the dot indicates derivative with
respect to time. Such a universe corresponds to an exponentially expanding
spacetime driven by a constant dark energy, $\Lambda$, in Einstein’s general
theory of relativity [6, 21]. Indeed, with $H=\text{const.}$, we see that

$\displaystyle{\frac{da}{a}=H\,dt\quad\Rightarrow\quad a\sim e^{Ht}}.$ (3)

Situating (idealised) habitable planets or cells in this universe for
colonisation, as in the immediately preceding section, implies that the
‘physical’ distances between the cells increase with time, even in the absence
of individual random motion. It is akin to two dots on the surface of a rubber
balloon ‘moving’ farther apart as the balloon is inflated. As such, the time
it would take for one civilisation to reach a nearby cell and colonise it, is
larger compared to that in a static universe, given the same propagation
speed, $v$.

In particular, let us consider two cells located at $\vec{r}_{i}$ and
$\vec{r}_{j}$, respectively, measured from the centremost cell. Noting that
our space is coordinatised by comoving coordinates indicated in (1), the time,
$t_{j}$, it takes to reach cell $j$ from cell $i$ starting at time, $t_{i}$,
is given by

$\displaystyle{|\vec{r}_{j}-\vec{r}_{i}|=\int_{t_{i}}^{t_{j}}dt\
\frac{v}{a}}.$ (4)

Needless to say, if $a$ is unity, it reduces to that of the case for static
universe in the immediately preceding section; i.e.,
$|\vec{r}_{j}-\vec{r}_{i}|=v(t_{j}-t_{i})$. Using the relation for the scale
factor given by (3) above, we find upon integration,

$\displaystyle{t_{j}=-\frac{1}{H}\ln\left(e^{-Ht_{i}}-\frac{H}{v}|\vec{r}_{j}-\vec{r}_{i}|\right)}.$
(5)

If both $|\vec{r}_{j}-\vec{r}_{i}|$ and $v$ are set to unity, for instance,
then for777We are omitting here explicit mention of units for brevity. As in
the previous section, $v$ is set to unity, leading to equality in distance and
time units. $H=0.1$ and $t_{i}=0$, we get $t_{j}\approx 1.05$ as opposed to
the case of static universe where $t_{j}=1.00$.

Figure 3: Behaviour of the number of cells occupied with time in a dark
energy-dominated universe, for varying values of the Hubble parameter;
$N=5000$ cells and 500 trials.

Figure 3 shows our simulation result for the number of colonised cells with
time in a dark energy-dominated universe characterised by a constant Hubble
parameter. Here, $N=5000$ cells and the number of trials is 500 for each
series111As in the case for static universe, with 500 trials, the standard
errors of the mean times to reach cells 2 to 5000, divided by the
corresponding mean times, are all below 0.5%; e.g., the average time to reach
cell 5000 is $\bar{t}=30.36$ units with a standard error of
$\sigma_{\bar{t}}=0.089$ corresponding to $\sigma_{\bar{t}}/\bar{t}\approx
0.29\%$., as in Sec. 2. With reference to the static universe ($H=0.00$), we
see a deviation to the right of the number of colonised cells with time. In
other words, it takes longer time to colonise a given number of cells from the
centremost cell due to the expansion of the Universe, and this time increases
with increasing $H$. Nonetheless, looking at the figure, the behaviour of $n$
with $t$ seems to still follow the Logistic Growth Model as in the static
universe. In fact, fitting this function to the data points, we find that the
coefficients of determination remain extremely closed to unity as we vary $H$
above from 0.00 to 0.06. For $H=0.06$ (see Fig. 4), where the deviation is
expected to be the highest among the three (non-static universe) series in
Fig. 3, we find $R^{2}=0.9992\approx 1$, suggesting an excellent fit.

Figure 4: Behaviour of the number of cells occupied with time in a dark
energy-dominated universe, for $H=0.060,\,N=5000$ cells, and 500 trials.
Figure 5: Behaviour of time to occupy all cells with respect to the Hubble
parameter, in a dark energy-dominated universe, for $N=5000$ cells and 500
trials.

Beyond $H=0.06$, however, we see an effective breakdown from the Logistic
Growth Model. This is because for larger values of $H$, all cells in our
sphere can no longer be colonised. This behaviour is apparent in Fig. 5
showing the time to occupy all cells, $T$, with respect to the Hubble
parameter222To make sense of the data here on a galactic scale, consider a
sphere the size of the Milky Way galaxy—radius of about $5\times
10^{20}\,\text{m}$—with only 5000 habitable planets as in our simulation. For
the expansion rate of $H=0.01$ with $v$ corresponding to one percent the speed
light, the total colonisation time of 17 units in the figure simply converts
to about 17 million years; i.e., 1 time unit $\approx$ 1 million years.. As
the expansion rate increases, $T$ correspondingly increases; that is, in a
seemingly linear manner for small $H$ but goes up fast beyond the linear
regime. Above $H=0.06$, there seems to be an asymptote where $T$ diverges. Our
curve fitting function for $T$ with parameters $t_{0}$ and $d$, borrowed from
(5), captures this behaviour with $R^{2}=0.9990$:

$\displaystyle{T=-\frac{1}{H}\ln\big{(}e^{-Ht_{0}}-Hd\big{)}}.$ (6)

For $H\approx 0$, its Taylor series expansion yields $T\sim H$ confirming the
mentioned linear behaviour. Beyond this is an asymptote given by
$e^{-Ht_{0}}-Hd$. For the results in Fig. 5, we have an effective asymptote at
$H=0.06367$.

The existence of a horizon is apparent in (5) describing the time, $t_{j}$, it
takes to reach cell $j$ from cell $i$ starting at time, $t_{i}$. Given a fixed
value of $H$, we find $t_{j}\rightarrow\infty$ as the argument of the
(natural) logarithmic function goes to zero; that is,

$\displaystyle{e^{-Ht_{i}}\rightarrow\frac{H}{v}|\vec{r}_{j}-\vec{r}_{i}|}.$
(7)

Early on, civilisation can propagate from the centremost cell to the
neighbouring cells because the inter-cell physical distances are still small.
But as time goes by, the universe expands, and these inter-cell distances
increase. Colonisation of new cells takes longer and longer until such time
that the next neighbouring cell to be colonised goes beyond the horizon; it
requires a propagation velocity higher than $v$ to be reached.

Certainly, one may consider a greater propagation velocity for the
civilisation to spread to all cells. Following this line of thinking, however,
will only increase the value of $H$ in Fig. 5 at which the total colonisation
time diverges. To see this, we recall the idea of a Hubble horizon [6, 21]
corresponding to a Hubble sphere with comoving radius defined as

$\displaystyle{R_{H}\equiv\frac{c}{aH}},$ (8)

where $c$ is the maximal causal velocity that is physically attainable; i.e.,
the velocity of light. In accord with Hubble’s law, the effective recession
velocity of the objects situated on the Hubble sphere relative to its centre
is $c$. It follows that the total colonisation time, even in our idealised
scenario where the spread of civilisation is made easy by simple propagation
to neighbouring cells, is bound to diverge, either with increasing $H$ or with
increasing radius of our chosen spherical section of a dark energy-dominated
universe. Furthermore, because of this divergence, the breakdown from the
Logistic Growth Model for the number of occupied cells with time seems
inevitable.

##  4 Matter-Dominated Universe

For a matter-dominated universe, the Hubble parameter is related to the scale
factor as $H^{2}\propto a^{-3}$, by virtue of the Friedmann equation [6, 21,
22]. With the definition given by $H\equiv\dot{a}/a$, we find $a\propto
t^{\frac{2}{3}}$. Following the same logic as in the immediately preceding
section, given this relationship between $a$ and $t$, we obtain the time to
reach cell $j$ from cell $i$ given by

$\displaystyle{t_{j}=\frac{2}{3H_{0}}\bigg{[}(1+\frac{3}{2}H_{0}t_{i})^{\frac{1}{3}}+\frac{H_{0}}{2v}|\vec{r}_{j}-\vec{r}_{i}|\bigg{]}^{3}-\frac{2}{3H_{0}}},$
(9)

where $H_{0}$ is the initial Hubble parameter. Needless to say, $H$ is no
longer constant unlike that of a dark energy-dominated universe. In our
simulation, $H_{0}$ is chosen to coincide with the corresponding constant $H$
for a dark energy-dominated universe, for the purpose of comparison.

Figure 6 shows our simulation results for the number of occupied cells with
time for different values of $H_{0}$. Similar to that for the dark energy-
dominated universe, the time it takes to occupy a given number of cells is
higher for higher initial value of the Hubble parameter. Furthermore, a short
visual inspection tells us that the graphs seem to follow the same pattern;
that is, a logistic growth function. Indeed, the coefficients of determination
in the curve fitting using this function, remain very close to unity for all
the series in the figure; for the highest $H_{0}$ in the graph,
$R^{2}=0.9994$, indicating an excellent fit. Having said this, we notice that
on average, the time it takes to occupy a certain number of cells is lower
compared to that of the dark energy-dominated universe. For instance, the
times to reach 4500 cells for $H=0.06$ are 11.5 and 12.6 (time units) for the
matter- and dark energy-dominated universes, respectively. This is justified
by the fact that although both universes start with the same value of the
Hubble parameter, that of the dark energy-dominated universe is constant while
the other one decreases with time; in particular, for the latter, $H\sim 1/t$.
Physically, beyond the initially set time and Hubble parameter, colonisation
is faster for a matter-dominated universe because it expands slower compared
to a dark energy-dominated universe.

Figure 6: Behaviour of the number of cells occupied with time in a matter-
dominated universe, for varying values of the initial Hubble parameter;
$N=5000$ cells and 500 trials.

Consistent with this observation, the time to occupy all cells in our
simulation is correspondingly lower for a matter-dominated universe; see Fig.
7. Furthermore, the figure shows that while $T$ for a dark energy-dominated
universe grows fast in the region $0.0<H\leq 0.06$ and even ‘blows up’ around
$H=0.06$, for a matter-dominated universe, $T$ seems to follow a linear
behaviour with $H_{0}$. In fact, as can be seen in Fig. 8, it crosses the
‘dangerous’ mark about $H=0.06$ for the dark energy-dominated universe without
problem; that is, as far as $H_{0}=1.0$ in our simulation. Our curve fitting
function with parameters $t_{0}$ and $d$, borrowed from (9), captures the
behaviour of $T$ with $H_{0}$ as

$\displaystyle{T=\frac{2}{3H_{0}}\bigg{[}(1+\frac{3}{2}H_{0}t_{0})^{\frac{1}{3}}+\frac{H_{0}d}{2v}\bigg{]}^{3}-\frac{2}{3H_{0}}}.$
(10)

For the result shown in Fig. 8, we have $R^{2}=1.0000$, signifying an
excellent fit. In agreement with our simulation result and the aforementioned
observations, for small values of $H_{0}$, the Taylor series expansion of (10)
tells us that $T\sim H_{0}$. Furthermore, for large values of $H_{0}$, the
time to colonise all cells scales as $H_{0}^{2}$; that is, a dominantly
quadratic behaviour. Lastly, in contrast to (6) for a dark energy-dominated
universe, the relationship given by (10) for a matter-dominated universe, has
no asymptote marking the divergence of $T$.

Figure 7: Comparison of the behaviour of the total times to occupy all cells
with respect to the (initial) Hubble parameter, in a matter-dominated universe
and dark energy dominated universe, for $N=5000$ cells and 500 trials. Figure
8: Behaviour of the time to occupy all cells with respect to the Hubble
parameter, in a matter-dominated universe, for $N=5000$ cells and 500 trials.

In spite of these observations, one may think more deeply if the time for a
civilisation to propagate to all cells can remain finite given a finite number
of cells but with large finite Hubble parameter. In hindsight, following the
same logic as in the previous section, the fact that a matter-dominated
universe also expands homogeneously and isotropically seems to point to the
existence of a Hubble horizon. Beyond the Hubble sphere, the cells recede from
the center at a velocity larger than that of light in accord with Hubble’s
law. Equation (8) for the Hubble horizon, applied for a matter-dominated
universe, seems to indicate an apparent horizon that could give rise to a
divergent $T$; indeed, a short inspection of (8) tells us that
$R_{H}\neq\infty$. However, we cut short this line of thinking for a divergent
$T$ by recalling that for a matter-dominated universe, $a\propto
t^{\frac{2}{3}}$. This certainly leads to a finite $R_{H}\sim t^{\frac{1}{3}}$
but the Hubble sphere is growing. This is in contrast to the case of dark
energy-dominated universe where the comoving Hubble sphere is shrinking. Now,
multiplying $R_{H}$ by $a$ to get the ‘proper’ Hubble radius and taking the
time derivative of the result yields $c=(aR_{H})\dot{}$. We find that although
there is an apparent Hubble horizon, the expansion rate of the Hubble sphere
relative to its centre is equal to $c$! Consequently, given enough time, cells
initially located beyond the horizon eventually enter the Hubble sphere
enabling effectively all cells in a spherical section of a matter-dominated
universe to be colonised.

##  5 Concluding Remarks

In the quest to understand the colonisation potential of an advanced
civilisation across the vast expanse of habitable planets in our universe, we
ponder in this work interesting questions from the perspective of
computational physics. These questions involve the behaviour of the number of
planets colonised with time, and the total colonisation time with respect to
the Hubble parameter describing the universe. Inspired by Percolation Theory,
we investigate these questions within an idealised scenario involving
homogeneous isotropic universe within which habitable planets are embedded. We
study three types of universes, namely, static, dark energy-dominated, and
matter-dominated.

We find that the growth in the number of colonised planets with time in a
spherical section of a static universe follows the Logistic Growth Model. Such
a behaviour confirms our expectation based on its other applications such as
the propagation of fire, spread of disease, and chemical reaction. In
contrast, the same might not be anticipated for a dark energy- and matter-
dominated universes because they are characterised by an expanding space.
Surprisingly, our results indicate that even with this expansion, the
behaviour of the number of colonised planets with time, fits well with the
logistic growth function, albeit with slower occupation or colonisation rate
and a caveat for a dark energy-dominated universe.

The case for a dark energy-dominated universe is characterised by a divergent
total colonisation time. While it behaves linearly for low values of the
Hubble parameter describing the universe, beyond some cut-off, it blows up
marking a breakdown from the Logistic Growth Model. The underlying reason is
the existence of a Hubble horizon corresponding to a shrinking Hubble sphere;
planets beyond this ‘moves’ faster than the speed of light preventing further
colonisation. For a matter-dominated universe, the total colonisation time
behaves linearly for small values of the Hubble parameter and quadratically
for large values of this parameter. Our simulation results suggest that the
total colonisation for a spherical section of this universe remains finite for
arbitrarily large values of the initial Hubble parameter. This is in spite of
the existence of an apparent Hubble horizon. We reason that while there is a
finite Hubble horizon for a matter-dominated universe, it is growing instead
of shrinking as in the case of a dark energy dominated universe. Since its
growth is faster than the rate of colonisation, the colonisation of all
planets in a spherical section of this universe seems to be always possible.

In spite of the limited ideal framework of this study, we find interesting
results about colonisation of planets in the universe from the perspective of
computational physics. However, our humble beginning leaves a lot of avenues
for future exploration. This may include factors involving (a) planet
habitability, (b) death or survival rate of civilisations, and (c) multiple
starting civilisations, amongst others. The first two are certainly related
and the addition of the third one requires more programming resources.
Focusing on the death rate alone for simplicity, we speculate that for a
constant probability of death upon arrival to previously uninhabited planets,
the first order behaviour, $T\sim H$, for small $H$, would remain the same,
albeit $n$ would grow slower with $t$, for both matter- and dark-energy
dominated universes333If the constant death rate corresponds to effectively
lowering the propagation velocity of civilisation in the curve fitting
functions for $T$ with respect to $H$ for both matter- and dark energy-
dominated universes, then the first order behaviour, $T\sim H$, for small $H$,
remains the same based on (6) and (10).. Spherical sections of both static and
matter-dominated universes would continue to be fully occupied given enough
time. However, for a dark-energy dominated universe, the breakdown from the
Logistic Growth Model for $n$ with respect to $t$ would be earlier given the
slower effective propagation rate of civilisation with a shrinking Hubble
sphere. Holding on tight with our effective extension of Percolation Theory,
we hope to gain deeper insight with the inclusion of these factors in the
study of the spread of civilisation in a homogeneously and isotropically
expanding Universe.

## Acknowledgement

CJC Jadrin would like to express his sincere gratitude to the Department of
Science and Technology Science Education Institute (DOST-SEI) for the
financial support during the conduct of this study.

## References

  * [1] Conselice C J, et al. 2016 The Evolution of Galaxy Number Density at z < 8 and Its Implications ApJ 830 2
  * [2] Lauer T R, et al. 2021 New Horizons Observations of the Cosmic Optical Background ApJ 906 2
  * [3] European Space Agency, gaia archive Accessed May 20, 2023, https://gea.esac.esa.int/archive/
  * [4] Lorenz R 2020 Maunder’s Work on Planetary Habitability in 1913: Early Use of the term ‘Habitable Zone’ and a ‘Drake Equation’ Calculation RNAAS 4 6
  * [5] Petigura E A, Howard A W and Marcy G W 2013 Prevalence of Earth-size planets orbiting Sun-like stars PNAS. 110 48
  * [6] Dodelson S and Schmidt F 2021, Modern Cosmology Elsevier Inc,
  * [7] Aghanim N et al. 2020 Planck 2018 results. VI. Cosmological parameters A&A 641
  * [8] Grøn Ø 2002 A new standard model of the universe Eur. J. Phys. 23 135
  * [9] Brin G D 1983 The Great Silence - The Controversy Concerning Extraterrestrial Intelligent Life Q. J. R. Astron. Soc. 24 283–309
  * [10] Annis J 1999 An astrophysical explanation for the ‘great silence J Br. Interplanet. Soc. 52 19-22
  * [11] Vakoch D 2011 The Drake Equation: Estimating the Prevalence of Extraterrestrial Life through the Ages Cambridge University Press
  * [12] Landis G A 1998 The Fermi Paradox: An Approach Based on Percolation Theory J. Br. Interplanet. Soc. 51 163-166
  * [13] Galera E, Galanti G R, and Kinouchi O 2019 Invasion percolation solves Fermi Paradox but challenges SETI projects Int. J. Astrobiol. 18 316-322
  * [14] Hair T W and Hedman A D 2013 Spatial dispersion of interstellar civilizations: a probabilistic site percolation model in three dimensions Int. J. Astrobiol. 12 45-52
  * [15] Lingam M 2016 Interstellar Travel and Galactic Colonization: Insights from Percolation Theory and the Yule Process Astrobiology 16 418-426
  * [16] Bezsudnov I and Snarskii A 2010 Where is everybody? – wait a moment … new approach to the Fermi paradox ArXiv e-prints : https://arxiv.org/abs/1007.2774
  * [17] Alinea A L, Kubota T, Nakanishi Y and Naylor W 2015 Adiabatic regularisation of power spectra in $k$-inflation JCAP 06
  * [18] Alinea A L 2016 Adiabatic regularization of power spectra in nonminimally coupled chaotic inflation JCAP 10
  * [19] Alinea A L and Kubota T 2018 Adiabatic regularization of the power spectrum in nonminimally coupled general single-field inflation Phys. Rev. D 97
  * [20] Alinea A L and Kubota T 2021 Transformation of primordial cosmological perturbations under the general extended disformal transformation Int. J. Mod. Phys. D 30
  * [21] Weinberg S Cosmology Oxford University Press
  * [22] Uzan J P and Lehoucq R 2001 A dynamical study of the Friedmann equations Eur. J. Phys 22 371
  * [23] Wiley K 2011 The Fermi Paradox, Self-Replicating Probes, and the Interstellar Transportation Bandwidth
  * [24] Borgue O and Hein A M 2021 “Near-term self-replicating probes - A concept design,” Acta Astronaut 187 546-556
  * [25] Rubakov V Aand Gorbunov D S 2017 Introduction to the Theory of the Early Universe: Hot big bang theory World Scientific
  * [26] Giordano N and Nakanishi H 2005 Computational Physics (2nd ed) Pearson
  * [27] Boudreau J F and Swanson E S, Applied Computational Physics Oxford University Press
  * [28] Sahimi M 1994 Application of Percolation Theory Taylor & Francis
  * [29] Basta M, Picciarelli V and Stella R 1994 An introduction to percolation Eur. J. Phys 15
  * [30] Paniagua J M, et al 1997 Percolation theory and the spin-1/2 Ising model for ferromagnets at low temperatures. A didactic application Eur. J. Phys 18
  * [31] Beer T 1990 Percolation Theory and Fire Spread Combustion Science and Technology
  * [32] Mello I, Squillante L, Gomes G and Seridonio A 2021 Epidemics, the Ising-model and percolation theory: A comprehensive review focused on Covid-19 Physica A 573
  * [33] Ren J and Zhang L 2017 Inhomogeneous Site Percolation on an Irregular Bethe Lattice with Random Site Distribution J. Stat. Phys. 168
  * [34] Tsoularis A and Wallace J 2002 Analysis of logistic growth models Math. Biosci. 179 21-55
  * [35] O’Neill M E 2014 PCG: A Family of Simple Fast Space-Efficient Statistically Good Algorithms for Random Number Generation Accessed May 20, 2023 urlhttps://www.pcg-random.org/pdf/hmc-cs-2014-0905.pdf
  * [36] Tsoularis A 2001 Analysis of Logistic Growth Models Res. Lett. Inf. Math. Sci, 2 23-46
  * [37] Cheng B and Wang Y 2020 A logistic model and predictions for the spread of the COVID-19 pandemic AIP Publishing
  * [38] Yin X and Zelenay P, Kinetic Models for the Degradation Mechanisms of PGM-Free ORR Catalysts ECS Transactions 85