# Exploiting Localized Transition Waves to Tune Sound Propagation in Soft
Materials

Audrey A. Watkins†, Austin Eichelberg† Osama R. Bilal<EMAIL_ADDRESS>
$^†$ Both authors contributed equally to this work. Department of Mechanical
engineering, University of Connecticut, Storrs, CT 06269, USA

###### Abstract

Programmable materials hold great potential for many applications such as
deployable structures, soft robotics, and wave control, however, the presence
of instability and disorder might hinder their utilization. Through a
combination of analytical, numerical, and experimental analyses, we harness
the interplay between instabilities, geometric frustration, and mechanical
deformations to control the propagation of sound waves within self-assembled
soft materials. We consider levitated magnetic disks confined by a magnetic
boundary in-plane. The assemblies can be either ordered or disordered
depending on the intrinsic disk symmetry. By applying an external load to the
assembly, we observe the nucleation and propagation of different topological
defects within the lattices. In the presence of instabilities, the defect
propagation gives rise to time-independent localized transition waves.
Surprisingly, in the presence of frustration, the applied load briefly
introduces deformation-induced order to the material. By further deforming the
lattices, new patterns emerge across all disk symmetries. We utilize these
patterns to tune sound propagation through the material. Our findings could
open new possibilities for designing exotic materials with potential
applications ranging from sound control to soft robotics.

††preprint: APS/123-QED

Patterns emerge within an assembly when the system elements are left to
interact based on their individual energy. When an external load is applied to
these self-assembled patterns, the elements keep reorienting themselves to
minimize their energy, resulting in different patternsTaheri  _et al._ (2015).
In the case of incompatible elements, geometric frustrationMellado  _et al._
(2012); Schönke  _et al._ (2015); Wang  _et al._ (2017); Kang  _et al._ (2014)
can arise, particularly, when the global energy minima of the system is
incompatible with the minima of the individual elementsDiep  _et al._ (2013).
Such geometric frustration can prohibit the self-assembly from achieving long-
range order. Furthermore, the presence of instabilities can influence the
emerging order (or disorder) within a tunable materialKochmann and Bertoldi
(2017); Li  _et al._ (2019); Goshkoderia  _et al._ (2020); Jin  _et al._
(2020); Haghpanah  _et al._ (2016); Florijn  _et al._ (2014, 2016). A valid
path for functional tunability of matter is external stimuli (such as heat or
magnetic fields) which can effectively apply an external load to tune material
properties Matar  _et al._ (2013); Bilal  _et al._ (2017a, b); Palermo  _et
al._ (2019); Wang  _et al._ (2018, 2020). For soft materials, however, this
usually translates to large deformation with the potential rise of
instabilities, disorder, incompatibility and geometric frustration. While such
materials hold great potential for many applications such as deployable
structures, soft roboticsChen  _et al._ (2018); Rothemund  _et al._ (2018);
Zareei  _et al._ (2020) and wave controlBertoldi and Boyce (2008); Wang  _et
al._ (2013); Shan  _et al._ (2014); Shim  _et al._ (2015); Raney  _et al._
(2016); Bilal  _et al._ (2017a, b); Li  _et al._ (2012), the presence of
frustration, instability and disorder within tunable soft materials might
hinder their utilization. Here, we analytically, numerically, and
experimentally show that in the presence of instabilities and disorder,
pattern evolution due to deformation can be harnessed. By deforming our
lattices, we control the nucleation and propagation of topological defects as
time-independent localized transition waves. Furthermore, we exploit the
evolving patterns to tune the propagation of sound waves within our material
systems.

Figure 1:  Concept. Minimum energy state for 1 disk (top) and 3 disks (middle)
confined by a fixed boundary with (a) two, (b) three and (c) four magnets per
disk. The potential energy as a function of the disk phase-angle for the
central disk is shown as a polar plot (bottom).

We consider the mechanics of levitated disks (lying in the horizontal plane)
with repulsive magnetic interactions within a fixed magnetic boundary. The
boundary magnets are equally spaced. The embedded disk magnets vary between
1-4, resembling the shape of a dot, line, triangle and a square (Fig.1). The
boundary magnets create energy potential-wells dictating the alignment of the
free-floating disks. When confined within an energy potential well, a single
disk translates and rotates freely to minimize its own energy. For a disk with
one magnet, the minimum energy state is at the center of the energy potential
well (i.e., the center of the four boundary magnets). In the case of two or
four magnets per disk, the disk has four degenerate minimum energy states
located at the center of the potential well with a 90∘ phase angle $\phi$
(Fig.1a,c). In the case of three magnets per disk, there exist multiple
minimum energy states that are off-center. Each of these off-center positions
has three minimum energy states with a 120∘ phase angle (Fig.1b). When
multiple disks are confined within the boundary, the self-alignment must take
into account the added nearest neighbour interaction. For both two and four
magnets per disk, the minimum energy states exist at the center of the
potential wells, resulting in perfectly aligned self-oriented lattices
(Fig.1a,c middle). In the case of three magnets per disk, the disks have no
clear state of self-alignment that minimizes both their individual and total
energies, hence the system is frustrated (Fig.1c). To experimentally validate
the potential energy calculations (Fig.1a-c), we enclose each disk type (i.e.,
with 1-4 magnets) within a fixed boundary atop an air bearing table. Both
examples of the numerical and the experimental self-aligned disks with 1-4
embedded magnets are included in the Supplementary Material. The disks start
in random positions relative to the boundary and to each other. Then, we
pressurize the air bearing, activating a layer of laminar air flow beneath the
surface of the disks to allow them to float freely in the horizontal plane
(i.e., levitating the disks similar to the arcade game “Air Hockey”).

To study the wave propagation characteristics of the assemblies, we first
consider a model of infinite periodic lattices with a single disk and its four
boundary magnets as the unit cell with nearest neighbour interaction (i.e.,
the immediately adjacent disks). Only the free floating disks are allowed to
move on both sides of the considered disk. Each disk has two degrees of
freedom in the $x$ and $y$ directions. The dispersion relation of the system
can be calculated using Bloch’s theoremBloch (1929) as:
$[-\omega^{2}\textbf{M}+\textbf{K}(\bm{\kappa})]\bm{\phi}=0,$ where $\omega$
is the frequency, $\kappa$ is the wavenumber, $\bm{\phi}=[u~{}v]^{T}$ is the
Bloch displacement vector in the $x$ and $y$ directions, $\bm{M}$ is the mass
matrix, and $\bm{K(\kappa)}$ is the stiffness matrixWatkins and Bilal (2020),
taking into consideration the static repulsive forces between the disks in the
equilibrium configurationJiao and Gonella (2020) (See Supp. Mat.). For disks
with three embedded magnets, we assume perfect alignment of the disks in order
to use Bloch’s theorem. The analytically computed dispersion curves show two
distinctive bands representing the longitudinal and shear motion in each of
the four lattices (Fig.2a-d). In addition, as the number of magnets per disk
increases, the frequency of the bands increases, due to the increase in
stiffness within the lattice.

Figure 2:  Lattice wave propagation. (a-d) Analytical (black) and experimental
(contour) dispersion curves for 1-4 magnets per disk. Numerical (orange) and
experimental (blue) transmission at a central disk within the assembly. (i-l)
The initial and final positions of ten disks with 1-4 embedded magnets. All
the lattices are perfectly aligned except for 3 magnets per disk in panel (k).

To numerically verify the infinite model, we consider a finite arrangement of
10 disks confined to a magnetic boundary (using Verlet methodPress  _et al._
(2007)). The disks can move freely in-plane. The simulations are initialized
with random disk positions. After the model reaches equilibrium, based on the
balance of repulsive forces, we excite the right most disk with a chirp signal
between 0.2 and 10 Hz. The fast Fourier transform (FFT) of the longitudinal
transmission ranges match closely with the analytical model (Fig.2e-h).
Remarkably, in the case of frustrated disks (i.e., with 3 magnets), the
simulated transmission through the disordered lattice closely matches the
dispersion curve calculated for a perfectly ordered lattice.

