# Transverse momentum distribution of charged particles and identified hadrons
in p–Pb collisions at the LHC with ALICE

for the ALICE Collaboration
Centro Studi e Ricerche e Museo Storico della Fisica “Enrico Fermi”, Rome,
Italy
Sezione INFN, Bologna, Italy
E-mail

###### Abstract:

Hadron production has been measured at mid-rapidity by the ALICE experiment at
the LHC in proton-lead (p–Pb) collisions at $\sqrt{s_{\rm NN}}$ = 5.02 TeV.
The transverse momentum ($p_{\rm T}$) distribution of primary charged
particles and of identified light-flavoured hadrons ($\pi^{\pm}$, K±,
K${}^{0}_{\rm S}$, p, $\bar{\rm p}$, $\Lambda$, $\bar{\Lambda}$) are presented
in this report. Charged-particle tracks are reconstructed in the central
barrel over a wide momentum range. Furthermore they can be identified by
exploiting specific energy loss (d$E$/d$x$), time-of-flight and topological
particle-identification techniques. Particle-production yields, spectral
shapes and particle ratios are measured in several multiplicity classes and
are compared with results obtained in Pb–Pb collisions at the LHC.

The measurement of charged-particle transverse momentum spectra and nuclear
modification factor RpPb indicates that the strong suppression of high-$p_{\rm
T}$ hadrons observed in Pb–Pb collisions is not due to initial-state effects,
but it is rather a fingerprint of jet quenching in hot QCD matter. The
systematic study of the hadronic spectral shapes as a function of the particle
mass and of particle ratios as a function of charged-particle density provides
insights into collective phenomena, as observed in Pb–Pb collisions. Similar
features, that could be present in high-multiplicity p–Pb collisions, will
also be discussed.

##  1 Introduction

High-energy heavy-ion (AA) collisions offer a unique possibility to study
hadronic matter under extreme conditions, in particular the deconfined quark-
gluon plasma which has been predicted by quantum chromodynamics (QCD) [1, 2,
3, 4]. The interpretation of the results depends crucially on the comparison
with results from smaller collision systems such as proton-proton (pp) or
proton-nucleus (pA). Proton-nucleus (pA) collisions are intermediate between
proton-proton (pp) and nucleus-nucleus (AA) collisions both in terms of system
size and number of produced particles. Comparing particle production in pp,
pA, and AA reactions is frequently used to separate initial state effects,
connected to the use of nuclear beams or targets, from final state effects,
connected to the presence of hot and dense matter. Moreover, pA collisions
allow for the investigation of fundamental properties of QCD; the $p_{\rm T}$
distributions and yields of particles of different mass at low and
intermediate momenta of $p_{\rm T}$ $\lesssim$ 3 ${\rm GeV}/c$ (where the vast
majority of particles is produced) can provide important information about the
system created in high-energy hadron reactions.

Previous results on identified particle production in pp and Pb–Pb collisions
at the LHC have been reported in [5, 6, 7, 8, 9, 10, 11]. Results on
transverse momentum distribution and nuclear modification factor of charged
particles in p–Pb collisions at $\sqrt{s_{\rm NN}}$ = 5.02 TeV have been
reported in [12]. In this paper we report on the measurement of $\pi^{\pm}$,
K±, ${\rm K}^{0}_{S}$, p($\rm\overline{p}$) and $\Lambda$($\bar{\Lambda}$)
production in p–Pb collisions at a nucleon-nucleon center-of-mass energy
$\sqrt{s_{\rm NN}}$ = 5.02 TeV.

##  2 Sample and Data analysis

The results presented here were obtained from a sample of the data collected
during the LHC p–Pb run at $\sqrt{s_{\rm NN}}$ = 5.02 TeV in the beginning of
2013. Due to the asymmetric beam energies for the proton and lead beams, the
nucleon-nucleon center-of-mass system was moving in the laboratory frame with
a rapidity of $y_{\rm NN}$ = $-0.465$ in the direction of the proton beam. A
detailed description of the ALICE apparatus can be found in [13] and a
description of the data-taking and trigger setup in minimum-bias trigger in
[14]. In order to study the multiplicity dependence, the selected event sample
was divided into seven event classes, based on cuts applied on the total
charge deposited in the VZERO-A scintillator hodoscope ($2.8<\eta_{\rm
lab}<5.1$, Pb beam direction).

The ALICE central-barrel tracking covers the full azimuth within $|\eta_{\rm
lab}|<0.9$. The tracking detectors are located inside a solenoidal magnet
providing a magnetic field of 0.5 T. The innermost barrel detector is the
Inner Tracking System (ITS). The Time Projection Chamber (TPC), the main
central-barrel tracking device, follows outwards. Finally the Transition
Radiation Detector (TRD) extends the tracking farther away from the beam axis.
Charged-hadron identification in the central barrel was performed with the
ITS, TPC [15] and Time-Of-Flight (TOF) [16] detectors [17]. Three approaches
were used for the identification of $\pi^{\pm}$, K±, and p($\bar{\rm p}$),
called “ITS standalone”, “TPC/TOF” and “TOF fits” and are described in details
in [8, 9]. Contamination from secondary particles was subtracted with a data-
driven approach, based on the fit to the transverse distance-of-closest
approach to the primary vertex (DCAxy) distribution with the expected shapes
for primary and secondary particles [8, 9].

The ${\rm K}^{0}_{S}$ and $\Lambda$($\bar{\Lambda}$) particles were identified
exploiting their “${\rm V}^{0}$” weak decay topology in the channels ${\rm
K}^{0}_{S}\to\pi^{+}\pi^{-}$ and
$\Lambda(\bar{\Lambda})\to\rm{p}\pi^{-}(\rm{\bar{p}}\pi^{+})$. The selection
criteria used to define two tracks as ${\rm V}^{0}$ decay candidates are
detailed in [6, 18]. The contribution from weak decays of the charged and
neutral $\Xi$ to the $\Lambda$($\bar{\Lambda}$) yield has been corrected
following a data-driven approach.

The study of systematic uncertainties follows the analysis described in [8, 9,
6, 18] and was repeated for the different multiplicity bins in order to
separate the sources of uncertainty which are dependent on multiplicity and
uncorrelated across different bins (depicted as shaded boxes in the figures).

##  3 Results

Figure 1: Ratios p/$\pi$ (left) and $\Lambda$/K${}_{\rm S}^{0}$ (right) as a
function of $p_{\rm T}$ in two multiplicity bins compared to results in Pb–Pb
collisions. The empty boxes show the total systematic uncertainty; the shaded
boxes indicate the contribution uncorrelated across multiplicity bins (not
estimated in Pb–Pb).

Figure 2: Ratios p/$\pi$ (left) and $\Lambda$/K${}_{\rm S}^{0}$ (right) as a
function of the charged-particle density d$N_{\rm ch}$/d$\eta$ in three
$p_{\rm T}$ intervals in p–Pb, Pb–Pb and pp collisions (pp only shown for
$\Lambda$/K${}_{\rm S}^{0}$). The dashed lines show the corresponding power-
law fit.

Figure 3: Exponent of the p/$\pi$ (left) and $\Lambda$/K${}_{\rm S}^{0}$
(right) power-law fit as a function of $p_{\rm T}$ in p–Pb, Pb–Pb and pp
collisions (pp only shown for $\Lambda$/K${}_{\rm S}^{0}$).

The $p_{\rm T}$ distributions of $\pi^{\pm}$, K±, ${\rm K}^{0}_{S}$,
p($\rm\overline{p}$) and $\Lambda$($\bar{\Lambda}$) in $0<y_{\rm CMS}\ <0.5$
are reported in [19] for different multiplicity intervals.
Particle/antiparticle as well as charged/neutral kaon transverse momentum
distributions are identical within systematic uncertainties. The $p_{\rm T}$
distributions show a clear evolution, becoming harder as the multiplicity
increases. The multiplicity dependence of the $p_{\rm T}$ spectral shape is
stronger for heavier particles, as evident when looking at the ratios ${\rm
K}/\pi$ = (K+\+ K-)/($\pi^{+}$\+ $\pi^{-}$), ${\rm p}/\pi$ = (p +
$\rm\overline{p}$)/($\pi^{+}$\+ $\pi^{-}$) and $\Lambda$/${\rm K}^{0}_{S}$ as
functions of $p_{\rm T}$, shown in Fig. 1 for the 0–5% and 60–80% event
classes. The ratios ${\rm p}/\pi$ and $\Lambda$/${\rm K}^{0}_{S}$ show a
significant enhancement at intermediate $p_{\rm T}$ $\sim 3$ ${\rm GeV}/c$,
qualitatively reminiscent of the one measured in Pb–Pb collisions [8, 9, 18].
The latter is generally discussed in terms of collective flow or quark
recombination [20, 21, 22]. A similar enhancement of the ${\rm p}/\pi$ ratio
in high-multiplicity d–Au collisions has also been reported for RHIC energies
[23].

It is worth noticing that the ratio ${\rm p}/\pi$ as a function of
$\mathrm{d}N_{\mathrm{ch}}/\mathrm{d}\eta$ in a given $p_{\rm T}$-bin follows
a power-law behavior: $\frac{\rm p}{\pi}\left(p_{\rm T}\right)=A(p_{\rm
T})\times\left[\mathrm{d}N_{\mathrm{ch}}/\mathrm{d}\eta\right]^{B(p_{\rm
T})}$. As shown in Fig. 2, the same trend is also observed in Pb–Pb
collisions. The exponent of the power-law function exhibits the same value in
both collision systems (Fig. 3, left). The same feature is also observed in
the $\Lambda$/${\rm K}^{0}_{S}$ ratio and this also holds in pp collisions
(Fig. 3, right).

##  4 Discussion

Figure 4: Results of blast-wave fits, compared to Pb–Pb data, pp data and MC
simulations from PYTHIA8 with and without color reconnection. Charged-particle
multiplicity increases from left to right. Uncertainties from the global fit
are shown as correlation ellipses for p–Pb and Pb–Pb data and with errors bars
for pp data. Figure 5: Pion, kaon, and proton transverse momentum
distributions in the 5-10% multiplicity class compared to the several models
(see text for details).

In heavy-ion collisions, the flattening of transverse momentum distribution
and its mass ordering find their natural explanation in the collective radial
expansion of the system [24]. This picture can be tested in a blast-wave model
[25] with a simultaneous fit to all particles. This parameterization assumes a
locally thermalized medium, expanding collectively with a common velocity
field and undergoing an instantaneous common freeze-out. The fit presented
here is performed in the same range as in [8, 9], also including ${\rm
K}^{0}_{S}$ and $\Lambda$($\bar{\Lambda}$). The results are reported Fig. 4.
Variations of the fit range lead to large shifts ($\sim 10\%$) of the fit
results (correlated across centralities), as discussed for Pb–Pb data in [8,
9]. As can be seen in Fig. 4, the parameters show a similar dependency with
event multiplicity as observed with the Pb–Pb data. Within the limitations of
the blast-wave model, this observation is consistent with the presence of
radial flow in p–Pb collisions. Under the assumptions of a collective
hydrodynamic expansion, a larger radial velocity in p–Pb collisions has been
suggested as a consequence of stronger radial gradients in [26]. On the other
hand it is worth noticing that very similar results are obtained when
performing the same study on pp spectra measured as a function of the event
multiplicity. Other processes not related to hydrodynamic collectivity could
also be responsible for the observed results. This is illustrated in Fig. 4,
which shows the results obtained by applying the same fitting procedure to
transverse momentum distributions from the simulation of pp collisions at
$\sqrt{s}$ = 7 TeV with the PYTHIA8 event generator (tune 4C) [27], a model
not including any collective system expansion. The fit results are shown for
PYTHIA8 simulations performed both with and without the color reconnection
mechanism [28, 29]. With color reconnection the evolution of PYTHIA8
transverse momentum distributions follows a similar trend as the one observed
for p–Pb, pp and Pb–Pb collisions at the LHC, while without color reconnection
it is not as strong. This generator study shows that other final state
mechanisms, such as color reconnection, can mimic the effects of radial flow
[30].