To experimentally verify the dispersion curves, we harmonically excite the
longitudinal mode of the floating self-aligned lattices (Fig.2i-l) using a
mechanical shaker with a chirp signal between 0.2 - 10 Hz and measure the
x-displacement of the disks through digital image correlation (DICe). The
longitudinal excitation takes place at the right most disk. The oscillatory
motion of the disks is processed using 1D and 2D fast Fourier transform,
resulting in a transmission spectrum and an experimental dispersion curve,
respectively (See Supp. Mat.).

The experimentally measured transmission ranges (1D-FFT) and dispersion curves
(2D-FFT) are in agreement with our numerical predictions within the
longitudinal mode (Fig.2e-h). Once more, in presence of frustration, the
experimentally measured transmission and dispersion through the disordered
lattice matches closely with a perfectly ordered lattice. In addition, we note
the emergence of an extra transmission band around 1 Hz for both measured
transmission and dispersion. This transmission anomaly dissipates as the wave
propagates further down the lattice due to the increased disorder. We note
that there is no experimental transmission correlating with the shear modes,
due to the absence of shear excitation.

A viable strategy for tuning wave propagation in soft material is applying an
external load to tune stiffnessBoechler  _et al._ (2011); Wang  _et al._
(2014); Florijn  _et al._ (2014); Babaee  _et al._ (2016); Amendola  _et al._
(2018), however, this usually tunes the material globally. An intriguing
feature of our material is the localized pattern transformation enabled by the
boundary-induced multi-stability of the potential wells. When a load is
applied at one end of the material, the disks initially experience an added
global stiffness. Once the applied load passes a threshold, the first disk
overcomes its boundary potential and snaps into the next potential well. As
the compression at the boundary increases, the migration of disks between
potentials creates a propagating defect, which gives rise to a stable
transition wave. If the compression pauses, the propagation of the defects
throughout the lattice stops (i.e., the transition wave is time independent).

Figure 3:  Time-independent transition wave. Numerically obtained energy per
unit cell as the lattice is compressed for a) one (b) two (c) three and (d)
four embedded magnets per disk.  Figure 4:  (a-d)Defect nucleation and
propagation. Step-wise deformation of ten levitated disks with 1-4 embedded
magnets per disk. The defects are marked with with the same color dashed box
as they propagate through the material. (e-h) The different types of defects
arising for various disk types. The perfectly aligned state is highlighted
with a green solid rectangle (top picture in each panel except for the
frustrated state in c). The normals are highlighted in yellow throughout all
panels.

We numerically simulate the compression of a 50 disk chain for each disk type
using the Verlet method (Fig.3). The disks start at random positions (and
phase angle in cases of more than one magnet per disk) within their potential
wells. After reaching equilibrium, we periodically increase the load applied
at the leftmost disk with intermittent pauses allowing the disks to
reassemble. In the case of one magnet per disk, the material starts at an
equilibrium state with all disks experiencing the same repulsion forces,
except those at the boundary (Fig.3a). As the load increases, the leftmost
disk escapes its potential well, snaps into the next well, and consequently
creates the first defect in the chain. This takes place once the disk
overcomes the energy barrier formed by the vertical pair of boundary magnets.
When the load increases further, the second disk snaps into the next potential
well and the defect propagates further down the chain as another defect
nucleates at the load interface. The transition wave propagating into the
chain creates a checker-like energy pattern (Fig.3a). The same phenomenon
arises when compressing a chain of disks with two and four embedded magnets
(Fig.3b,d). In the case of two-magnets per disk, the energy pattern is
distinctly different from the previous pattern (for 1 magnet per disk) and the
transition wave propagates further down the chain at the same rate of
compression as panel a. However, the transition wave does not reach the end of
the chain in either case. At four-magnets per disk, the transition wave
reaches the end of the chain at ($\approx$ 0.25m) compression, significantly
changing the emerging energy pattern. In the case of 3 magnets per disk, the
equilibrium positions before deformation have a high energy state due to the
frustration-induced disorder. This state is elucidated by the thin green line
at top of Figure 3c, which is consistent with our previous experimental
observation (Fig.2k). Surprisingly, the frustration-induced disorder vanishes
after applying a small load to the chain. This compression ($\approx$ 0.01m)
causes all the disks to align perfectly with a minimum energy state (a thin
horizontal dark blue line in Figure 3c). The self-alignment, or this emergent
deformation-induced order happens every time the chain is compressed by a full
unit cell length (i.e., 0.028m), appearing as flat blue lines in Figure 3c.

To experimentally verify our numerical observations, we consider a chain of 10
free-floating disks within the same magnetic boundary (Fig.2). We apply a load
at the leftmost disk, at equilibrium (Fig.4a-d), causing the disk to snap into
the next potential well. This nucleates the first defect (highlighted with
dashed red boxes in Figure4a-d). In the case of 3 magnets per disk, we
experimentally observe the emergent deformation-induced order. All the ordered
lattices are highlighted green in Figure 4a-d. As the load increases, the
defect migrates further down the chain and a new defect arises (highlighted
with dashed orange boxes in Figure 4a-d). As more defects nucleate, a pattern
of alternating defect and defect-free cells emerges. This is congruent with
the energy patterns observed numerically (Fig.3).

We observe different defect topologies depending on the disk type. For disks
with one embedded magnet, we identify one defect type as: two disks, aligned
horizontally, confined within the same potential well (Fig.4e). In the case of
two magnets per disk, we observe two distinct defects with two disks per
potential well. The normal to the line connecting the centers of the two
defect disks points either left or right (Fig.4f). In the case of three
magnets per disk, we observe a combination of these three defect topologies
(i.e., flat, pointing left or pointing right) (Fig.4g). In the case of four
magnets per disk, the defect is composed of three disks, not two, sharing two
adjacent potential wells. There is either a horizontal alignment of the three
disks or a mix of left and right pointing normals (Fig.4h).

The deformation-induced pattern transformation displays a distinct topological
signature with different disk couplings. This translates to a change in
material properties as patterns evolve. We experimentally test the dynamical
characteristics of the emerging patterns by exciting three of the newly formed
lattices for each of the disk types (Fig.5). For reference, we first
harmonically excite the undeformed self-aligned chain of ten magnets with a
chirp signal between 0.2-35 Hz and record the transmission through the lattice
(1D-FFT). The most left disk is then compressed until a new pattern emerges
across the ten disks. The chain is then excited with the same chirp signal.
The experiment is repeated for two emerging patterns (panels (ii,iii) in
Fig.5) in addition to the pattern of the undeformed disks (panels (i) in
Fig.5). The basic repeating pattern in each case is also included
(Fig.5-insets).