The $p_{\rm T}$ distributions in the 5-10% bin are compared in Fig. 5 with
calculations from the DPMJET [31], Kraków [32] and EPOS LHC 1.99 v3400 [33]
models. The transverse momentum distributions in the 5-10% multiplicity class
are compared to the predictions by Kraków for $11\leq N_{\rm part}\leq 17$,
since the $\mathrm{d}N_{\mathrm{ch}}/\mathrm{d}\eta$ from the model matches
best with the measured value in this class. DPMJET and EPOS events have been
selected according to the charged particle multiplicity in the VZERO-A
acceptance in order to match the experimental selection. DPMJET distributions
are softer than the measured ones and the model overpredicts the production of
all particles for $p_{\rm T}$ lower than about 0.5–0.7 ${\rm GeV}/c$ and
underpredicts it at higher momenta. At high-$p_{\rm T}$, the $p_{\rm T}$
spectra shapes of pions and kaons are rather well reproduced for momenta above
1 and 1.5 ${\rm GeV}/c$ respectively. Final state effects may be needed in
order to reproduce the data. In fact, The Kraków model reproduces reasonably
well the shape of pions and kaons below transverse momenta of 1 ${\rm GeV}/c$
where hydrodynamic effects are expected to dominate. For higher momenta, the
observed deviations for pions and kaons could be explained in a hydrodynamic
framework as due to the onset of a non-thermal component. EPOS can reproduce
the pion and proton distributions within 20% over the full measured range,
while larger deviations are seen for kaons and lambdas. It is interesting to
notice that when final state interactions are disabled in EPOS, the
description of many pp and p–Pb observables worsens significantly [33].

##  5 Conclusions

We presented a comprehensive measurement of $\pi^{\pm}$, K±, ${\rm
K}^{0}_{S}$, p($\rm\overline{p}$) and $\Lambda$($\bar{\Lambda}$) in p–Pb
collisions at $\sqrt{s_{\rm NN}}$ = 5.02 TeV at the LHC. The transverse
momentum distributions show a clear evolution with multiplicity, similar to
the pattern observed in high-energy pp and heavy-ion collisions, where in the
latter case the effect is usually attributed to collective radial expansion.
Models incorporating final state effects give a better description of the
data.

## References

  * [1] N. Cabibbo and G. Parisi,  Phys. Lett. B59, 67 (1975).
  * [2] E. V. Shuryak,  Phys. Lett. B78, 150 (1978).
  * [3] L. D. McLerran and B. Svetitsky,  Phys. Lett. B98, 195 (1981).
  * [4] E. Laermann and O. Philipsen,  Ann. Rev. Nucl. Part. Sci. 53, 163 (2003).
  * [5] ALICE Collaboration, K. Aamodt et al.,  Eur. Phys. J C71, 1655 (2011).
  * [6] ALICE Collaboration, K. Aamodt et al.,  Eur. Phys. J. C71, 1594 (2011).
  * [7] ALICE Collaboration, B. Abelev et al.,  Phys. Lett. B712, 309 (2012).
  * [8] ALICE Collaboration, B. Abelev et al.,  Phys. Rev. Lett. 109, 252301 (2012).
  * [9] ALICE Collaboration, B. Abelev et al.,  (2013), hep-ex/1303.0737.
  * [10] CMS Collaboration, S. Chatrchyan et al.,  Eur. Phys. J. C72, 2164 (2012).
  * [11] CMS Collaboration, V. Khachatryan et al.,  JHEP 1105, 064 (2011).
  * [12] ALICE Collaboration, B. Abelev et al.,  Phys. Rev. Lett. 110, 082302 (2012).
  * [13] ALICE Collaboration, K. Aamodt et al.,  JINST 3, S08002 (2008).
  * [14] ALICE Collaboration, B. Abelev et al.,  Phys. Rev. Lett. 110, 032301 (2013).
  * [15] J. Alme et al.,  Nucl. Instrum. Meth. A622, 316 (2010).
  * [16] A. Akindinov et al.,  Eur. Phys. J. Plus 128, 44 (2013).
  * [17] ALICE Collaboration,  Performance of the ALICE Experiment at CERN LHC,  in preparation.
  * [18] ALICE Collaboration,  (2013), nucl-ex/1307.5530.
  * [19] ALICE Collaboration, B. B. Abelev et al.,  (2013), nucl-ex/1307.6796.
  * [20] R. Fries, B. Muller, C. Nonaka, and S. Bass,  Phys. Rev. Lett. 90, 202303 (2003).
  * [21] P. Bozek,  (2011), nucl-th/1111.4398.
  * [22] B. Muller, J. Schukraft, and B. Wyslouch,  Ann. Rev. Nucl. Part. Sci. 62, 361 (2012).
  * [23] PHENIX Collaboration, A. Adare et al.,  (2013), nucl-ex/1304.3410.
  * [24] U. W. Heinz,  Concepts of heavy ion physics,  CERN-2004-001-D, 2004.
  * [25] E. Schnedermann, J. Sollfrank, and U. W. Heinz,  Phys. Rev. C48, 2462 (1993).
  * [26] E. Shuryak and I. Zahed,  (2013), hep-ph/1301.4470.
  * [27] R. Corke and T. Sjostrand,  JHEP 1103, 032 (2011).
  * [28] P. Z. Skands and D. Wicke,  Eur. Phys. J. C52, 133 (2007).
  * [29] H. Schulz and P. Skands,  Eur. Phys. J. C71, 1644 (2011).
  * [30] A. Ortiz, P. Christiansen, E. Cuautle, I. Maldonado, and G. Paic,  Phys. Rev. Lett. 111, 042001 (2013).
  * [31] S. Roesler, R. Engel, and J. Ranft,  p. 1033 (2000), hep-ph/0012252.
  * [32] P. Bozek,  Phys. Rev. C85, 014911 (2012).
  * [33] T. Pierog, I. Karpenko, J. Katzy, E. Yatsenko, and K. Werner,  (2013), hep-ph/1306.0121.


# Scaling, saturation, and upper bounds in the failure of topologically
interlocked structures

Shai Feldfogel  Institute for Building Materials, ETH Zurich, Switzerland
Konstantinos Karapiperis  Department of Mechanical and Process Engineering,
ETH Zurich, Switzerland  Jose Andrade  Department of Mechanical and Civil
Engineering, Caltech, Pasadena, California, USA  David S. Kammer  Institute
for Building Materials, ETH Zurich, Switzerland

###### Abstract

Topological Interlocking Structures (TIS) have been increasingly studied in
the past two decades. However, some fundamental questions concerning the
effects of Young’s modulus and the friction coefficient on the structural
mechanics of the most common type of TIS application - centrally loaded slabs
- are not yet clear. Here, we present a first-of-its-kind parametric study
based on the Level-Set-Discrete-element-Method that aims to clarify how these
two parameters affect multiple aspects of the behavior and failure of
centrally-loaded TIS slabs. This includes the evolution of the structural
response up to and including failure, the foremost structural response
parameters, and the residual carrying capacity. We find that the structural
response parameters in TIS slabs scale linearly with Young’s modulus, that
they saturate with the friction coefficient, and that the saturated response
provides an upper-bound on the capacity of centrally loaded TIS slabs reported
in the literature. This, together with additional findings, insights, and
observations, comprise a novel contribution to our understanding of the
interlocked structural form.

##  1 Introduction

Topological Interlocking Structures (TIS) are assemblies of un-bonded building
blocks that obtain their structural integrity through the blocks’ interlocking
geometries and through the contact and friction interactions that develop at
their interfaces. The most common and widely studied TIS application is TIS
slabs, under quasi-statically-applied indentation loads [1, 3, 4, 5, 6, 7, 8,
9, 10, 2]. Fig 1(b) illustrates a centrally loaded TIS slab. Fig 1(f-k)
illustrates through sequential snapshots the typical response of such slabs as
observed in experiments and Fig. 1(e) depicts a typical load-deflection curve.

As seen in Fig. 1(f-h), under small loads, the blocks largely stick against
one another and the slab deforms as a whole. We refer to this type of
structural action as stick-governed [11, 12]. As the load increases, slipping
mechanisms initiate near the loaded block and gradually become more
noticeable, see Fig. 1(i,j). Eventually, block/s start to fall-off from the
assembly, and the load drops to zero which mark the structural failure, see
Fig. 1(k). We refer to this type of structural action as slip-governed. The
stick-governed initial response and the slip-governed failure depicted in Fig.
1(f-k), is, by far, the most prevalent one observed in experiments [1, 13, 14,
15, 16, 5, 17]. Therefore, sound understanding of the stick governed regime,
and most importantly, of the slip-governed one which precedes the fall-off of
blocks is key to the study and design of TIS slabs.

Figure 1: Examined configuration - (a) a typical internal block and cross-
sections in the xz and yz planes; (b) examined slab configuration with
boundary and loading conditions. The force P exerted by the indenter on the
slab and the corresponding indenter displacement $\delta$ are indicated by a
yellow arrow in the -z direction; (c) a typical strip showing that $M=6$; (d)
a typical load-deflection curve indicating in red the peak load and the
loading energy; and (f-k) sequential snapshots illustrating the typical
response of TIS slabs to centrally applied loads.

Since TIS slabs transmit loads through contact and friction forces between the
blocks, the main material parameters that govern the interplay between the
stick- and slip-regimes, and thus the behavior and failure of TIS slabs as a
whole, are Young’s modulus $E$ and the interfaces’ ability to resist sliding
[18, 16, 19]. $E$ is important mainly because it governs the slab’s in-plane
compressive stiffness which in turn governs the in-plane forces, the global
stiffness, and the peak load. The interfaces’ ability to resist sliding,
commonly represented by the friction coefficient $\mu$, is important because
it determines the transition to and the evolution of the slip-mechanisms that
eventually lead to blocks’ fall-off and structural failure.

The importance of $E$ and $\mu$ in governing the structural action of TIS
slabs raises fundamental questions that are at focus of the present study: How
does the global structural response parameters, namely the peak load, the
loading energy, and the initial stiffness scale with $E$ and $\mu$? How do $E$
and $\mu$ affect the extent of slip-mechanism-induced damage in terms of the
number of blocks that fall-off and the commensurate residual load carrying
capacity? Is there a saturation level of $E$ and $\mu$ beyond which they do
not further contribute to the carrying capacity or affect the nature of the
failure mechanism? Is there an upper bound to the carrying capacity, and if
so, how is related to $E$ and $\mu$?

As we show next, these fundamental questions have not yet been systematically
addressed in the literature111We attribute this somewhat surprising fact to
the difficulty of capturing and quantifying the slip-governed response of TIS
using the Finite Element Method (FEM) – the most widely-used structural
analysis tool.. The main objective of this paper is to address these open
questions and thereby to provide new fundamental insights into the behavior
and failure of TIS slabs.

Heretofore, the majority of the studies that considered the effects of $E$ and
$\mu$ on the behavior of TIS slabs focused on the stick-governed regime. Most
of these studies employed the Finite Element Method (FEM) [12, 11, 16, 18],
others used an analytical model [19] based on the thrust-line analogy [14, 19,
20], and yet another took a direct experimental approach [11]. These studies
generally concur that the initial global stiffness and the peak load in the
stick-governed regime increase approximately linearly as a function of $E$.
Concerning the effects of $\mu$, it was found experimentally that the peak
load [11] and the initial global stiffness increased while the loading energy
decreased as a function of $\mu$ [12].

By considering a wider range of $\mu$’s than in the aforementioned
experimental studies, Khandelwal et al. [19] found, using FEM analyses, that
the initial global stiffness grows as a function of $\mu$ until, beyond a
certain value of $\mu$, it becomes saturated, meaning it does not further
grow. Dugue et al. [18] observed a similar stiffness saturation. In both these
studies, the phenomenon of saturation beyond a certain $\mu$ was not explored
in terms of the peak load and the loading energy, which are arguably more
significant response parameters, or in terms of the type of failure mechanism.
While the studies that focused on the stick-governed regime provide data on
how $E$ and $\mu$ affect the response parameters on TIS, they do not address
the slip-governed regime which is the one that features much more prominently
in TIS failure. Even within the stick-goverened regime, the treatment is
partial in that it only considered a relatively small ranges of $E$ and $\mu$
and the latter’s effects on only few of the response parameters of interest.

A couple of studies considered the effects of $E$ and $\mu$ in more realistic
context of TIS’ slip-governed failure. Feng et al. [16] studied the slip-
governed response of a TIS slab to a central impact load. They found based on
FEM analyses that all the aforementioned response parameters generally grow
with $\mu$’s in the range 0.1-1. However, more quantitative information on
each of these parameters scales with $\mu$, as well as the effect of $E$ was
not discussed. Also, the more commonly occurring scenario of quasi-static
loading, which is likely to have a different dependence on $E$ and $\mu$, was
not considered.

Koureas et al. [21] studied the effects $E$ and $\mu$ in the more commonly-
considered case of quasi-static loading in the context of a TIS-like beam
configuration. This study showed that the peak load, loading energy, and
initial global stiffness scale linearly with $E$ and increase sub-linearly as
a function of $\mu$ until they saturate for large $\mu$’s. It also showed that
this saturated structural response corresponds to a failure mechanism that is
stick-governed. While providing new observations and interesting clues
relating to the open question that are at the focus of the present study, the
limitation of [21] in the context of the present study is that it did not
consider TIS slabs but a simpler beam configuration. The more complex two-way
structural action of TIS slabs does not allow to apply the beam-based results
in [21] to TIS slabs.

In summary, while some of the effects of $E$ and $\mu$ on the behavior of TIS
slabs have been sporadically addressed in the literature, the existing studies
focused on the pre-failure stick-governed regime, they considered a relatively
small range of $E$ and $\mu$, and they only looked at one of the structural
response parameters at a time. The one study that more fully addressed the
effects of $E$ and $\mu$ within the slip-governed regime was limited to TIS-
like beams and did not address the richer two-way action of TIS slabs, which
are the more common TIS application. The main objectives of this study are to
help bridge this gap and to shed new light on the effects of $E$ and $\mu$ on
the behavior and failure of TIS slabs. To obtain these objective, we a perform
a numerical parametric study using the Level-Set-Discrete-Element-Method (LS-
DEM). A numerical approach allows to isolate the effects of $E$ and $\mu$, as
well as to explore a wide and independent range of these parameters.

##  2 Methodology - an LSDEM approach

LSDEM, originally developed for granular mechanics applications [22], has
recently been adapted for structural analysis of TIS [23], and validated
against experimental results of TIS slabs [24]. It’s main advantage compared
with FEM is that it is well-equipped to capture the slip-governed response of
TIS while closely estimating the load-deflection curve. The concepts
underlying the application of LS-DEM to structural analysis of TIS have been
detailed elsewhere [23], so, we will mention here only the most essential
features.

LS-DEM’s application to TIS involves a seeming contradiction. On the one hand,
TIS are (by definition) infinitely strong and stiff when the blocks’ $E$ is
infinite. On the other hand, our LS-DEM model operates under the assumption of
non-breakable rigid-body blocks. Therefore, on its face, modeling TIS with LS-
DEM should always yield infinite global stiffness and carrying capacity. The
fortunate reason that this is not the case is that the contact between blocks
is enforced in LS-DEM in a penalty sense, which invariably involves small
penetrations between blocks. The cumulative result of these penetrations is
effective global compressibility and commensurate finite global stiffness and
carrying capacity [24].

The main concept underlying our model is to harness this built-in mechanism of
effective global compressibility and to relate it to the true elastic
compressibility of a given TIS slab. Accordingly, we consider the penetration
penalty parameter in our model $k_{n}$ as a proxy of the true $E$ and we
calculate it as:

$k_{n}=E\cdot\frac{M}{L}$ (1)

where $M$ is the number of interfaces in a single row of blocks and $L$ is the
side length of the considered slab. Using a $k_{n}$ determined from Eq. (1),
it has been shown in [24] that the LS-DEM model has the same effective in-
plane stiffness as a TIS slab with an elastic modulus $E$. In what follows, we
shall adhere to Eq. (1) so that $E$ and $k_{n}$ are always directly relatable
by the constant factor $M/L$, which is always known for a given slab. Further
details concerning our LS-DEM model for TIS – additional assumptions, the
rational behind it, and its ability to capture the experimentally observed
behavior and failure of TIS – are available in [24].

##  3 Results

###  3.1 Examined configuration and numerical modeling

The examined configuration follows the experimental set-up in [13] and it is
depicted in Figure 1(a,b). In the slab’s plan view shown in Fig. 1(c), red
line segments indicate the six interfaces in a single row of blocks, so that
$M/L=6/50$ [1/mm]. An illustrative $P-\delta$ curve indicating the main
response parameters we examine is shown in Fig. 1(d).

We chose the set-up described in Figure 1(a,b) for the $E$-$\mu$ parametric
investigation for two reasons. First, square slabs under central indentation
loading is the most common set-up studied in the TIS literature. Second, among
many possible block geometries, planar-faced blocks with an angle of
inclination of 5∘ have been shown to lead to the highest carrying capacities
among a variety of different block geometries [21]. In that sense, they
provide an upper-bound on the block-geometry degree of freedom. Additional
details on the configuration and on the LS-DEM simulations are available in
[24, 23].

Regarding the numerical modeling, the slab was quasi-statically loaded by
prescribing a constant velocity to a 2.5 mm spherical indenter. A rate of 6
mm/sec was found to be low enough so as to not introduce any observable
dynamical effects. For our analyses, we first applied gravity load until the
assembly reached a relaxed state, which we then took as the initial
configuration for the indentation phase. The material density and the friction
coefficient were taken equal to 2500 kg/m3 and $\mu$=0.23, respectively. The
LS-DEM parameters we used for the analyses were 0.6 mm for the surface
discretization and 0.25 mm for the volume discretization. They were chosen
based on the convergence study in ([23]).

###  3.2 Parametric study and discussion

To isolate the effects of $E$ and $\mu$ on the structural response we varied
each of them at a time while keeping the other one constant. We first consider
the effects of $E$ in the range 3.3-75 GPa with a constant $\mu=0.23$ in
3.2.1. Next, we consider the effects of $\mu$ in the range 0.19-9 with a
constant $E=4.2$ GPa in 3.2.2. The constant values of $E$ and $\mu$ are
representative ones for polymers and ceramics typically used to manufacture
TIS blocks [1, 13, 11, 25, 17, 26]. The explored range of the parameters
covers a wide range of materials including concrete (30 GPa) and Aluminum (70
GPa), which have also been used to manufacture TIS blocks [27, 11].

####  3.2.1 The effects of Young’s modulus $E$ \- scaling

Fig. 2 depicts the effects of $E$ on the $P-\delta$ curves and on the examined
response parameters - peak load, initial stiffness, and loading energy. The
structural response is illustrated in Fig. 3 for several $E$’s through
snapshots at constant $\delta$ intervals.

The $P-\delta$ curves for $E\geq 4.2$ GPa exhibits three main response phases,
see Fig. 2a: (I) the load increases linearly up to $\delta\approx 0.4$ mm;
(II) the stiffness gradually decreases until the peak load is reached at
$\delta\approx 1.7-1.8$ mm; and (III) the stiffness becomes negative and
continues to decrease until an abrupt load drop at $\delta=3.6$ mm brings the
load back to zero. Juxtaposing these phases with the evolution of the
structural structural response depicted in Fig. 3, we identify phase (I) with
a mostly stick-governed regime, phase (II) with the initiation and early
evolution of slipping mechanisms at the interfaces of the loaded central
block, and phase (III) with the response becoming slip-governed and the
eventual fall-off of the central block, causing the abrupt load drop at
$\delta$=3.6 mm.

Figure 2: The effects of $E$ on the mechanical behavior and failure of TIS
slabs: (a) Load-deflection curves for various $E$’s (in GPa); (b) peak load
vs. $E$; (c) global initial stiffness vs. $E$; and (d) loading energy vs. $E$.
The dashed and dashed-dotted lines were based on linear interpolation from the
five lowest and highest $E$’s, respectively Figure 3: The effects of $E$ on
the structural response through sequential snapshots: (a) $E=3.3$ GPa; (b)
$E=4.2$ GPa; (c) $E=7.5$ GPa; and (d) $E=75$ GPa.

The similarlity in failure modes for $E\approx 3.7\pm 0.5$ GPa suggests that,
above a threshhold $E$ (which in our case is somewhere between 3.3 and 4.2
GPa), the failure mode converges to a slip-governed one where only the loaded
block ends up falling-off, and where the rest of the structure remains
standing. Furthermore, the higher $E$ is, the less deformed the assembly
remains after the central block falls-off. We note that, although the failure
and the regime that precedes it are slip-governed, the peak loads occur either
entirely within the stick-governed regime or only slightly after the onset of
slip-mechanisms, see the $\delta=1.6$ mm snapshots in Fig. 3.

The behavior in the $E$=3.3 GPa case is somewhat different. Here, the
$P-\delta$ curve is more bell-shaped, it does not exhibit a load drop, and the
maximal deflection is notably smaller at $\delta\approx 2.7$ mm (Fig. 2a).
This corresponds to a response that is entirely stick-governed and it leads to
a collapse-like failure mode where all the blocks end up falling-off (Fig.
3a). The marked difference between the failed configuration for $E\geq 4.2$
GPa and the $E$=3.3 GPa case, together with the one-block-fall-off failure for
$E\geq 4.2$ GPa, speaks to the importance of $E$ in terms of residual carrying
capacity. Here too, the peak loads occurs at $\delta=1.6$ mm, within the
stick-governed regime.

Regarding the effects of $E$ on the response parameters, both the peak load
and the loading energy scale linearly with E, see Fig. 2(b,d). By comparing
Fig. 3(b-d) to Fig. 3(a) at $\delta$=1.6 mm, we also observe that the linear
scaling of the peak load is unaffected by whether slip-mechanisms at the
interfaces of the loaded block have initiated when the peak load is reached or
not. This can be explained through the idea of thrust-lines that transmit
loads through the slab via an internal ”truss” whose depth is the effective
thickness of the slab, that is, the nominal thickness of the slab minus the
deflection under the load. Specifically, the transmission of the load to the
supports through the interfacial contact forces does not depend on whether
there is stick or slip at the interfaces of the loaded block. Rather, it
depends on the magnitude of the deflection of the loaded block, which defines
the slab’s effective thickness. This is evidenced by the similarity in the
distribution of contact forces at the same prescibed displacement across the
range of examined $E$’s, compare Fig. 4(a-c) at common $\delta$’s. We also
note the ratio between the maximal contact force and $E$ is more or less the
same for the entire range of examined $E$’s. The arch-like force-chains
through which the load is transmitted to the supports are shown from two side
views in Fig. 4(d,e).

Figure 4: The evolution of contact forces at the interfaces as a function of
$E$ through sequential snapshots taken normal to the slab’s surface. The
contact forces are represented through cylinders whose axes indicates the
forces’ direction and whose radii indicate their magnitude: (a) $E=3.3$ GPa;
(b) $E=7.5$ GPa; and (c) $E=75$ GPa. Figures (d) and (e) depict the arch-like
internal force-chains through two side views for $E=75$ GPa at $\delta$=0 mm.

The linear scaling of the loading energy is non-trivial and noteworthy because
in ordinary structures and materials larger $E$ are typically correlated with
decreased energy absorption capabilities.