With the one-magnet per disk arrangement, the reference pattern has a single
longitudinal transmission band (Fig.5a(i)). This single band turns into two
distinct bands when the pattern is a defect cell sandwiched between two
defect-free cells (Fig.5a(ii)). As the load increases and the pattern turns
into a single defect-free cell and a defect cell, both bands widen and shift
to higher frequencies (Fig.5a(iii)). The same phenomenon takes place for disks
with two embedded magnets. The reference assembly with perfectly aligned disks
transitions into a defect sandwiched between two defect-free cells. The
pattern then transforms into a single defect-free and a defect cell, with the
emergence of two wider transmission bands at higher frequencies
(Fig.5b(i-iii)). In the three magnets per disk case, the initial pattern is
frustrated with no long range order. As observed earlier (Fig.2g), two
transmission bands are present (Fig.5c(i)). Once the lattice is compressed,
the deformation-induced order emerges and the low frequency disorder-induced
band disappears (Fig.5c(ii)). As the load increases, the pattern transforms
further into a defect cell sandwiched between two defect-free cells
(Fig.5c(iii)). Finally, At four magnets per disk, the lattice transforms from
aligned disks in a 0∘ phase angle into two neighbouring lattices with 0∘ and
45∘ phase angles (Fig.5d(ii)). The difference in transformation for this disk
type is a combination of both the number of disks within the defect (i.e., 3
instead of 2) and the small number of disks considered in the experiment
(i.e., only 10 disks). As the load increases, two transmission bands emerge as
the pattern transforms into an alternation of two defect-free cells and one
defect cell (Fig.5d(iii)).

Figure 5:  Tunable wave propagation. (a-d) Experimentally measured
transmission in a lattice with 1-4 embedded magnets per disk. (i) Transmission
through the uncompressed configuration for each disk type. (ii-iii)
Transmission within the lattice after deformation. The insets show the
corresponding basic building block in each case.

In summary, we consider free-floating magnetic disks within a magnetic
boundary with periodic potentials. The levitated disks perfectly align in the
case of 1, 2 and 4 magnets per disk. In the case of 3 magnets per disk, the
global lattice assembly shows no long range order due to the presence of
geometric frustration. Despite the clear disorder, the transmission within the
lattice closely resembles that of a perfectly ordered lattice. For all
different disk types, we further transform the self-alignment of the disks by
applying a load at one end of each emerging lattice. We observe the nucleation
and propagation of different defect types (depending on the number of magnets
per disk) within the lattices, giving rise to time-independent stable
transition waves due to instabilities. Remarkably, in the presence of
frustration (i.e., 3 magnets per disk), the applied load briefly introduces a
stable order to the assembly. By further deforming the lattices, new patterns
emerge across all disk types. We harness these emerging patterns to tune the
wave propagation characteristics of the soft material. Our findings shed light
on the understanding of deformation-induced pattern transformations,
particularly the rise of frustration-induced disorder and deformation-induced
order. The nature of the localized and stable defect propagation, in addition
to the nonlinear coupling potentials, can be harnessed in designing the next
generation of soft materials for applications ranging from sound manipulation
to soft robotics.

###### Acknowledgements.

This work is supported by the University of Connecticut’s start-up fund (ORB).

## References

  * Taheri  _et al._ (2015) S. M. Taheri, M. Michaelis, T. Friedrich, B. Förster, M. Drechsler, F. M. Römer, P. Bösecke, T. Narayanan, B. Weber, I. Rehberg, _et al._ , Proceedings of the National Academy of Sciences 112, 14484 (2015).
  * Mellado  _et al._ (2012) P. Mellado, A. Concha, and L. Mahadevan, Physical review letters 109, 257203 (2012).
  * Schönke  _et al._ (2015) J. Schönke, T. M. Schneider, and I. Rehberg, Physical Review B 91, 020410 (2015).
  * Wang  _et al._ (2017) P. Wang, Y. Zheng, M. C. Fernandes, Y. Sun, K. Xu, S. Sun, S. H. Kang, V. Tournat, and K. Bertoldi, Physical Review Letters 118, 084302 (2017).
  * Kang  _et al._ (2014) S. H. Kang, S. Shan, A. Košmrlj, W. L. Noorduin, S. Shian, J. C. Weaver, D. R. Clarke, and K. Bertoldi, Physical review letters 112, 098701 (2014).
  * Diep  _et al._ (2013) H. Diep _et al._ ,  _Frustrated spin systems_ (World Scientific, 2013).
  * Kochmann and Bertoldi (2017) D. M. Kochmann and K. Bertoldi, Applied mechanics reviews 69 (2017).
  * Li  _et al._ (2019) J. Li, T. D. Pallicity, V. Slesarenko, A. Goshkoderia, and S. Rudykh, Advanced Materials 31, 1807309 (2019).
  * Goshkoderia  _et al._ (2020) A. Goshkoderia, V. Chen, J. Li, A. Juhl, P. Buskohl, and S. Rudykh, Physical Review Letters 124, 158002 (2020).
  * Jin  _et al._ (2020) L. Jin, R. Khajehtourian, J. Mueller, A. Rafsanjani, V. Tournat, K. Bertoldi, and D. M. Kochmann, Proceedings of the National Academy of Sciences 117, 2319 (2020).
  * Haghpanah  _et al._ (2016) B. Haghpanah, L. Salari-Sharif, P. Pourrajab, J. Hopkins, and L. Valdevit, Advanced Materials 28, 7915 (2016).
  * Florijn  _et al._ (2014) B. Florijn, C. Coulais, and M. van Hecke, Physical review letters 113, 175503 (2014).
  * Florijn  _et al._ (2016) B. Florijn, C. Coulais, and M. van Hecke, Soft matter 12, 8736 (2016).
  * Matar  _et al._ (2013) O. B. Matar, J. Vasseur, and P. A. Deymier, in  _Acoustic Metamaterials and Phononic Crystals_ (Springer, 2013) pp. 253–280.
  * Bilal  _et al._ (2017a) O. R. Bilal, A. Foehr, and C. Daraio, Proceedings of the National Academy of Sciences 114, 4603 (2017a).
  * Bilal  _et al._ (2017b) O. R. Bilal, A. Foehr, and C. Daraio, Advanced materials 29, 1700628 (2017b).
  * Palermo  _et al._ (2019) A. Palermo, Y. Wang, P. Celli, and C. Daraio, Physical Review Applied 11, 044057 (2019).
  * Wang  _et al._ (2018) Y. Wang, B. Yousefzadeh, H. Chen, H. Nassar, G. Huang, and C. Daraio, Physical review letters 121, 194301 (2018).
  * Wang  _et al._ (2020) Y.-F. Wang, Y.-Z. Wang, B. Wu, W. Chen, and Y.-S. Wang, Applied Mechanics Reviews 72 (2020).
  * Chen  _et al._ (2018) T. Chen, O. R. Bilal, K. Shea, and C. Daraio, Proceedings of the National Academy of Sciences 115, 5698 (2018).
  * Rothemund  _et al._ (2018) P. Rothemund, A. Ainla, L. Belding, D. J. Preston, S. Kurihara, Z. Suo, and G. M. Whitesides, Science Robotics 3 (2018).
  * Zareei  _et al._ (2020) A. Zareei, B. Deng, and K. Bertoldi, Proceedings of the National Academy of Sciences 117, 4015 (2020).
  * Bertoldi and Boyce (2008) K. Bertoldi and M. C. Boyce, Physical Review B 78, 184107 (2008).
  * Wang  _et al._ (2013) P. Wang, J. Shim, and K. Bertoldi, Physical Review B 88, 014304 (2013).
  * Shan  _et al._ (2014) S. Shan, S. H. Kang, P. Wang, C. Qu, S. Shian, E. R. Chen, and K. Bertoldi, Advanced Functional Materials 24, 4935 (2014).
  * Shim  _et al._ (2015) J. Shim, P. Wang, and K. Bertoldi, International Journal of Solids and Structures 58, 52 (2015).
  * Raney  _et al._ (2016) J. R. Raney, N. Nadkarni, C. Daraio, D. M. Kochmann, J. A. Lewis, and K. Bertoldi, Proceedings of the National Academy of Sciences 113, 9722 (2016).
  * Li  _et al._ (2012) J. Li, J. Shim, J. Deng, J. T. Overvelde, X. Zhu, K. Bertoldi, and S. Yang, Soft Matter 8, 10322 (2012).
  * Bloch (1929) F. Bloch, Journal of Physics 52, 555 (1929).
  * Watkins and Bilal (2020) A. A. Watkins and O. R. Bilal, Frontiers in Materials 7, 410 (2020).
  * Jiao and Gonella (2020) W. Jiao and S. Gonella, Physical Review B 102, 054304 (2020).
  * Press  _et al._ (2007) W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery,  _Numerical recipes 3rd edition: The art of scientific computing_ (Cambridge university press, 2007).
  * Boechler  _et al._ (2011) N. Boechler, J. Yang, G. Theocharis, P. Kevrekidis, and C. Daraio, Journal of Applied Physics 109, 074906 (2011).
  * Wang  _et al._ (2014) P. Wang, F. Casadei, S. Shan, J. C. Weaver, and K. Bertoldi, Physical review letters 113, 014301 (2014).
  * Babaee  _et al._ (2016) S. Babaee, N. Viard, P. Wang, N. X. Fang, and K. Bertoldi, Advanced Materials 28, 1631 (2016).
  * Amendola  _et al._ (2018) A. Amendola, A. Krushynska, C. Daraio, N. M. Pugno, and F. Fraternali, International Journal of Solids and Structures 155, 47 (2018).