Unlike the peak load and loading energy, the initial global stiffness scales
bilinearly with $E$ (Fig. 2c). The bilinearity kink occurs at the intersection
of the dashed line, which was obtained from linearly interpolating the
smallest five $E$’s, and the dahsed-dotted line, which was obtained based on
the five highest $E$’s. We attribute the bilinearity-kink to slipping that
occurs for $E\geq 20$ GPa at the very first loading stages and that lowered
the initial effective thickness, and, commensurately, the global stiffness of
the slab. That this effect indeed was at play is supported by the fact that,
per a given $\delta$, more slipping occurs for higher $E$, as can be seen by
comparing each of the columns of Fig. 3. Specifically, at $\delta$=0.8, mm we
see that early slipping occurs for $E$=75 GPa but not for $E=7.5<20$ GPa. The
reason why more slipping occurs for higher $E$’s in the first place is that
the ratio $\mu/E$, which represents a normalized resistance to sliding, is
inversely related to $E$.

####  3.2.2 The effects of the friction coefficient $\mu$ \- saturation

Considering the effects of $\mu$ on the $P-\delta$ curves and the response
parameters depicted in Fig. 5(a), the two most noticeable differences compared
with the case of increasing $E$’s are that (1) the pattern of the $P-\delta$
curves changes for larger $\mu$’s, and, most notably (2) the curves eventually
converge both in shape and in value for $\mu\approx 9$. In all cases, the
stiffness gradually decreases until the peak load is reached, and is followed
immediately by an abrupt load drop. However, whereas the curves for $\mu\leq
1$ are bell-like until the first load drop, for $\mu\geq 3$ they gradually
become more linear, until, as $\mu$ approaches 9 at $\delta\approx 3.2$ mm,
the steep load drop occurs. As a results of this change in the shape of the
curves, the peak load occurs at larger $\delta$’s for larger $\mu$’s, unlike
the case with$E\geq 4.2$ GPa, where it is always around $\delta\approx
1.7-1.8$ mm.

Regarding the deflection at which the load-drop occurs for $\mu\geq 3$, it
converges to the thickness of the plate $h=3.18$ mm. This specific value is
explainable in terms of the thrust-line TIS analogy [19] by the vanishing of
the effective thickness and the zero structural depth at $\delta=h$. What is
interesting here is that, contrary to what the thrust-line model predicts, the
slab still maintains a degree of carrying capacity past the load drop at
$\delta=h$, as evidenced by the stair-case shapes observed for $\mu\geq 1$. We
attribute this ability of TIS slabs to sustain loads for $\delta>h$ to the
transmission of loads through slip-mechanisms from central rows/columns of
blocks to ones closer to the boundary, a feature not shared by one-way-bending
TIS beams.

Figure 5: The effects of $\mu$ on the mechanical behavior and failure of TIS
slabs: (a) Load-deflection curves for various $\mu$’s; (b) peak load vs.
$\mu$; (c) global initial stiffness vs. $\mu$; and (d) loading energy vs.
$\mu$. The dashed and dashed-dotted lines were based on linear interpolation
from the five lowest and highest $\mu$’s, respectively Figure 6: The effects
of $\mu$ on the structural response through sequential snapshots: (a)
$\mu=0.19$; (b) $\mu=0.21$; (c) $\mu=1$; and (d) $\mu=9$. The diagonal strips
indicate the onset of slip-mechanisms and the configuration at which the peak
load occurs. The response converges for high $\mu$’s.

Examining the structural response snapshots in Fig. 6, we make the following
observations:

  1. 1.

Similarly to the effect of $E$, there exist a threshhold $\mu=0.2$ that
separates two very different failed configurations. For $\mu\geq 0.21$, the
response converges in the sense that only the central block falls-off while
the rest of the blocks stick and the structures remains standing. For
$\mu<0.2$, the response is mostly stick-governed and the failure is collapse-
like mode where all blocks end up falling-off.

  2. 2.

The global deformations that occur during loading are recovered to a degree
that depends on how high $\mu$ is. In the limit $\mu=9$, the entire global
deformation at $\delta=3.2$ mm is recovered at $\delta=4.8$ mm (Fig. 6d). This
speaks to $\mu$’s importance to the TIS’ residual carrying capacity and
global-deformation-recoverability upon un-loading.

  3. 3.

The larger $\mu$ is, the longer the stick-regime is in terms of $\delta$. This
can be seen in the later onset of slip-mechanism for larger $\mu$’s in the
right diagonal strip in Fig. 6. Recalling that the peak loads occur either
entirely within stick-regime or very early after the initiation of slip-
mechanisms, the prolonged stick-regime for larger $\mu$’s explains the larger
and later-occurring peak load in Fig. 5(a). This is evidenced by the fact that
the $\delta$’s that just precede the onset of slip-mechanisms (shown in the
left diagonal strip in Fig. 6) are very close to the ones where the peak loads
occur for the different $\mu$’s, see Fig. 5(a).

Turning to $\mu$’s effects on the response parameters, we observe that,
following an initial linear increase as a function of $\mu$, all the response
parameter plateau or saturation for large $\mu$’s, see Fig. 5(b-d). Noting in
Fig. 5(a) that the higher peak load value for $\mu=9$ actually occurs before
$\delta=3.2$ mm and within the pre-failure oscillations that occur for the
larger $\mu$’s, the peak loads saturates at $\mu\approx 7$. The global
stiffness saturates at $\mu=2$ and the loading energy saturates at about
$\mu=4-5$.

The small deviations from the saturated loading energy between $\mu=4-9$ are
attributed to the post-peak-load ”staircases” in the $P-\delta$ curves
discussed above in relation to Fig. 5(a). These ”staircases”, which we
attribute to the release of load for central rows of blocks and their
transmission to adjacent ones closer to the boundary, have a minor effect of
the loading energy response parameter and no effect on the other ones.

The saturation of response parameters with $\mu$ has been previously observed
in other contexts: saturation of the initial global stiffness in TIS slabs
[19], saturation of all response parameters in TIS-like beams [21], and
strength saturation in granular systems [28, 29]. Here, we show, for the first
time, that not only the stiffness but all aspects of the structural response
of TIS slabs saturate as a function of $\mu$. This includes the peak load, the
loading energy, the $P-\delta$ curve, the evolution of the structural response
up to and including the failure and post-failure configuration.

In the context of $\mu$ saturation, we note that it occurs at unrealistically
high values of $\mu$. This, however, does not mean that the attainability of a
saturated response is similarly unrealistic. The reason is that the mechanical
meaning of $\mu$ saturation is the maximal possible suppression/delaying of
slip mechanisms. This is obtained through maximization of the interface
resistance to sliding, of which $\mu$ is but a simplistic and rather crude
simplification used for modeling purposes. The ability to increase the
interfacial resistance to sliding through surface architecture, see Djumas
[17], allows to reach the saturated capacities with realistic $\mu$’s [21, 30]

####  3.2.3 $\mu$ saturation = upper bound on carrying capacity

To further explore the significance of the saturated structural response shown
in Fig. 2, we compared the normalized saturated load-deflection curve and the
saturated carrying capacity obtained with $\mu=9$ with all the available
experimental results for the carrying capacity of (non-prestressed) centrally
loaded TIS slabs we found in the literature 7. we found that the normalized
experimental capacities always fall within the envelope that is the normalized
saturated response we found in the present study. Since each of the
experimental references was taken as the maximal experimental capacity
obtained from experimental set-ups reported in that paper, the four reference
points actually represent around 60 experimental slabs made with a wide range
of materials (glass, ceramics, polymers, Alumina-Silicate), slab lineal
dimensions and thicknesses, block shapes (planar and curved-faced), and block
types (dense and hollow). The fact that in such a wide range of experimental
benchmarks, the normalized experimental capacities invariably fall under the
saturated response envelope suggests that the latter provides an upper bound
on the capacity of TIS slabs. This underscores not only the validity of
capacity saturation, but it also speaks to its important physical and
structural significance.

Figure 7: The saturated response is an upper bound on the carrying capacity of
TIS slabs: The experimental peak loads in KH2012[14], KR2012[15], M2016 [31],
and M2019[13] all fall within the envelope of the normalized saturated load-
deflection curve obtained with $\mu=9$. The normalized load $P^{*}$ is taken
as $P/E\cdot L^{2}$ and the normalized deflection $\delta^{*}$ is taken as
$\delta/h$ where $E$,$L$, and $h$ are, respectively, the slabs Young’s
modulus, side length, and thickness.

##  4 Conclusion

We have presented a systematic study on the various effects of the foremost
material properties of Topologically Interlocked Structures (TIS), namely
Young’s modulus $E$ and the friction coefficient $\mu$, on the structural
response of the most common type of TIS application - centrally loaded slabs.

Specifically, we have examined how the load-deflection curve, the structural
response parameters - peak load, initial global stiffness, and loading energy,
the structural response, the type of failure, and the post-failure
configuration depend on $E$ and $\mu$ across a wide and physically
representative range of these two parameters.

We have addressed these fundamental, and previously unaddressed, questions
using a recently developed modeling approach for TIS [24] which is based on
the Level-Set-Discrete-Element-Method, a uniquely-equipped tool to model the
interface-mechanics-governed behavior of TIS slabs.

The main conclusions from our $E-\mu$ parametric study are:

  1. 1.

the peak load and the loading energy scale linearly with $E$, while the
initial global stiffness scales bi-linearly with $E$. The maximal deflection
converges for large $E$’s.

  2. 2.

All response parameters scale initially linearly for small $\mu$’s but later
converge, or saturate, for larger $\mu$’s, as does the load-deflection curve.

  3. 3.

Above threshold values of $E$ and $\mu$, the failure mechanism converges to a
slip-governed one where only the central block falls off and where the
residual global deformations become minimal.

  4. 4.

Below these thresholds, a collapse-like failure where all blocks end up
falling occurs at a relatively small deflection.

  5. 5.

For large $\mu$’s, the response converges to one where a major load drop
occurs when the deflection equals the thickness of the slab. Following this
load drop, the structure still has carrying capacity beyond this deflection
and smaller, staircase-like load drops follow.

  6. 6.

The saturated structural response provides an upper bound on the capacity of
TIS slabs.

##  5 Acknowledgement

Shai Feldfogel was a Swiss Government Excellence Scholarship holder for the
academic years 2021-2022 (ESKAS No. 2021.0165).

## References

  * [1] A.v. Dyskin et al.  “Fracture Resistant Structures Based on Topological Interlocking with Non-planar Contacts”  In _Advanced Engineering Materials_ 5.3, 2003, pp. 116–119  DOI: 10.1002/adem.200390016
  * [2] Yuri Estrin et al.  “Architecturing materials at mesoscale: some current trends”  In _Materials Research Letters_ 9.10, 2021, pp. 399–421  DOI: 10.1080/21663831.2021.1961908
  * [3] A.. Dyskin, Y. Estrin, A.. Kanel-Belov and E. Pasternak  “Topological interlocking of platonic solids: A way to new materials and structures”  In _Philosophical Magazine Letters_ 83.3, 2003, pp. 197–203  DOI: 10.1080/0950083031000065226
  * [4] A.. Dyskin et al.  “The principle of topological interlocking in extraterrestrial construction”  In _Acta Astronautica_ 57.1, 2005, pp. 10–21  DOI: 10.1016/j.actaastro.2004.12.005
  * [5] A. Molotnikov et al.  “Percolation mechanism of failure of a planar assembly of interlocked osteomorphic elements”  In _Engineering Fracture Mechanics_ 74.8, 2007, pp. 1222–1232  DOI: 10.1016/j.engfracmech.2006.07.012
  * [6] M. Carlesso et al.  “Enhancement of sound absorption properties using topologically interlocked elements”  In _Scripta Materialia_ 66.7, 2012, pp. 483–486  DOI: 10.1016/j.scriptamat.2011.12.022
  * [7] Mateus Carlesso et al.  “Improvement of sound absorption and flexural compliance of porous alumina-mullite ceramics by engineering the microstructure and segmentation into topologically interlocked blocks”  In _Journal of the European Ceramic Society_ 33.13, 2013, pp. 2549–2558  DOI: 10.1016/j.jeurceramsoc.2013.05.006
  * [8] A.. Dyskin, Yuri Estrin and E. Pasternak  “Topological Interlocking Materials”  In _Architectured Materials in Nature and Engineering: Archimats_ , Springer Series in Materials Science, 2019, pp. 23–49  URL: https://doi.org/10.1007/978-3-030-11942-3_2
  * [9] Arcady V. Dyskin, Elena Pasternak and Yuri Estrin  “Mortarless structures based on topological interlocking”  In _Frontiers of Structural and Civil Engineering_ 6.2, 2012, pp. 188–197  DOI: 10.1007/s11709-012-0156-8
  * [10] Yuri Estrin, Vinayak R. Krishnamurthy and Ergun Akleman  “Design of architectured materials based on topological and geometrical interlocking”  In _Journal of Materials Research and Technology_ 15, 2021, pp. 1165–1178  DOI: 10.1016/j.jmrt.2021.08.064
  * [11] S. Schaare et al.  “Point loading of assemblies of interlocked cube-shaped elements”  In _International Journal of Engineering Science_ 46.12, 2008, pp. 1228–1238  DOI: 10.1016/j.ijengsci.2008.06.012
  * [12] A. Autruffe et al.  “Indentation Behaviour of Interlocked Structures Made of Ice: Influence of the Friction Coefficient”  In _Advanced Engineering Materials_ 9.8, 2007, pp. 664–666  DOI: 10.1002/adem.200700111
  * [13] Mohammad Mirkhalaf, Amanul Sunesara, Behnam Ashrafi and Francois Barthelat  “Toughness by segmentation: Fabrication, testing and micromechanics of architectured ceramic panels for impact applications”  In _International Journal of Solids and Structures_ 158, 2019, pp. 52–65  DOI: 10.1016/j.ijsolstr.2018.08.025
  * [14] S. Khandelwal, T. Siegmund, R.. Cipra and J.. Bolton  “Transverse loading of cellular topologically interlocked materials”  In _International Journal of Solids and Structures_ 49.18, 2012, pp. 2394–2403  DOI: 10.1016/j.ijsolstr.2012.04.035
  * [15] T. Krause et al.  “Mechanical Properties of Topologically Interlocked Structures with Elements Produced by Freeze Gelation of Ceramic Slurries”  In _Advanced Engineering Materials_ 14.5, 2012, pp. 335–341  DOI: 10.1002/adem.201100244
  * [16] Yuezhong Feng, Thomas Siegmund, Ed Habtour and Jaret Riddick  “Impact mechanics of topologically interlocked material assemblies”  In _International Journal of Impact Engineering_ 75, 2015, pp. 140–149  DOI: 10.1016/j.ijimpeng.2014.08.003
  * [17] Lee Djumas, George P. Simon, Yuri Estrin and Andrey Molotnikov  “Deformation mechanics of non-planar topologically interlocked assemblies with structural hierarchy and varying geometry”  In _Scientific Reports_ 7.1, 2017, pp. 11844  DOI: 10.1038/s41598-017-12147-3
  * [18] M. Dugué, M. Fivel, Y. Bréchet and R. Dendievel  “Indentation of interlocked assemblies: 3D discrete simulations and experiments”  In _Computational Materials Science_ 79, 2013, pp. 591–598  DOI: 10.1016/j.commatsci.2013.07.014
  * [19] S. Khandelwal, T. Siegmund, R.. Cipra and J.. Bolton  “Scaling of the Elastic Behavior of Two-Dimensional Topologically Interlocked Materials Under Transverse Loading”  In _Journal of Applied Mechanics_ 81.3, 2013  DOI: 10.1115/1.4024907
  * [20] M. Short and T. Siegmund  “Scaling, Growth, and Size Effects on the Mechanical Behavior of a Topologically Interlocking Material Based on Tetrahedra Elements”  In _Journal of Applied Mechanics_ 86.11, 2019  DOI: 10.1115/1.4044025
  * [21] Ioannis Koureas, Mohit Pundir, Shai Feldfogel and David S. Kammer  “On the failure of beam-like topologically interlocked structures”  In _International Journal of Solids and Structures_ 259, 2022, pp. 112029  DOI: 10.1016/j.ijsolstr.2022.112029
  * [22] Reid Kawamoto, Edward Andò, Gioacchino Viggiani and José E. Andrade  “Level set discrete element method for three-dimensional computations with triaxial case study”  In _Journal of the Mechanics and Physics of Solids_ 91, 2016, pp. 1–13  DOI: 10.1016/j.jmps.2016.02.021
  * [23] Shai Feldfogel, Konstantinos Karapiperis, Jose Andrade and David S. Kammer  “A discretization-convergent Level-Set-DEM” arXiv:2209.15431 [cs, math]  arXiv, 2022  DOI: 10.48550/arXiv.2209.15431
  * [24] Shai Feldfogel, Konstantinos Karapiperis, Jose Andrade and David S. Kammer  “Failure of topologically interlocked structures – a Level-Set-DEM approach”  arXiv, 2022  DOI: 10.48550/arXiv.2210.14155
  * [25] Ahmed S. Dalaq and Francois Barthelat  “Strength and stability in architectured spine-like segmented structures”  In _International Journal of Solids and Structures_ 171, 2019, pp. 146–157  DOI: 10.1016/j.ijsolstr.2019.04.012
  * [26] Lee Djumas, Andrey Molotnikov, George P. Simon and Yuri Estrin  “Enhanced Mechanical Performance of Bio-Inspired Hybrid Structures Utilising Topological Interlocking Geometry”  In _Scientific Reports_ 6.1, 2016, pp. 26706  DOI: 10.1038/srep26706
  * [27] A. Rezaee Javan et al.  “The impact behaviour of plate-like assemblies made of new interlocking bricks: An experimental study”  In _Materials & Design_ 134, 2017, pp. 361–373  DOI: 10.1016/j.matdes.2017.08.056
  * [28] Konstantinos Karapiperis et al.  “Investigating the incremental behavior of granular materials with the level-set discrete element method”  In _Journal of the Mechanics and Physics of Solids_ 144, 2020, pp. 104103  DOI: 10.1016/j.jmps.2020.104103
  * [29] Stefan Luding  “Anisotropy in cohesive, frictional granular media”  In _Journal of Physics: Condensed Matter_ 17.24, 2005, pp. S2623  DOI: 10.1088/0953-8984/17/24/017
  * [30] Ioannis Koureas, Mohit Pundir, Shai Feldfogel and David S. Kammer  “Beam-like topologically interlocked structures with hierarchical interlocking”, 2022  URL: https://arxiv-export3.library.cornell.edu/abs/2212.06454
  * [31] M. Mirkhalaf, J. Tanguay and F. Barthelat  “Carving 3D architectures within glass: Exploring new strategies to transform the mechanics and performance of materials”  In _Extreme Mechanics Letters_ 7, Mechanics in Extreme Manufacturing, 2016, pp. 104–113  DOI: 10.1016/j.eml.2016.02.016


11institutetext:  School of Computer Science, Carleton University, Ottawa,
Ontario, Canada. 22institutetext: Research supported in part by NSERC
Discovery grant.

# Bike Assisted Evacuation on a Line of Robots with Communication Faults

Khaled Jawhar 11 Evangelos Kranakis 1122

###### Abstract

Two autonomous mobile robots and a non-autonomous one, also called bike, are
placed at the origin of an infinite line. The autonomous robots can travel
with maximum speed $1$. When a robot rides the bike its speed increases to
$v>1$, however only exactly one robot at a time can ride the bike and the bike
is non-autonomous in that it cannot move on its own. An Exit is placed on the
line at an unknown location and at distance $d$ from the origin. The robots
have limited communication behavior; one robot is a sender (denoted by S) in
that it can send information wirelessly at any distance and receive messages
only in F2F (Face-to-Face), while the other robot is a receiver (denoted by R)
in that it can receive information wirelessly but can send information only
F2F. The bike has no communication capabilities of its own. We refer to the
resulting communication model of the ensemble of the two autonomous robots and
the bike as S/R.

Our general goal is to understand the impact of the non-autonomous robot in
assisting the evacuation of the two autonomous faulty robots. Our main
contribution is to provide a new evacuation algorithm that enables both robots
to evacuate from the unknown Exit in the S/R model. We also analyze the
resulting evacuation time as a function of the bike’s speed $v$ and give upper
and lower bounds on the competitive ratio of the resulting algorithm for the
entire range of possible values of $v$.

Key words and phrases. Faults, Line, Robots, Search, Receiver, Sender, S/R
Communication Model.

##  1 Introduction

Evacuation (also known as group search) is similar to search except that it
involves many entities which cooperate in order to find an unknown exit. There
are plenty of applications related to search and evacuation in a distributed
system such as data mining, crawling, and surveillance which makes it an area
of interest. Many researchers have studied the linear search and evacuation
problems with one or more robots moving at different speeds. Most
investigations in this domain highlight algorithms that achieve the best
possible upper and lower bounds. There are several factors that may affect the
solution of the linear search and evacuation problems. For example, the
robot(s) may know how far the exit is from the origin but may not be aware of
the direction (to the left or right of the origin), the robot(s) may know the
direction but not the distance from the origin to the exit, or in a more
complicated case the distance and the direction are unknown to the robots. In
our present study, the distance and the direction are unknown and there is a
communication fault affecting both robots.

###  1.1 Preliminaries, Notation and Terminology

The system consists of two autonomous mobile agents (referred to as robots)
and non-autonomous one (referred to as bike). The two robots and the bike
together are sometimes referred to as the ensemble. The two robots have
different identities and can see each other when they are situated at the same
position. The search domain is the infinite line and it is bidirectional in
that the robots can move in either direction on the line without this
affecting their speeds. The two autonomous mobile agents can move around on
their own with maximum speed $1$ or ride the bike with speed $v>1$. The robots
can change direction at any time as specified by an algorithm; further, this
change in direction of a robot is instantaneous, with or without the bike). An
autonomous robot riding the bike is also referred to as biker, otherwise it
also referred to as hiker.

An evacuation algorithm is given by the complete description of the
trajectories traced by the two bikers until they both find the exit and
evacuate. We are interested in evacuation algorithms which achieve the best
competitive ratio which is defined as the optimal ratio of the evacuation time
achieved by an algorithm which enables both autonomous robots to evacuate from
the exit when its location is unknown divided by the evacuation time needed if
all the robots know where the exit is.

The bike is not autonomous and cannot move and/or communicate on its own and
thus plays only the role of assistant in order to speed up the search. The
autonomous robot using the bike has an advantage in that it can move with
speed $v>1$ which is of course faster than its walking speed $1$. However, in
our model the bike is also a limited resource in that it can be used by only
one robot at a time. There is definitely an advantage in the robots sharing
the bike since this will ultimately improve the overall evacuation time and
hence also the competitive ratio of the ensemble.

An important aspect in our algorithms will be “bike switching”, by which we
mean changing the rider of the bike. We will assume throughout the paper that
bike switching between robots is instantaneous and at no time cost. Note that
the robots may recognize the presence of the bike when they are at the same
location as the bike.

We employ a model referred to as S/R model (this model was first proposed on
the infinite line by [14] and on a unit disk by [18]) in which the robots have
certain limited communication behavior and may communicate throughout the
execution of the algorithms in the following way. F2F (resp. Wireless)
communication can take place only when the robots are co-located (resp. at any
distance). The sender can only send information wirelessly and receive
information Face-To-Face (F2F) while the receiver can only receive information
wirelessly and send information F2F. A typical communication exchange may
involve, e.g., “exit is found”, “bike released”, “switch bike”, etc. Note that
the robots are endowed with pedometers and have computing abilities so that
they can deduce the location of the other robot and/or the bike from relevant
communications exchanged. Moreover, their ability to communicate F2F never
fails. The bike itself has no communication capabilities on its own. In a way
we think of the S/R model as describing a type of fault and in the sequel, we
refer to the ensemble of the two robots and the bike as the S/R
(communication) model.