Supporting Information:
Exploiting Localized Transition Waves to Tune Sound Propagation in Soft
Materials

Figure S1:  Minimum energy states. Equilibrium position and orientation of
disks within a fixed magnetic boundary with (a) one, (b) two, (c) three, and
(d) four embedded magnets per disk. Polar plot of the energy per disk as a
function of the disk phase angle superimposed on the corresponding disks.

Minimum energy states: The equilibrium positions of one, two, and three disks
within a magnetic boundary are determined both numerically and experimentally.
Each disk is first positioned at a random location within its own potential
well inside the boundary. The disks move in both $x$ and $y$ directions and
rotate freely to find their lowest energy state. The numerically calculated
minimum energy state for each case is presented in Figure S1. After the disk
reaches its equilibrium position, the disk is rotated 360∘ and the change of
energy as a function of orientation angle is overlaid on top of each disk as a
polar plot:

$U=-\int_{0}^{\theta_{d}}r\times\sum_{n=1}^{N}\left\\{\Sigma F_{y}\cdot
cos(\theta_{d}+\frac{2\pi n}{N})-\Sigma F_{x}\cdot sin(\theta_{d}+\frac{2\pi
n}{N})\right\\}d\theta_{d},$ (1)

where $\theta_{d}$ is the phase angle of the disk and $\Sigma F_{x}$ and
$\Sigma F_{y}$ are the sum of forces in the $x$ and $y$ directions,
respectively, on magnet $n$.

Figure S2:  Experimental disk minimum energy states. Equilibrium position and
orientation of disks within a fixed magnetic boundary with (a) one, (b) two,
(c) three, and (d) four embedded magnets per disk.

In the case of one magnet per disk, the minimum energy state is at the center
of the potential well for one, two, or three disks. Given the symmetry of the
disk, the energy minima is independent from the disk’s phase angle (Fig. S1a).
In the case of two magnets per disk, the minimum energy state is located at
the middle of the potential well with a 90∘ (vertical) phase angle. The same
vertical orientation represents the equilibrium positions for one, two, and
three disks, separately, within the magnetic boundary (Fig. S1b). It is
important to note the change in the energy polar plot when all three disks are
present within the boundary (Fig. S1b-right). In the case of a single disk
with three embedded magnets, the minimum energy state is not located at the
middle of the potential well, but is rather slightly higher or lower than the
center point along the y-axis. When two disks with three embedded magnets are
positioned within the boundary, they retain an opposite phase angle in order
to minimize their own energy. By adding a disk at the central potential well,
(i.e., when three disks are enclosed by the boundary), there exists no unique
phase angle or location for all disks to minimize their own energy as well as
the total energy of the assembly, hence the system is geometrically frustrated
(Fig. S1c). In the case of four magnets per disk, the minimum energy state is
located at the middle of the potential well with a 0∘ (plus sign) phase angle.
The same plus sign orientation represents the equilibrium positions for one,
two, and three disks, separately, within the magnetic boundary (Fig. S1d). The
energy polar plots of all disks do not change when the third disk is added
within the boundary (Fig. S1d-right), in contrast to the case of two magnets
per disk.

Experimental methods: The fixed boundaries are laser cut out of a transparent
acrylic plastic sheet using a Full-spectrum 24 pro-series laser cutter. The
spacing between the embedded magnets (from magnet center to center) in the
boundary is 28 mm in both the $x$ and $y$ directions. The width of the channel
in the $y$ direction is 20 mm. The disks are fabricated out of VeroBlack
material using an Object 3D printer (Connex 500). The fabricated disks have a
radius $r_{disk}=7.5$ mm. Both the boundary and the disk magnets have a radius
of $r_{magnet}=1.5$ mm. The disks are placed within the boundary so that the
magnets are oriented with their north pole facing upwards. The fixed boundary
magnets have the same magnetic orientation. To minimize friction, we attach a
glass cover slide at the bottom of each disk and float them on an air bearing
(New way S1030002).

To experimentally validate the potential energy calculations in Figure S1, we
enclose one, two, and three disks of each disk type (i.e., with 1-4 magnets)
within a fixed boundary. The disks start in random positions relative to the
boundary and to each other. Then, we activate a layer of laminar air flow
beneath the surface of the disks to allow them to float freely (i.e.,
levitating the disks similar to the arcade game “Air Hockey”). As predicted
numerically, the disks with 1, 2, and 4 embedded magnets self-align in a
perfectly ordered lattice (Fig. S2).

To experimentally determine the dynamical properties of our lattices, we
excite the rightmost end of the lattice with a mechanical shaker (Brüel and
Kjær 4180) and a function generator (Keysight Technologies 33512B). The
excitation takes place as a chirp signal, ranging from 0.2 Hertz - 35 Hertz.
The motion of the disks is captured using a computer vision camera (Blackfly S
USB3) and the resulting images are analyzed using the digital image
correlation software (DICe). A fast Fourier transform (FFT) is applied to the
displacement signals and the transmitted frequency ranges are observed. The
FFT graphs depict the ranges of propagation and attenuation of mechanical
waves through the lattice. We obtain our experimental dispersion curves by
exciting the lattices using a broad frequency spectrum. Then track the motion
of discrete points along the wave propagation direction (in our case the
lattice -made out of disks- is already discrete). The resulting matrix (which
includes the displacement time-history of each particle) can be transformed to
the presented dispersion curves by applying the fft2 command in Matlab.

Figure S3:  Disk Unit Force Vectors. Location and coupling of the magnetic
forces used in determining the stiffness matrix of a lattice with disks
containing 1, 2, 3, or 4 embedded magnets.

Dispersion Analysis: The dispersion relation of each system is calculated by
solving the eigenvalue problem

$[-\omega^{2}\textbf{M}+\mathbf{K}(\bm{\kappa})]\bm{\phi}=0$ (2)

where $\omega$ is the frequency, $\kappa$ is the wavenumber,
$\bm{\phi}=[u~{}v]^{T}$ is the Bloch displacement vector in the $x$ and $y$
directions $\bm{M}$ is the mass matrix:

$\bm{\mathbf{M}}=\begin{bmatrix}m_{i}&0\\\ 0&m_{i}\end{bmatrix}$ (3)

where $m_{i}$ is the mass of the disk with $i$ embedded magnets ($i=1,...,4$).
The stiffness matrix is a function of wave number as:

$\mathbf{\textbf{K}\left(\kappa\right)}=\sum_{n=1}^{N}\Biggl{(}\sum_{\alpha=1}^{2N}\biggl{[}f_{,d}(d_{n,\alpha})\textbf{e}_{n,\alpha}\otimes\textbf{e}_{n,\alpha}\left[cos(\kappa\cdot\textbf{R}_{n,\alpha})-1\right]\\\
+\frac{f(d_{n,\alpha})}{(d_{n,\alpha})}(\textbf{I}-\textbf{e}_{n,\alpha}\otimes\textbf{e}_{n,\alpha})\left[cos(\kappa\cdot\textbf{R}_{n,\alpha})-1\right]\biggr{]}\\\
-\sum_{\beta=1}^{4}\biggl{[}f_{,d}(d_{n,\beta})\textbf{e}_{n,\beta}\otimes\textbf{e}_{n,\beta}+\frac{f(d_{n,\beta})}{(d_{n,\beta})}(\textbf{I}-\textbf{e}_{n,\beta}\otimes\textbf{e}_{n,\beta})\biggr{]}\Biggr{)}$
(4)

where $\alpha$ is the index of the magnets in neighboring unit cells (e.g.
$j+1$, $j-1$), $\beta$ is index of the boundary magnets, $N$ is the number of
magnets in a single disk, $n$ is the index of the magnet in unit cell $j$,
$\bm{R_{n,\alpha}}$ is the unit vector between magnet $n$ and $\alpha$ and
$\otimes$ is the dyadic product. Figure S3 shows a visualization of the unit
vectors between unit cells for the different disk types. The repulsive forces
between the magnets are modeled as an inverse power law $f(d)=Ad^{\gamma}$
with $f_{,d}(d)$ representing the function’s first derivative and $\gamma=-4$.
In the instance of dipole-dipole interaction,
$A=3\mu\beta^{2}/4\pi=9.617\times 10^{-11}$ where $\mu$ is the permeability of
air and $\beta$ is the magnetic moment Watkins and Bilal (2020).

Numerical Analysis. The compression of the lattice and the transformation of
the patterns within the lattice are also considered numerically. To simulate
the nucleation and propagation of defects in the different lattices, four
cases of 10 disks confined to a boundary potential are considered. The
equation of motion for each disk is integrated in time using the Verlet method
Press  _et al._ (2007). A small amount of damping is introduced to each of the
disks while an increment of 1 mm of compression over a 5 second period occurs
over 10 seconds intervals; the additional 5 seconds are included to allow for
the disks to stabilize. Depending on both the number and arrangement of
magnets within each disk, the configuration of disks in the lattice can
restructure themselves into a new ordered pattern based on the ratio of
elements to potential wells (Fig.S4).

Figure S4:  Numerical Defect Nucleation and propagation. Compression of 10
disks confined to a fixed magnetic boundary with (a) one, (b) two, (c) three,
and (d) four embedded magnets per disk.

In the case of one magnet per disk, the uncompressed system is initially
ordered. As the compression begins, defects of two disks per potential well
form at the left end and travel right. This behavior continues and eventually
fills the entire system with an alternating pattern of one and two disks per
potential well (Fig S4 a). For disks with multiple magnets, the color of the
disk corresponds to its phase angle (Fig. S4 b-d). In the case of two magnets
per disk, all disks align vertically (90∘ phase angle) shortly after the
simulation starts. As the disks are compressed from the right, a defect of two
disks occupying the same potential well nucleates. Similar to our observation
in experiments, the defect propagates down the lattice and new ones emerge at
the compression front. The defect types that emerge are similar to the ones
identified experimentally in the main text of the paper (i.e., right or left
leaning normals to the line connecting the disk centers). The defect cells can
be easily identified with the brown phase color of the two disks within the
same well. In the case of three magnets per disk, the uncompressed system is
already disordered (note the change in color between both ends of the lattice)
and only becomes ordered after a slight compression (Fig. S4c). Although the
uncompressed system is not ordered in the same way as the systems with one,
two, and four magnets, the creation and propagation of defects still occur in
a similar manner. The defect types within the numerical simulations are also
similar to those observed in the experiments. Alternating orange and blue
disks within the same well show a flat defect and gray off-center disks show
the right or left defects (Fig. S4c). In the case of four magnets per disk, an
ordered self-assembly is a natural, unforced occurrence, similar to the
assembly of disks with 1 and 2 embedded magnets. As the system is compressed,
defects form and travel down the length of the lattice. Unlike the instances
of 1, 2, and 3 magnets per disk, the defects that travel through the lattice
in the case of 4 embedded magnets consist of three disks per defect. The
defect travels through the lattice in a similar manner, with more defects
emerging as the level of compression is increased (Fig. S4d).

Figure S5:  Time-independent transition wave. Numerically obtained energy per
unit cell as the lattice is compressed for (a) one (b) two (c) three and (d)
four embedded magnets per disk.

To test the generality of the nucleation and propagation of defects apart from
boundary (edge) effects, we simulate the deformation of longer lattices
confined to the boundary potentials. We consider a lattice made out of 50
disks for each disk type (i.e., 1-4 embedded magnets). To visualize defects in
a much larger lattice, we use the energy per disk within a unit cell to locate
defects within the lattice. Unit cells with a relatively higher energy
typically contain more than 1 disk (i.e., a defect exists). We represent the
energy per unit cell (i.e., the four boundary magnets) as:

$E_{UC}=\sum_{\beta=1}^{4}\sum_{k=1}^{D}\sum_{j=1}^{N}|{Ad_{\beta
k_{j}}^{\gamma}}|$ (5)

where $D$ is the number of free-floating disks within the unit cell
boundaries, $N\in[{1,2\dots 4}]$ is the number of magnets per disk, $d_{\beta
k_{j}}$ is the distance between the position of magnet $j$ in disk $k$ and
magnet $\beta$ in the fixed boundary, $A=3\mu\beta^{2}/4\pi=9.617\times
10^{-11}$, where $\mu$ is the permeability of air and $\beta$ is the magnetic
moment. The compression of the lattices is considered for the length of 15
unit cells (i.e., 420 mm). The energy level changes are plotted in Figure 3.
Supplementary Figure S5 gives a closer look at what the energy evolution mean
physically. Each red horizontal line corresponds to the respective set of
thirteen unit cells within the lattice (S5(i-iii)). In each case, as the
defect begins to travel down the lattice, the calculated energy corresponds to
the number of disks within a single potential well and the disk’s phase
angles. The blue coloring of the energy graphs corresponds to the ordered
alignment of the disks. As the disks begin to rotate, changing their phase
angles, the energy differences are apparent as the blue turns to a green,
orange, or yellow color. Supplementary Figure S5(c) shows the rightmost end of
the lattice containing 3 embedded magnets going through an alternating energy
cycle as the disks rotate in and out of their perfectly ordered state.


# PiTSBiCG: Parallel in Time Stable Bi-Conjugate Gradient Algorithm

Mohamed Kamel RIAHI  Department of Applied Mathematics, Khalifa University, PO
Box 127788, Abu Dhabi, UAE<EMAIL_ADDRESS>
###### Abstract

This paper presents a new algorithm for the parallel in time (PiT) numerical
simulation of time dependent partial/ordinary differential equations. We
propose a reliable alternative to the well know parareal in time algorithm, by
formulating the parallel in time problem algebraically and solve it using an
adapted Bi-Conjugate gradient stabilized method. The proposed Parallel in time
Stable Bi-Conjugate algorithm (PiTSBiCG) has a great potential in stabilizing
the parallel resolution for a variety of problems. In this work, we describe
the mathematical approach to the new algorithm and provide numerical evidences
that show its superiority to the standard parareal method.

###### keywords:

Parallel in time algorithm , BiCGStab , parareal , Acceleration, Parallel
computing , Numerical Simulation of PDEs.

††journal:

##  1 Introduction