Throughout the paper we will be using S and R to denote the sender and
receiver respectively (which we often refer to as bikers) and $B$ to denote
the bike (non-autonomous robots). The robots are equipped with pedometers
which can be used either when walking or riding the bike and are identical in
all their capabilities (locomotion). The origin of the real line will be at
the point $x=0$ on the $x$-axis and this will also be the starting location of
the robots and the bike. The adversary may place the exit at either of the
points $\pm d$, where $d>0$ will denote the unknown distance of the exit from
the origin. In addition, $v>1$ will denote the speed that a biker can attain
when riding the bike.

The following lemma is shown in [20] can be proved easily using the standard
analysis of bike-sharing and will be used extensively in the sequel.

###### Lemma 1

Assume two robots of maximum speed $1$ are sharing a bike of speed $v>1$.
Together the ensemble can cover a segment of length $d$ in time
$\frac{d(v+1)}{2v}$. In fact the ensemble travels with speed equal to
$\frac{2v}{v+1}$. Moreover, these time and speed values are optimal.

(Proof of Lemma 1 can be found in the appendix.)

###  1.2 Related work

The first search problems in the literature concerned a single robot and were
proposed for a continuous infinite line by Beck [5] and Bellman [6] with the
focus on stochastic search models and their analysis. Assuming that the
distance and the direction to the exit are unknown, they proposed an optimal
algorithm with competitive ratio $9$. Additional research in that field can be
found in [1, 21]. There have been many variants of the search problem such as
having static or moving targets, involving multiple robots with limited
communication behavior, or robots with differing speeds, etc.

The search and evacuations problems were studied in environments with multiple
distinct speed robots [4, 13, 16]. Evacuation problems have been studied
extensively as well as in various domains including disks, triangles, and
circles [8, 10, 12]. Evacuation algorithms for two robots on a line without
bike assistance are known for two robots without limited communication
behavior and with different speeds $v$ and $1$ such that $v<1$; for example,
there is an evacuation algorithm in the F2F model with optimal competitive
ratio of $\frac{1+3v}{1-v}$, for $v\leq\frac{1}{3}$, and $9$ otherwise, see
[4]. The same paper also considers the wireless communication model. Our work
differs from [4] in two aspects: first we are using a mixed communication
model S/R, and second a bike (non-autonomous robot) is present which is shared
by the two autonomous robots and has the effect of increasing the overall
evacuation time of the ensemble. Additional studies on linear search can be
found in the works of [2, 3] and under various models of linear search
concerning search cost and robot communication in [7, 9, 17].

An interesting variant of the original linear search problem is concerned with
the case where some of the robots are faulty (crash or byzantine). There are
several research works on this theme. The two main papers in this line of
research are [16] for crash-faulty robots and [13] for Byzantine-faulty
robots. Moreover, [15] studies search and evacuation with a near majority of
faulty agents (e.g., three robots at most one of which is byzantine
faulty)–which also exhibits the worst-case evacuation time. In addition, the
S/R mixed communication model (also considered in our current paper) was
considered as a way to model group search in the presence of robots with
faulty communication capabilities and was first introduced for an infinite
line by [14]; it was shown that there is an evacuation algorithm with
competitive ratio of $3+2\sqrt{2}$ and this is optimal; however, this paper
does not include the concept of the non-autonomous robot considered here.

We should also mention that the S/R faulty communication model considered in
our present paper was also studied for a unit disk in the recent paper [18].
Additional research on evacuation from a unit disk in the presence of a
(byzantine) faulty robot was first studied in [11] and further elaborated for
crash and byzantine faults in [19].

Bike-assisted search and evacuation on a line was first considered for two
robots without faults in [20] and soon thereafter [18] for a unit disk. To the
best of our knowledge the present study of Bike-assisted search and evacuation
on a line for two robots and a bike in the S/R model has not been considered
in the search literature.

###  1.3 Outline and results of the paper

In the sequel, in Section 2 we give the main evacuation algorithm in three
parts: Subsubsection 2.1.1 includes the algorithm when $v\in[1,3]$, Subsection
2.1.2 the algorithm when $v\in[3,10]$, and Subsection 2.2 the algorithm when
$v\in[10,\infty)$, while in Subsection 2.3 we plot the evacuation time for the
entire range of the bike’s speed $v$ so that we can gauge the performance of
the main algorithm. Table 1 displays in detail the competitive ratios in each
of the algorithms for the upper bound cases, respectively, based on the speed
$v$ of the biker.

Table 1: Upper bound on the competitive ratio for the three different algorithms; the first column gives the algorithm, the second column the speed $v$ of the biker for which the respective upper bound on the competitive ratio in the third column is valid. Algorithm | Bike Speed | Competitive Ratio
---|---|---
Algorithm 1 | $v\in[1,3]$ | $\frac{2v}{v+1}\left(\frac{2v+\frac{-(7v+v^{2})+v\sqrt{v^{2}+30v+97}}{2v+6}}{v\frac{-(7v+v^{2})+v\sqrt{v^{2}+30v+97}}{2v+6}}\right)$
Algorithm 2 | $v\in[3,10]$ | $\frac{2v}{v+1}\left(1+\frac{1}{v}+\frac{v^{2}-3v-2-\sqrt{v^{4}+18v^{3}-7v^{2}+4v+4}}{4v(1-v)}\right)$
Algorithm 3 | $v\in[10,\infty]$ | $\frac{2v}{v+1}\left(\frac{9}{v}+\frac{1}{2}-\frac{1}{2v^{2}}\right)$

We prove a lower bound in Section 3, and in Table 2 we display the competitive
ratios in each of the algorithms for the lower bound cases, respectively,
based on the interval to which the speed $v$ of the biker belongs to.

Table 2: Lower bound on the competitive ratio for the evacuation of the two robots for the given range on the speed $v$ of the biker. Bike Speed $v$ | Competitive Ratio
---|---
$v\leq 3$ | $\frac{v^{2}+2v-3}{v^{2}-1}$
$v>3$ | $\frac{6}{v+1}$

Finally we conclude in Section 4 with a summary and a discussion of additional
ideas for research.

Due to space limitations all missing proofs can be found in the appendix.

##  2 Evacuation in the S/R Model

In this section we describe the algorithms in order to obtain upper bounds for
two robots in the S/R model. We also analyze the resulting competitive ratios.
Recall that the competitive ratio of an algorithm is defined as the evacuation
time of the algorithm divided by the shortest time required by both robots to
reach the exit if they know where the exit is.

In order to reach the exit in shortest time possible, there are two cases to
consider. In the first case, the two robots choose to move in opposite
direction and the one which finds the exit first will notify the other one
(see Subsection 2.1). In the second case, the sender takes the bike and uses a
Zig-Zag strategy and the receiver imitates the sender using its maximum unit
speed (see Subsection 2.2). Let us study the algorithms proposed in each of
the two cases.

###  2.1 Robots move in opposite direction

Let us consider the case that the sender $S$ and the receiver $R$ are moving
in opposite direction. To this end we will “reduce” the speeds of the sender
and the receiver and our analysis will show the choice that will yield the
shortest evacuation time. Let us assume that $S$ moves with speed $0\leq
u_{1}\leq 1$ and $R$ takes the bike and moves with speed $0\leq u_{2}\leq v$.
As already proposed, the objective is to find the optimal speeds $u_{1}$ and
$u_{2}$ that should be used by $S$ and $R$, respectively, in order for the
resulting evacuation time to be minimal. There are two cases to take into
account depending on which of the two autonomous robots finds the exit first:

##### Case 1:

$S$ finds the exit first.

As soon as S reaches the exit, it communicates with R which is on the other
side to proceed to the exit. Based on Fig. 1,

Figure 1: Depicted is Case 1 where two robots S and R are moving in opposite
direction such that S finds the exit first.

S needs time $\frac{d}{u_{1}}$ to reach the exit. When S reaches the exit, R
would be at distance $\frac{du_{2}}{u_{1}}$ on the other side and it needs
time $\frac{du_{2}}{vu_{1}}$ to go back to the origin and additional time
$\frac{d}{v}$ to go from the origin to the exit. Thus taking this into account
the total evacuation time will be as follows:

$\displaystyle{\cal E}_{1}=\frac{d}{u_{1}}+\frac{du_{2}}{vu_{1}}+\frac{d}{v}$
(1)

##### Case 2:

$R$ finds the exit first.

Based on Fig. 2,

Figure 2: Depicted is Case 2 where two robots S and R are moving in opposite
direction and such that R finds the exit.

when $R$ reaches the exit, then $S$ is at distance $\frac{du_{1}}{u_{2}}$ to
the left of the origin. In this case $R$ needs to go back to the left to bring
$S$ toward the origin. When $R$ reaches the origin, $S$ will be at distance
$\frac{du_{1}}{u_{2}}+\frac{du_{1}}{v}$ to the left of the origin. At that
point if $S$ meets $R$ at distance $x$ from its current position, then we have
the following:

$\displaystyle\frac{x}{u_{1}}=\frac{du_{1}}{u_{2}v}+\frac{du_{1}}{v^{2}}+\frac{x}{v}$
$\displaystyle\implies
x\left(\frac{v-u_{1}}{u_{1}v}\right)=\frac{du_{1}v+du_{1}u_{2}}{u_{2}v^{2}}$
$\displaystyle\implies x=\frac{du_{1}^{2}v+du_{1}^{2}u_{2}}{(v-u_{1})u_{2}v}$

After $S$ and $R$ meet, they will be away from the origin by the following
distance:

$\displaystyle\frac{du_{1}}{u_{2}}+\frac{du_{1}}{v}+x$
$\displaystyle=\frac{du_{1}}{u_{2}}+\frac{du_{1}}{v}+\frac{du_{1}^{2}v+du_{1}^{2}u_{2}}{(v-u_{1})u_{2}v}=\frac{du_{1}v^{2}+dvu_{1}u_{2}}{vu_{2}(v-u_{1})}$

Thus $S$ and $R$ will be away from the exit by

$d+\frac{du_{1}v^{2}+dvu_{1}u_{2}}{vu_{2}(v-u_{1})}=\frac{du_{1}v+du_{2}v}{u_{2}(v-u_{1})}.$

Using Lemma 1, the time needed to share the bike at distance
$\frac{du_{1}v+du_{2}v}{u_{2}(v-u_{1})}$ from the exit will be as follows:

$\displaystyle\frac{(du_{1}v+du_{2}v)(v+1)}{2vu_{2}(v-u_{1})}=\frac{du_{1}v^{2}+du_{1}v+du_{2}v^{2}+du_{2}v}{2u_{2}v(v-u_{1})}$

Now we can find the evacuation time ${\cal E}_{2}$ as follows:

$\displaystyle{\cal E}_{2}$
$\displaystyle=\frac{d}{u_{2}}+\frac{d}{v}+\frac{du_{1}}{u_{2}v}+\frac{du_{1}}{v^{2}}+\frac{du_{1}^{2}v+du_{1}^{2}u_{2}}{(v-u_{1})u_{2}v^{2}}+\frac{du_{1}v^{2}+du_{1}v+du_{2}v^{2}+du_{2}v}{2u_{2}v(v-u_{1})}$
$\displaystyle=\frac{2dv^{2}+2du_{2}v+du_{1}v^{2}+du_{1}v+du_{2}v^{2}+dvu_{2}}{2u_{2}v(v-u_{1})}.$
(2)

Since we do not know where the exit is situated, we should consider the worst
case scenario of the evacuation time, denoted by ${\cal E}$, to be the maximum
of ${\cal E}_{1}$ and ${\cal E}_{2}$, as given in Equations (1) and (2), as
follows:

$\displaystyle{\cal E}$ $\displaystyle=\max\\{{\cal E}_{1},{\cal E}_{2}\\}$
(3)

In order to minimize the total evacuation time, we need to know for what value
of $v$ is each of ${\cal E}_{1}$ and ${\cal E}_{2}$ minimized.