Parallelization across the time directions has been proposed in the early
sixsty[1], see also [2] and references therein. This new fashion of
parallelization has been brought up again, to the attention of researchers and
engineers with the first appearance of the so-called parareal algorithm by
Lions et al [3]. Since then, many variants of the of the parareal algorithm
have been proposed, tested, and analyzed thoroughly. Motivated by the
flexibility that such algorithm proposes, many applications have attracted the
use of the parallel in time method; Among which, Fluid mechanics [4],
Neutronics [5], Optimal control [6, 7], Quantum physics [8] finance [9] etc…
Interested reader my refer to the recent review [1] and references therein.

The parareal algorithm benefits from a predictor-corrector scheme to gain
acceleration in a shorter time compared to serial computation. In such scheme
the prediction is performed by a computationally cheap coarse solver
$\mathcal{G}_{\triangle}$ while the correction is performed using an accurate
high ordered but computationally expensive fine solver
$\mathcal{F}_{\triangle}$. If the problem is formulated over an interval
$[0,T]$ and $(T_{n})_{n=0}^{N}$ is a set of times in this interval, then the
parareal in time method aims at building a sequence $(\lambda_{n}^{k})_{k}$
for each time $T_{n}$ such that, as $k$ goes to infinity, $\lambda_{N}^{k}$
converges to the fine solution at time $T_{n}$. As it will be outlined in
section 2, the sequence $(\lambda_{n}^{k})_{k}$ is defined through a recursion
formula involving predictions of $\mathcal{G}_{\triangle}$ and corrections of
$\mathcal{F}_{\triangle}$. The method decomposes the time domain in the sense
that it allows to divide the propagation of the fine solver over $[0,T]$ into
propagation on sub-intervals that can be run concurrently and simultaneously
on several processors.

##  2 Parallel in time settings

The main goal behind time parallelization algorithms is to be able to solve a
given time dependent problem over a predefined set of sub-intervals. Ideally,
one should be able to solve only on these sub-intervals. Unfortunately, and
because the sequential nature of the time-evolution problem, this is not
possible unless at least these sub-intervals inter-communicate information.

For a positive time $T>0$ and a bounded domain $\Omega$ with Lipschitz
boundary $\partial\Omega$, we consider the following Cauchy problem

$\begin{cases}\dfrac{\partial}{\partial
t}y(t,x)-\mathcal{A}(t,y(t,x))=f(t,x)&\hbox{ in }[0,T]\times\Omega\\\
y(0,x)=y_{0}.&\hbox{ at }\\{0\\}\times\Omega\end{cases}$ (1)

for which $\mathcal{A}(t,\bullet(t,x))$ stands for a differential linear
operator, and $f(t,x)$ is a given source term. The Cauchy problem (1) is
supplemented with appropriate boundary condition at $\partial\Omega$. It is
assumed that (1) is well posed and has a unique solution on its computational
domain. Therefore, it generates a semi-group of evolution that we shall denote
by $\mathcal{F}_{\triangle}$, in such a way $y(t+\Delta
t,x)=\mathcal{F}_{\triangle}(y(t,x))$ represents the solution at time
$t+\Delta t$ from a given initial condition $y(t,x)$ at time $t$. Without loss
of generality, we assume that the time interval is equality split into $N$
sub-intervals. The aim of parallel in time computing is solving the problem
(1) on sub-intervals $[T_{n},T_{n+1}]$ for $n=0,\cdots,N-1$, with
$0=T_{0}<T_{1}<\cdots<T_{n}<T_{n+1}:=T_{n}+\Delta t<\cdots<T_{N}:=T$. To this
end, we attribute $y_{n}(t,x)$ to the sub-interval $[T_{n},T_{n+1}]$, which
represents solution to

$\begin{cases}\dfrac{\partial}{\partial
t}y_{n}(t,x)-\mathcal{A}(t,y_{n}(t,x))=f_{n}(t,x)&\hbox{ in
}[T_{n},T_{n+1}]\times\Omega\\\ y_{n}(T_{n},x)=\lambda_{n}&\hbox{ on
}\\{T_{n}\\}\times\Omega.\end{cases}$ (2)

For a given initial condition, $\lambda_{n}$, and source term
$f_{n}=f_{|_{[T_{n},T_{n+1}]}}$. Problems (2) are independent and are, hence,
solvable in a parallel fashion. Although, compared to the sequential case of
(1), the parallel solutions $(y_{n})_{n}$ has to satisfy the continuity
condition, which writes

$\lambda_{n+1}=\mathcal{F}_{\triangle}\left(y(T_{n},x)\right)=\mathcal{F}_{\triangle}\left(\lambda_{n})\right),\quad\forall
n\in\\{0,\cdots,N-1\\}$ (3)

By writing the collection of the initial conditions $({\bf\lambda}_{n})_{n}$
in a vector representation as such
$\boldsymbol{\Lambda}=({\bf\lambda}_{0},\cdots,{\bf\lambda}_{n},\cdots.{\bf\lambda}_{N-1})^{t}$,
the continuity conditions is therefore satisfied by the solution to the
following algebraic linear system.

$\underbrace{\begin{bmatrix}I&0&\cdots&\cdots&0\\\
-\mathcal{F}_{\triangle}&I&0&\cdots&\vdots\\\ 0&\ddots&\ddots&\ddots&\vdots\\\
\vdots&\ddots&-\mathcal{F}_{\triangle}&I&0\\\
0&\cdots&0&-\mathcal{F}_{\triangle}&I\end{bmatrix}}_{{\mathbfcal{F}}}\underbrace{\begin{bmatrix}\lambda_{0}\\\
\vspace{0.1in}\vdots\\\ \lambda_{n}\\\ \vdots\\\
\lambda_{N-1}\end{bmatrix}}_{\boldsymbol{\Lambda}}=\underbrace{\begin{bmatrix}\lambda_{0}\\\
\vspace{0.05in}\vdots\\\ \vdots\\\ \vdots\\\ 0\end{bmatrix}.}_{\mathbfcal{B}}$
(4)

In [3] the parareal algorithm updates the initial conditions using a
predictor-corrector scheme as follows

$\lambda_{n+1}^{k+1}=\mathcal{G}_{\triangle}(\lambda_{n}^{k+1})+\mathcal{F}_{\triangle}(\lambda_{n}^{k+1})-\mathcal{G}_{\triangle}(\lambda_{n}^{k+1})$
(5)

where $\mathcal{F}_{\triangle}$ and $\mathcal{G}_{\triangle}$ are respectively
a fine an coarse solver for the time-evolution problem at hand. It has been
shown at the early stage of the parareal algorithm that the scheme (5) is
algebraically equivalent to

$\boldsymbol{\Lambda}^{k+1}=\boldsymbol{\Lambda}^{k}+\mathbfcal{G}^{-1}\left(\mathbfcal{B}-\mathbfcal{F}\cdot\boldsymbol{\Lambda}^{k}\right)$
(6)

In the present work, we shall investigate the algebraic structure of the non-
symmetric operator-block linear system arising from the time domain
decomposition. We, indeed, propose new and robust time parallel algorithm that
outperform the parareal algorithm. The new method has also the capability of
being applied to any evolution equation, as it is the case for the parareal
method. On the other hand, our algorithm performs well with large sub-domains
compared to the parareal method.

##  3 PiTSBiCG method

The approach we follow in developing our numerical algorithm is based on the
nature of the block-non-symmetric linear system (4). It is well-known in the
computational linear algebra literature (see for instance [10]) that the most
appropriate and efficient linear solvers are BiCG/QMR and their enhanced
variant such as GMRes. Our focus goes for the BiCG which uses the Lanczos
Biorthogonalization procedure, and process the solution through projection
upon the following Krylov vector space

${\bf\mathcal{K}}_{m}=\left\\{v_{1},\mathbfcal{F}\cdot
v_{1},\cdots,\mathbfcal{F}^{m-1}\cdot v_{1}\right\\},$

and orthogonal to

${\bf\mathcal{L}}_{m}=\left\\{w_{1},(\mathbfcal{F}^{T})\cdot
w_{1},\cdots,(\mathbfcal{F}^{T})^{m-1}\cdot w_{1}\right\\},$

It is worth mentioning that, minimal residual algorithms, such as QMR and
GMRES, that are based on matrix algebra factorization, are not suitable to
solve (4) as of the operators-block-structure of the problem. In the solution
process of (4) we only consider the matrix-by-vector product., where the
operator block structures are neither assembled nor stored. Furthermore,
algorithms that do not use factorization techniques can, actually, be used to
solve the parallel in time problem (4). Although, this may potentially engage
operator-block-structure matrix transpose calculation in a least square CG-
based programming [11, 12], or with the use of the BiCG method. Nonetheless,
the CGS algorithm which is designed to avoid evaluating the transpose of the
linear system my also face accuracy challenges as such in the matrix linear
algebra, where the residual may present unacceptable high variations hence
affects the outcome in the iterative process, where substantial buildup of
rounding error may often be observed. The BiCGStab combines residual formulas
from both BiCG and CGS and has been shown to be effective in both classical
and block versions [13]. The later observation encouraged us to consider such
algorithm for the time parallel settings.

In the sequel, we adapt the well known linear algebra iterative solver
BiCGStab to the time parallelisation setting and consider solving the
“fictive111fictive: means that the linear system is neither assembled nor
stored.” linear system $\mathbfcal{F}\cdot\boldsymbol{\Lambda}=\mathbfcal{B}$,
which we precondition using the left-preconditioner operator-block matrix
$\mathbfcal{G}^{-1}\approx\mathbfcal{F}^{-1}$. Indeed, the inverse of the
matrix $\mathbfcal{F}$ writes simply

$\mathbfcal{F}^{-1}:=\begin{bmatrix}I&0&\cdots&\cdots&0\\\
\mathcal{F}_{\triangle}&I&0&\cdots&\vdots\\\
\mathcal{F}_{\triangle}^{2}&\ddots&\ddots&\ddots&\vdots\\\
\vdots&\ddots&\mathcal{F}_{\triangle}&I&0\\\
\mathcal{F}_{\triangle}^{N}&\cdots&\mathcal{F}_{\triangle}^{2}&\mathcal{F}_{\triangle}&I\end{bmatrix},\,\begin{array}[]{c}\text{
hence }\\\ \text{approximated}\\\
\text{by}\end{array},\quad\mathbfcal{G}^{-1}:=\begin{bmatrix}I&0&\cdots&\cdots&0\\\
\mathcal{G}_{\triangle}&I&0&\cdots&\vdots\\\
\mathcal{G}_{\triangle}^{2}&\ddots&\ddots&\ddots&\vdots\\\
\vdots&\ddots&\mathcal{G}_{\triangle}&I&0\\\
\mathcal{G}_{\triangle}^{N}&\cdots&\mathcal{G}_{\triangle}^{2}&\mathcal{G}_{\triangle}&I\end{bmatrix}$
(7)

Therefore, in practice we solve

$\mathbfcal{G}^{-1}\cdot\mathbfcal{F}\cdot\boldsymbol{\Lambda}=\mathbfcal{G}^{-1}\cdot\mathbfcal{B}$
(8)

without of course assembling this linear system. We rather use the matrix-by-
vector output products, following the steps of Algorithm 1.

Input: Ininital guess
$\boldsymbol{\Lambda}^{0}=({\bf\lambda}_{0},\cdots,{\bf\lambda}_{N-1})^{T}$,

tolerance $\epsilon$, restart tolerance $\epsilon_{0}$

1
$\mathbfcal{R}^{0}\leftarrow\mathbfcal{G}^{-1}(\mathbfcal{B}-\mathbfcal{F}\cdot\boldsymbol{\Lambda}^{0})$

2 $\tilde{\mathbfcal}{R}\leftarrow\mathbfcal{R}^{0}$

3 $\mathbfcal{P}^{0}\leftarrow\mathbfcal{R}^{0}$

4 $k\leftarrow 0$

5 while  _$\|\mathbfcal{R}^{k}\|_{N,\infty} >\epsilon$_ do

6
$\mathbfcal{D}^{k}\leftarrow\mathbfcal{G}^{-1}\cdot\mathbfcal{F}\cdot\mathbfcal{P}^{k}$

7
$\alpha_{k}=\langle\mathbfcal{R}^{k},\tilde{\mathbfcal}{R}\rangle_{\left(L^{2}(\Omega)\right)^{N}}/\langle\mathbfcal{D},\tilde{\mathbfcal}{R}\rangle_{\left(L^{2}(\Omega)\right)^{N}}$

8 $\mathbfcal{S}^{k}\leftarrow\mathbfcal{R}^{k}-\alpha_{k}\mathbfcal{D}^{k}$

9
$\mathbfcal{K}^{k}\leftarrow\mathbfcal{G}^{-1}\cdot\mathbfcal{F}\cdot\mathbfcal{S}^{k}$

10
$\omega_{k}=\langle\mathbfcal{K}^{k},\mathbfcal{S}^{k}\rangle_{\left(L^{2}(\Omega)\right)^{N}}/\langle\mathbfcal{K}^{k},\mathbfcal{K}^{k}\rangle_{\left(L^{2}(\Omega)\right)^{N}}$

11
$\boldsymbol{\Lambda}^{k+1}\leftarrow\boldsymbol{\Lambda}^{k}+\alpha_{k}\mathbfcal{P}^{k}+\omega_{k}\mathbfcal{S}^{k}$

12
$\mathbfcal{R}^{k+1}\leftarrow\mathbfcal{S}^{k}-\omega_{k}\mathbfcal{K}^{k}$

13
$\beta_{k}=\frac{\alpha_{k}}{\omega_{k}}\langle\mathbfcal{R}^{k+1},\tilde{\mathbfcal}{R}\rangle_{\left(L^{2}(\Omega)\right)^{N}}/\langle\mathbfcal{R}^{k},\tilde{\mathbfcal}{R}\rangle_{\left(L^{2}(\Omega)\right)^{N}}$

14
$\mathbfcal{P}^{k+1}\leftarrow\mathbfcal{R}^{k+1}+\beta_{k}\left(\mathbfcal{P}^{k}-\omega_{k}\mathbfcal{D}^{k}\right)$

15 if
_$\left|\langle\mathbfcal{R}^{k+1},\tilde{\mathbfcal}{R}\rangle_{\left(L^{2}(\Omega)\right)^{N}}\right|\leq\epsilon_{0}$_
then

16 $\tilde{\mathbfcal}{R}\leftarrow\mathbfcal{R}^{k+1}$

17 $\mathbfcal{P}^{k}\leftarrow\mathbfcal{R}^{k+1}$

18 end if

19 $k\leftarrow k+1$

20

21 end while

Result: $\boldsymbol{\Lambda}^{\star}=\boldsymbol{\Lambda}^{k+1}$

Algorithm 1 PiTSBiCG algorithm