Using Equation (1), we see that ${\cal E}_{1}$ is minimized if $u_{1}=1$ and
$u_{2}=0$ since we assume that the exit is on the side of the sender $S$. In
this case ${\cal E}_{1}$ will be as follows: ${\cal
E}_{1}=d\cdot\frac{1+v}{v}$.

Using Equation (2), we see that ${\cal E}_{2}$ is minimized if $u_{2}=v$ and
$u_{1}=0$ since we assume that the exit is on the side of the receiver $R$. In
this case ${\cal E}_{2}$ will be as follows: ${\cal
E}_{2}=d\cdot\frac{v^{3}+5v^{2}}{2v^{3}}=d\cdot\frac{v+5}{2v}$. Thus

${\cal E}_{1}-{\cal
E}_{2}=d\cdot(\frac{1+v}{v}-\frac{v+5}{2v})=d\cdot\frac{v-3}{2v}.$

We conclude that ${\cal E}_{1}\leq{\cal E}_{2}$, if $v\leq 3$. Thus we have
two cases to consider based on the speed of the biker.

##### Case 1:

$v\leq 3$
We have ${\cal E}_{1}<{\cal E}_{2}$ and thus the maximum evacuation time is
${\cal E}=\max\\{{\cal E}_{1},{\cal E}_{2}\\}={\cal E}_{2}$ and in this case
${\cal E}_{2}$ is minimized if we set $u_{2}=v$.

##### Case 2:

$3\leq v\leq 10$
(Note that if $v>10$, then the zig-zag algorithm discussed later in section
2.2 performs better.) We have ${\cal E}_{2}<{\cal E}_{1}$ and thus the maximum
evacuation time is ${\cal E}=\max\\{{\cal E}_{1},{\cal E}_{2}\\}={\cal E}_{1}$
and in this case ${\cal E}_{1}$ is minimized when we set $u_{1}=1$.

####  2.1.1 Algorithm 1: $1\leq v\leq 3$

#####

In the first algorithm, the receiver takes the bike and moves in one direction
with its maximum speed, while the sender moves with optimal speed $0\leq
u_{1}\leq 1$ in the other direction. If the sender finds the exit first, then
it informs the receiver and stays at the exit, while if the receiver finds the
exit first, then it changes direction to catch up with the sender and then
both robots ride-share the bike toward the exit. The first algorithm performs
better when the speed of the bike satisfies $1\leq v\leq 3$. If we assume the
sender moves to the left with speed $0\leq u_{1}\leq 1$ and the receiver to
the right with maximum speed $v$.

The algorithm will be as follows:

Algorithm 1 OppDirectionWithBikerMovingAtMaxSpeed ($S$ source, $D$
destination)

1:Sender moves to the left with speed $u_{1}$

2:Receiver takes the bike and moves to the right with speed $v$

3:if Receiver finds the Exit first then

4: it changes direction and catches up to the Sender;

5: The two robots change direction and ride-share the bike towards the Exit;

6:else

7: if Sender finds the Exit first then

8: it informs the Receiver that the Exit has been found but stays at the Exit;

9: The Receiver rides the bike to the Exit;

We prove the following result.

###### Theorem 2.1

The optimal competitive ratio of Algorithm 1 is

$\frac{2v}{v+1}\left(\frac{2v+\frac{-(7v+v^{2})+v\sqrt{v^{2}+30v+97}}{2v+6}}{v\frac{-(7v+v^{2})+v\sqrt{v^{2}+30v+97}}{2v+6}}\right)$

###### Proof

We have that $S$ moves in one direction with speed $0\leq u_{1}\leq 1$ and $R$
moves in the other direction with speed $v$. Thus we will have the following
equations:

$\displaystyle{\cal E}_{1}=\frac{2d}{u_{1}}+\frac{d}{v}$ $\displaystyle{\cal
E}_{2}=\frac{dv^{2}+5dv+du_{1}v+du_{1}}{2v(v-u_{1})}$

In order to determine the optimal speed $u_{1}$ for the sender, we should set
the evacuation times in both cases above to be equal, namely ${\cal
E}_{2}={\cal E}_{1}$. This will give us the following:

$\displaystyle\frac{d(v^{2}+5v+u_{1}v+u_{1})}{2v(v-u_{1})}=\frac{2dv+du_{1}}{u_{1}v}\implies(3+v)u_{1}^{2}+(7v+v^{2})u_{1}-4v^{2}=0$

This is a quadratic equation of $u_{1}$ and it has two solutions. Keeping the
positive root we get

$u_{1}=\frac{-(7v+v^{2})+v\sqrt{v^{2}+30v+97}}{2v+6}.$

We expect the optimal speed $u_{1}$ to be less than $1$ when the speed $v$ is
between $1$ and $3$. If we substitute $u_{1}$ in the evacuation time, denoted
by ${\cal E}$, calculated in any of the two cases we get the following

$\displaystyle{\cal
E}=d\left(\frac{2v+\frac{-(7v+v^{2})+v\sqrt{v^{2}+30v+97}}{2v+6}}{v\frac{-(7v+v^{2})+v\sqrt{v^{2}+30v+97}}{2v+6}}\right)$

Thus, using Lemma 1, the competitive ratio will be

$\frac{2v}{v+1}\left(\frac{2v+\frac{-(7v+v^{2})+v\sqrt{v^{2}+30v+97}}{2v+6}}{v\frac{-(7v+v^{2})+v\sqrt{v^{2}+30v+97}}{2v+6}}\right).$

This completes the proof of Theorem 2.1. ∎

####  2.1.2 Algorithm 2: $3\leq v\leq 10$

#####

In the second algorithm, the sender moves in one direction with unit speed
while the receiver takes the bike and moves in the other direction with
optimal speed $1<u_{2}<v$. If the sender finds the exit first, then it informs
the receiver and stays at the exit, while if the receiver finds the exit
first, then it changes direction to catch up with the sender and then both
robots ride-share the bike toward the exit. The second algorithm performs
better when the speed $3\leq v\leq 10$ as we will see later in figure 4. We
assume the sender moves to the left with unit speed and the receiver to the
right with speed $1\leq u_{2}\leq v$.

The main algorithm will be as follows:

Algorithm 2 OppDirectionWithWalkerMovingAtMaxSpeed ($S$ source, $D$
destination)

1:Sender moves to the left with unit speed

2:Receiver takes the bike and moves to the right with speed $u_{2}$

3:if Receiver finds the Exit first then

4: it changes direction and catches up to the Sender;

5: The two robots change direction and ride-share the bike towards the Exit;

6:else

7: if Sender finds the Exit first then

8: it informs the Receiver that the Exit has been found but stays at the Exit;

9: The Receiver rides the bike to the Exit;

We prove the following result.

###### Theorem 2.2

The optimal competitive ratio of Algorithm 2 is:

$\displaystyle\frac{2v}{v+1}\left(1+\frac{1}{v}+\frac{v^{2}-3v-2-\sqrt{v^{4}+18v^{3}-7v^{2}+4v+4}}{4v(1-v)}\right)$

###### Proof

$S$ moves in one direction with unit speed and $R$ moves in the other
direction with speed $3\leq u_{2}\leq v$. We will have the following:

$\displaystyle{\cal E}_{1}=\frac{dv+du_{2}+d}{v}\mbox{ and }{\cal
E}_{2}=\frac{3dv+du_{2}v+3du_{2}+d}{2u_{2}(v-1)}$

The optimal speed $u_{2}$ is obtained by setting the evacuation time in the
two cases above to be equal, namely ${\cal E}_{2}={\cal E}_{1}$. Thus we have
the following:

$\displaystyle\frac{u_{2}+v+1}{v}=\frac{3v^{2}+3u_{2}v+u_{2}v^{2}+v}{2u_{2}v(v-1)}$

If we multiply out we will end up having the following quadratic equation in
the variable $u_{2}$: $(2-2v)u_{2}^{2}+(3v-v^{2}+2)u_{2}+3v^{2}+v=0.$ Keeping
the positive root gives:

$\displaystyle
u_{2}=\frac{v^{2}-3v-2-\sqrt{v^{4}+18v^{3}-7v^{2}+4v+4}}{4(1-v)}$

Thus the evacuation time will be

$\displaystyle
d\left(1+\frac{1}{v}+\frac{v^{2}-3v-2-\sqrt{v^{4}+18v^{3}-7v^{2}+4v+4}}{4v(1-v)}\right)$

Using Lemma 1, it follows that the competitive ratio will be

$\displaystyle\frac{2v}{v+1}\left(1+\frac{1}{v}+\frac{v^{2}-3v-2-\sqrt{v^{4}+18v^{3}-7v^{2}+4v+4}}{4v(1-v)}\right)$

This completes the proof of Theorem 2.2. ∎

###  2.2 Algorithm 3: $10<v$

In the next Algorithm 3, the sender takes the bike and uses a doubling
strategy to search for the exit and moves a distance $2^{k}$ during the $k$-th
iteration. The receiver also uses a doubling strategy but since it is moving
with unit speed it will try to stay as close as possible to the biker. This
can be achieved by having the receiver moves a distance $\frac{2^{k}}{v}$
during the $k^{th}$ iteration, since moving any further will cause the sender
to be farther away from the receiver during the $(k+1)^{st}$ iteration.

In the third algorithm the Receiver imitates the Sender and as it will be seen
in the sequel it performs better when the speed $v$ satisfies $v>10$. The main
algorithm is as follows:

Algorithm 3 (ReceiverImitateSender)

1:for $k\leftarrow 1$ to $\infty$ do

2: if k is odd (resp.even) then

3: Sender takes the bike and moves right (resp. left) a distance $2^{k}$
unless the exit is found;

4: Receiver moves right (resp. left) a distance $\frac{2^{k}}{v}$;

5: if exit is found by sender then

6: Communicate with receiver;

7: Sender moves back $\frac{d}{2}-\frac{d}{2v}$ to leave the bike for the
receiver and then returns to exit;

8: Receiver continues toward the exit after picking up the bike;

9: Quit;

10: Sender turns; then moves left (resp. right), returns to the origin;

11: Receiver turns; then moves left (resp. right), returns to the origin;

###### Theorem 2.3

The competitive ratio for Algorithm 3 is at most

$\left(\frac{2v}{v+1}\right)\left(\frac{9}{v}+\frac{1}{2}-\frac{1}{2v^{2}}\right).$

###### Proof

In this algorithm the sender uses a doubling strategy with maximum speed $v$.
The receiver will follow the sender but will move $\frac{2^{k}}{v}$ in each
iteration instead of $2^{k}$. The sender will reach the exit first then it
will communicate with the receiver to proceed to the exit. The sender will go
back distance $\frac{d}{2}-\frac{d}{2v}$ to drop off the bike so that the
receiver can pick it up on its way to the exit. We will justify why the sender
needs to move $\frac{d}{2}-\frac{d}{2v}$ after reaching the exit to leave the
bike for the receiver. After the sender reaches the exit, there is no benefit
to stay at the exit with the bike since the receiver which is moving with unit
speed can benefit from the bike to reach the exit faster.

The key to finding the distance $x$ which is the distance between the exit and
the point where the bike is dropped off is to have the sender drop it off at a
point such that when it goes back to the exit it will reach the exit at the
same time as the receiver.

Figure 3: Graph depicts the minimal distance $x$ away from the origin at which
the bike is dropped off by the sender and allows the receiver to take the bike
and reach the exit sooner.

Based on Fig. 3, if we recall that $d$ is the distance from the origin to the
exit and $x$ is the distance from the exit to the point where the sender drops
off the bike, then we have the following:

$\displaystyle d-x+\frac{x}{v}=\frac{d}{v}+\frac{x}{v}+x$
$\displaystyle\implies x=\frac{d}{2}-\frac{d}{2v}$

This will guarantee that when the sender drops off the bike at distance $x$,
it will reach the exit at the same time as the receiver.

Assume that the exit is found during the $k^{th}$ iteration, then
$2^{k-2}<d\leq 2^{k}$. We can calculate the evacuation time, denoted by ${\cal
E}$, as follows:

$\displaystyle{\cal E}$ $\displaystyle=\frac{2\cdot 2^{0}}{v}+\frac{2\cdot
2^{1}}{v}+\dots+\frac{2\cdot 2^{k-1}}{v}+d-x+\frac{x}{v}$
$\displaystyle=\frac{2(2^{k}-1)}{v}+d-x+\frac{x}{v}$
$\displaystyle=\frac{2^{k+1}}{v}-\frac{2}{v}+d-\frac{d}{2}+\frac{d}{2v}+\frac{d}{2v}-\frac{d}{2v^{2}}$
$\displaystyle=\frac{2^{k+1}}{v}-\frac{2}{v}+\frac{d}{2}+\frac{d}{v}-\frac{d}{2v^{2}}$

Combining similar terms above we obtain

$\displaystyle{\cal E}$ $\displaystyle\leq
2^{3}\cdot\frac{2^{k-2}}{v}-\frac{2}{v}+\frac{d}{2}+\frac{d}{v}-\frac{d}{2v^{2}}$
$\displaystyle\leq\frac{8d}{v}-\frac{2}{v}+\frac{d}{2}+\frac{d}{v}-\frac{d}{2v^{2}}$
$\displaystyle\leq\frac{9d}{v}+\frac{d}{2}-\frac{d}{2v^{2}}-\frac{2}{v}$
$\displaystyle\leq\frac{9d}{v}+\frac{d}{2}-\frac{d}{2v^{2}}$

Thus, using Lemma 1, the resulting competitive ratio will be upper bounded by
the quantity
$\left(\frac{2v}{v+1}\right)\left(\frac{9}{v}+\frac{1}{2}-\frac{1}{2v^{2}}\right).$
This completes the proof of Theorem 2.3. ∎

###  2.3 Performance of the Upper bound algorithms

The performance of the three algorithms is compared in Fig. 4. As it can be
seen in this figure, the performance of the algorithm depends on the speed $v$
of the non-autonomous robot.

Figure 4: Graph showing the performance of the 3 algorithms by showing how the
competitive ratio fluctuates in terms of the speed. Note in the pictures above
Algorithms 1 is depicted in the range $[1,3]$ and Algorithms 2 and 3 are
depicted in the range $[1,+\infty)$.

From the plots we can confirm that Algorithm 1 performs best if the speed is
less than 3, while Algorithm 2 performs best if the speed is between $3$ and
$10$. Finally, Algorithm 3 performs best if the speed is larger than 10 and in
this case the competitive ratio converges to $1$ as $v\to\infty$, which is
also the best that could be achieved if both autonomous robots have full
knowledge in that they know where the exit is.

##  3 Lower Bound

Finally, we prove a lower bound.

###### Theorem 3.1

The lower bound for the competitive ratio is

$\left\\{\begin{array}[]{cl}\frac{6}{v+1}&\mbox{ if $v\leq 3$}\\\
\frac{v^{2}+2v-3}{v^{2}-1}&\mbox{ if $v>3$}\end{array}\right.$

(Proof of Theorem 3.1 can be found in the appendix.)

##  4 Conclusion

The purpose of our investigations was to study the limits of communication for
robots with mixed faults. In particular, in this paper we gave evacuation
algorithms and studied the upper and lower bounds for evacuation in the S/R
communication model where the faults involve a sender and a receiver seeking
to evacuate through an unknown exit situated on an infinite line.

In particular, in Section 3 we presented a lower bound and in subsection 2.3
we presented the asymptotic behaviour of the main algorithm. However the
possibility of improving the lower bound remains an open problem.

An extension of the problem already analyzed could also be considered for the
case of multiple (more than two) autonomous robots in a setting where the
autonomous robots may suffer any of S/R, crash, and byzantine faults. One
could also consider the case of multiple non-autonomous robots. The
possibility of studying other search domains (geometric star graph, disk,
polygon, etc) would also be interesting.

## References

  * [1] R. Ahlswede and I. Wegener.  Search problems.  Wiley-Interscience, 1987.
  * [2] R. Baeza Yates, J. Culberson, and G. Rawlins.  Searching in the plane.  Information and Computation, 106(2):234–252, 1993.
  * [3] R. Baeza-Yates and R. Schott.  Parallel searching in the plane.  Computational Geometry, 5(3):143–154, 1995.
  * [4] E. Bampas, J. Czyzowicz, L. Gasieniec, D. Ilcinkas, R. Klasing, T. Kociumaka, and D. Pajak.  Linear search by a pair of distinct-speed robots.  Algorithmica, 81(1):317–342, 2019.
  * [5] A. Beck.  On the linear search problem.  Israel Journal of Mathematics, 2(4):221–228, 1964.
  * [6] R. Bellman.  An optimal search.  SIAM Review, 5(3):274–274, 1963.
  * [7] P. Bose and J.-L. De Carufel.  A general framework for searching on a line.  Theoretical Computer Science, pages 703:1–17, 2017.
  * [8] S. Brandt, F. Laufenberg, Y. Lv, D. Stolz, and R. Wattenhofer.  Collaboration without communication: Evacuating two robots from a disk.  In International Conference on Algorithms and Complexity, pages 104–115. Springer, 2017.
  * [9] M. Chrobak, L. Gasieniec, G. T., and R. Martin.  Group search on the line.  In SOFSEM, pages 164–176. Springer, 2015.
  * [10] J. Czyzowicz, S. Dobrev, K. Georgiou, E. Kranakis, and F. MacQuarrie.  Evacuating two robots from multiple unknown exits in a circle.  Theoretical Computer Science, 709:20–30, 2018.
  * [11] J. Czyzowicz, K. Georgiou, M. Godon, E. Kranakis, D. Krizanc, W. Rytter, and M. Włodarczyk.  Evacuation from a disc in the presence of a faulty robot.  In International Colloquium on Structural Information and Communication Complexity, pages 158–173. Springer, 2017.
  * [12] J. Czyzowicz, K. Georgiou, and E. Kranakis.  Group search and evacuation.  In Distributed Computing by Mobile Entities, pages 335–370. Springer, 2019.
  * [13] J. Czyzowicz, K. Georgiou, E. Kranakis, D. Krizanc, L. Narayanan, J. Opatrny, and S. Shende.  Search on a line by byzantine robots.  In ISAAC, pages 27:1–27:12, 2016.
  * [14] J. Czyzowicz, R. Killick, E. Kranakis, D. Krizanc, L. Narayanan, J. Opatrny, D. Pankratov, and S. Shende.  Group evacuation on a line by agents with different communication abilities.  ISAAC 2021, pages 57:1–57:24, 2021.
  * [15] J. Czyzowicz, R. Killick, E. Kranakis, and G. Stachowiak.  Search and evacuation with a near majority of faulty agents.  In SIAM Conference on Applied and Computational Discrete Algorithms (ACDA21), pages 217–227. SIAM, 2021.
  * [16] J. Czyzowicz, E. Kranakis, D. Krizanc, L. Narayanan, and J. Opatrny.  Search on a line with faulty robots.  Distributed Computing, 32(6):493–504, 2019.
  * [17] E. D. Demaine, S. P. Fekete, and S. Gal.  Online searching with turn cost.  Theoretical Computer Science, 361(2):342–355, 2006.
  * [18] K. Georgiou, N. Giachoudis, and E. Kranakis.  Evacuation from a disk for robots with asymmetric communication.  In 33rd International Symposium on Algorithms and Computation (ISAAC 2022). Schloss Dagstuhl-Leibniz-Zentrum für Informatik, 2022.
  * [19] K. Georgiou, E. Kranakis, N. Leonardos, A. Pagourtzis, and I. Papaioannou.  Optimal circle search despite the presence of faulty robots.  In International Symposium on Algorithms and Experiments for Sensor Systems, Wireless Networks and Distributed Robotics, pages 192–205. Springer, 2019.
  * [20] K. Jawhar and E. Kranakis.  Bike assisted evacuation on a line.  In SOFSEM (47th International Conference on Current Trends in Theory and Practice of Computer Science). Springer, LNCS, 2021.
  * [21] L. Stone.  Theory of optimal search.  Academic Press New York, 1975.

##  Appendix 0.A Proof of Lemma 1

###### Proof

(Lemma 1) Based on Fig. 5, both robots are initially situated at the origin
which is at distance $d$ from the exit.

Figure 5: Graph depicts the movement of two robots sharing the bike and
starting at the origin in order to cover distance $d$.

The objective of both robots is to ride share the bike and arrive at the exit
at the same time. If we assume the first robot drops off the bike at distance
$x$ then we should have the following:

$\displaystyle\frac{x}{v}+d-x=x+\frac{d-x}{v}.$

Solving the last equation in $x$ we derive that $x=\frac{d}{2}$. Therefore the
time it takes for the two robots to find the Exit is
$\frac{x}{v}+d-x=\frac{d(v+1)}{2v}$. For a proof of the lower bound see [20].
∎

##  Appendix 0.B Proof of Theorem 3.1

###### Proof

(Theorem 3.1) The Exit is placed at one of the endpoints $\pm d$ while the
robots and the bike start at the origin $O$.
In order to find the lower bound we consider the worst case scenario which is
achieved when the biker takes the bike and reaches $-d$ and the exit is at
$d$. Depending on the speed of the biker, there are two cases to consider as
follows:

##### Case 1:

$v\geq 3$

As we can see in Fig. 6, if the biker is fast then it can reach the exit
before the other robot and this occurs if the speed $v\geq 3$.

Figure 6:  To achieve the optimal lower bound, both robots move in opposite
direction. When the biker reaches $-d$ in time $\frac{d}{v}$, the adversary
places the exit at $d$. At time $\frac{d}{v}$, the hiker would be at point
$\frac{d}{v}$.

In this case the biker needs to share the bike with the other robot. The
closest point $P$ to the exit that the other robot can reach is $\frac{d}{v}$
to the right of the origin. In order to find the evacuation time, we assume
that the biker meets the other robot at distance $x$ to the right of point
$P$. In order to achieve minimal evacuation time, we have the following:

$\displaystyle\frac{d}{v}+\frac{d}{v^{2}}+\frac{x}{v}+d-x-\frac{d}{v}=x+\frac{1}{v}(d-\frac{d}{v}-x)$
$\displaystyle\implies\frac{d}{v}+\frac{d}{v^{2}}+\frac{x}{v}+d-\frac{d}{v}-x=x+\frac{d}{v}-\frac{d}{v^{2}}-\frac{x}{v}$
$\displaystyle\implies 2x-\frac{2x}{v}=\frac{2d}{v^{2}}+d-\frac{d}{v}$
$\displaystyle\implies x=\frac{2d+dv^{2}-dv}{2v(v-1)}$

As a consequence, if ${\cal E}$ denotes the evacuation time then we have

$\displaystyle{\cal E}$
$\displaystyle=\frac{d}{v}+x+\frac{d}{v}-\frac{x}{v}-\frac{d}{v^{2}}$
$\displaystyle=\frac{d}{v}+\frac{d(v^{2}-v+2)}{2v(v-1)}+\frac{d}{v}-\frac{d}{v^{2}}-\frac{dv^{2}-dv+2d}{2v^{3}-2v^{2}}$
$\displaystyle=\frac{d(v^{2}+2v-3)}{2v(v-1)}$

##### Case 2:

$v\leq 3$

In this case the biker cannot arrive to the exit before the other robot and
therefore the minimal evacuation time will be $\frac{3d}{v}$

Using Lemma 1, we conclude that a lower bound on the competitive ratio,
denoted by $CR$, will satisfy the following inequalities:

##### Case 1

: $v\geq 3$

$\displaystyle\frac{v^{2}+2v-3}{v^{2}-1}\leq\frac{d(v^{2}+2v-3)}{2v(v-1)}\frac{2v}{d(v+1)}\leq
CR$ (4)

##### Case 2

: $v\leq 3$

$\displaystyle\frac{6}{v+1}\leq\frac{3d}{v}\frac{2v}{d(v+1)}\leq CR$ (5)

The lower bound follows from Inequalities (4) and (5). ∎