The pseudo-code is depicted in Algorithm 1, in which we consider the PiTSBiCG
to solves the parallel across the time direction linear system in its
preconditioned version (8). The structure of the algorithm is very much close
to the linear algebraic version of the BiCGStab well-known algorithm, although
we emphasize the block-operator structure issued from the parallel across the
time direction problem. Assume we have $N$ sub-domain for the parallel in time
computation, the PiTSBiCG generates $L^{2}(\Omega)$ solutions through its
iterative process that are defined at the breakpoints. Algebraic-wise, the set
of all solutions produced, whether they are initial conditions, residuals, or
directions they are piledup into a vector form element of the vector space
$\left(L^{2}(\Omega)\right)^{N+1}$. Such vector space is endowed with the
following block-structured inner product
$\langle\cdot,\cdot\rangle_{\left(L^{2}(\Omega)\right)^{N}}=\sum_{n=0}^{N-1}\langle\cdot,\cdot\rangle_{L^{2}(\Omega)}$
. It is worth noting that the inverse symbol applied for the Matrix-operator
$\mathbfcal{G}$ is described as a preconditioner at lines 6 and 9, which means
that we solve the preconditioning resolution sequentially. This is, permitted
as the semi-group operators $\mathcal{G}_{\triangle}$ are made, by
construction, coarse and non-expansive. One can see from (7) that such inverse
can easily be solved following forward substitution and this is nothing but a
sequential resolution. We here recall that none of the operator-matrices is
stored, and only a matrix-by-vector product is involved in our algorithm. We
present in the sequel applications of our method in a variety of examples and
compare its performance with the one of the parareal algorithm.

##  4 Numerical evidences

We consider the following advection-diffusion-reaction equation

$\begin{cases}\dfrac{\partial}{\partial t}y(t,x)+{\bf u}\cdot\nabla
y(t,x)+\nabla\cdot(\mu(t,x)\nabla y)+r(t,x)y(t,x)=f(t,x)&\hbox{ in
}[0,T]\times\Omega\\\ y(t,x)=g(t,x)&\hbox{ in }[0,T]\times\partial\Omega\\\
\dfrac{\partial y}{\partial\vec{n}}(t,x)=h(t,x)&\hbox{ in
}[0,T]\times\partial\Omega\\\ y(0,x)=y_{0}&\hbox{ at }\\{0\\}\times\Omega\\\
\end{cases}$ (9)

  * 1.

Diffusion equation: in (9) we consider ${\bf u}=0$, $r(t,x)=0$, and $\mu=1$.

  * 2.

Diffusion reaction equation: in (9) we consider ${\bf u}=0$, $r(t,x)=1.5$, and
$\mu=1$.

  * 3.

Advection diffusion reaction: in (9) we consider ${\bf u}=[-y,x]$. We put
$r(t,x)=5\cdot 10^{-1}$ and reduce the diffusion term $\mu=1\cdot 10^{-1}$.

| Iterations | Matrix-vector multiplication
---|---|---
Diffusion |  |
Diffusion Reaction |  |
Advection Reaction Diffusion |  |
Figure 1: Residual decreasing in values versus number of iterations (left
column) and versus Matrix-by-vector multiplication (right column), Diffusion
PDE case (first row), Diffusion reaction PDE case (second row), and Advection
reaction diffusion PDE case (third row).

The results of our numerical simulations are depicted in Figure (1). We used
the Finite Element method for the space discretize the handled PDEs with step
size $h=2\cdot 10^{-2}$, where the computational domain $\Omega$ is taken to
be a unit square centered at the origin. We also used Backward Euler time
marching for the time discretization with step size $dt=1.10^{-3}$, where the
total physical time $T=6.4$. The initial condition $y_{0}$ is taken as a
Gaussian. These settings applies for all numerical simulation we have
conducted, using FreeFem++ [14] software.

For each type of the PDEs in (9), we consider $2^{n}$ sub-domains with
$n=2,3,\cdots,6$. These sub-domains are allocated to $2^{n}$ CPUs in a Linux
workstation. The results are presented in term of the decreasing in value of
the residual of the preconditioned linear system versus the iterations (left)
and versus the matrix-by-vector multiplication (right). The results show the
clear superiority of the proposed algorithm compared to the plain version of
the parareal algorithm. The clear out-performance is shown in both the number
of iterations and in terms of the operations of matrix-by-vector
multiplication. We note here that PiTSBiCG makes two iterations (Bi-
directions) compared to the parareal algorithm, this means if one wants to
evaluate the iterations that ensure decent in value of the residual has to
multiply by two the iteration of PiTSBiCG. For this reason, it is more
appropriate to analyze the performance in terms of the Matrix-by-vector
product.

Future consideration, as extension of the presented results, includes the
analysis of the algorithm and exploit its applicability in real-world
engineering’s problems involving differential equation whether ordinary,
partial or fractional.

##  5 Conclusion

This paper presented a novel parallel in time algorithm based on simple
investigation of the equivalent algebraic structure of the initial condition
system of equations. We showed that our method outperforms the well-known
parareal method and presented several numerical examples to support our claim.

## References

  * [1] M. J. Gander, 50 years of time parallel time integration, Contributions in Mathematical and Computational Sciences 9 (2015) 69–113.  doi:10.1007/978-3-319-23321-5_3.
  * [2] K. Burrage, Parallel and sequential methods for ordinary differential equations, Clarendon Press, 2001.
  * [3] J.-L. Lions, Y. Maday, G. Turinici, Résolution d’edp par un schéma en temps ”pararéel”, Comptes Rendus de l’Académie des Sciences - Series I - Mathematics 332 (7) (2001) 661–668.  doi:https://doi.org/10.1016/S0764-4442(00)01793-6.
  * [4] P. F. Fischer, F. Hecht, Y. Maday, A parareal in time semi-implicit approximation of the navier-stokes equations 40 (2005) 433–440.  doi:https://doi.org/10.1007/3-540-26825-1_44.
  * [5] A.-M. Baudron, J.-J. Lautard, Y. Maday, M. K. Riahi, J. Salomon, Parareal in time 3d numerical solver for the lwr benchmark neutron diffusion transient model, Journal of Computational Physics 279 (2014) 67–79.  doi:http://dx.doi.org/10.1016/j.jcp.2014.08.037.
  * [6] Y. Maday, M.-K. Riahi, J. Salomon, Parareal in Time Intermediate Targets Methods for Optimal Control Problems, Springer Basel, Basel, 2013, pp. 79–92.  doi:10.1007/978-3-0348-0631-2_5.
URL https://doi.org/10.1007/978-3-0348-0631-2_5

  * [7] M. J. Gander, F. Kwok, J. Salomon, Paraopt: A parareal algorithm for optimality systems, SIAM Journal on Scientific Computing 42 (5) (2020) A2773–A2802.  doi:10.1137/19M1292291.
URL https://doi.org/10.1137/19M1292291

  * [8] M. Riahi, J. Salomon, S. Glaser, D. Sugny, Fully efficient time-parallelized quantum optimal control algorithm, Physical Review A 93 (4) (2016).  doi:10.1103/PhysRevA.93.043410.
  * [9] G. Bal, Y. Maday, A “parareal” time discretization for non-linear pde’s with application to the pricing of an american put 23 (2002) 189–202.  doi:https://doi.org/10.1007/978-3-642-56118-4_12.
  * [10] Y. Saad, Iterative methods for sparse linear systems, Society for Industrial and Applied Mathematics, 2003.
  * [11] An adjoint based implementation of the parareal algorithm, Procedia Computer Science 9 (2012) 1021–1029, proceedings of the International Conference on Computational Science, ICCS 2012.  doi:https://doi.org/10.1016/j.procs.2012.04.110.
  * [12] F. Chen, J. S. Hesthaven, Y. Maday, A. S. Nielsen, An adjoint approach for stabilizing the parareal method, Tech. rep. (2015).
  * [13] V. Simoncini, A stabilized qmr version of block bicg, SIAM Journal on Matrix Analysis and Applications 18 (2) (1997) 419–434.  doi:10.1137/S0895479894264673.
  * [14] F. Hecht, New development in freefem++, J. Numer. Math. 20 (3-4) (2012) 251–265.
URL https://freefem.org/