Powders & Grains 2017

11institutetext: Institut fÃr Theoretische Physik I, UniversitÃt Erlangen-
NÃrnberg, StaudtstraÃe 7, 91058 Erlangen, Germany 22institutetext: School of
Engineering and Information Technology, Murdoch University, 90 South Street,
Murdoch, WA 6150, Australia 33institutetext: Institute for Multiscale
Simulation, NÃgelsbachstrasse 49b, 91052 Erlangen, Germany

# Pomelo, a tool for computing Generic Set Voronoi Diagrams of Aspherical
Particles of Arbitrary Shape

Simon Weis 11<EMAIL_ADDRESS>Philipp W. A. SchÃnhÃfer
SW and PS have contributed equally to the work described in this project11 2 2
<EMAIL_ADDRESS>Fabian M. Schaller 11 Matthias SchrÃter 33 Gerd
E. SchrÃder-Turk 11 2 2

###### Abstract

We describe the development of a new software tool, called "Pomelo", for the
calculation of Set Voronoi diagrams. Voronoi diagrams are a spatial partition
of the space around the particles into separate Voronoi cells, e.g. applicable
to granular materials. A generalization of the conventional Voronoi diagram
for points or monodisperse spheres is the Set Voronoi diagram, also known as
navigational map or tessellation by zone of influence. In this construction, a
Set Voronoi cell contains the volume that is closer to the surface of one
particle than to the surface of any other particle. This is required for
aspherical or polydisperse systems.

Pomelo is designed to be easy to use and as generic as possible. It directly
supports common particle shapes and offers a generic mode, which allows to
deal with any type of particles that can be described mathematically. Pomelo
can create output in different standard formats, which allows direct
visualization and further processing. Finally, we describe three applications
of the Set Voronoi code in granular and soft matter physics, namely the
problem of packings of ellipsoidal particles with varying degrees of particle-
particle friction, mechanical stable packings of tetrahedra and a model for
liquid crystal systems of particles with shapes reminiscent of pears.

$d(K,\vec{x})$

$\vec{x}$

K

$K$

$d(K,\vec{x})$

Figure 1: Left: The Voronoi Diagram of a system of monodisperse spheres.
Center: The Voronoi Diagram is not suitable for a system of bidisperse
spheres, as Voronoi cells are overlapping with particles. Right: The Set
Voronoi Diagram of a mixture of differently shaped objects.

The analysis of geometries and structures on a micro scale level is an
important aspect of granular and soft matter physics to attain knowledge about
many interesting properties of particle packings, including contact numbers,
anisotropy, local volume fraction, etc. RefHecke ; RefAste ; RefWang . A well-
established concept is the so called Voronoi Diagram. Here, the system is
investigated by dividing the space into separate cells in respect to the
positions of the center of the particles. A cell assigned to a certain
particle is defined as the space (or region of space) that contains all the
volume closer to the center of this specific particle than to any other one
(see figure 1 left). This partition of space, however, only yields precise
results for monodisperse spheres as the construction fails otherwise due to
morphological properties of the objects. For nonspherical or polydisperse
particles the classical Voronoi diagram is of limited usefullnes, as shown in
figure 1 (center) for a system of bidisperse spheres. A generalized version of
the Voronoi Diagram, the Set Voronoi Diagram RefSchaller2013 , also known as
navigational map RefLuchnikov  or tessellation by zone of influence RefPreteux
, has to be applied. In this case the cells contain all space around the
particle which is closer to the particle’s surface than to the surface of any
other particle. Figure 1 (right) shows the Set Voronoi Diagram of a mixture of
differently shaped particles.
Here we introduce a software tool called Pomelo, which calculates Set Voronoi
Diagrams based on the algorithm described in RefSchaller2013 . Pomelo is
particularly versatile due to its ability to generically handle any arbitrary
shape which can be described mathematically. For instance objects neither have
to be convex nor simply connected. In the first step of the algorithm (figure
2 left) a triangulation of the particle’s bounding surface is generated to
sample its shape. Pomelo offers functionality for some common particle shapes,
but the user is can also specify generic particle shapes. After the
discretisation of the surface the system is tessellated by calculating the
classical Voronoi diagram of all surface points of this triangulation. This is
shown in figure 2 (center). To get the Set Voronoi tessellation of the system,
cells belonging to points on the same particle surface are merged to a single
cell in the last step (figure 2 right). The resulting partition represents the
Set Voronoi Diagram of the system. With this algorithm, systems and mixtures
of particles of any arbitrary shape can be treated.

Figure 2: Sketch of the Set Voronoi algorithm from left to right. Generating
points on particle surface. Calculating Voronoi Diagram of surface points.
Merging cells belonging to the same particle.

##  1 Pomelo

Pomelo is a generic Set Voronoi tool written in c++11 and licensed under GPL3.
Pomelo can be downloaded, see ref. RefPomeloDownload , as well as all
instructions regarding setting up, building and using Pomelo. The system
requirements are g++ 4.9.2 or clang++ 3.5.0-10 or any higher version.

While Pomelo can directly handle common particle shapes (mono- and
polydisperse spheres, tetrahedra, ellipsoids and spherocylinders), it also
provides a generic mode. The latter works for any shape which surface can be
described mathematically. The following two sections will describe how to use
Pomelo in both cases. To use Pomelo in generic mode (see section 1.2), lib-lua
5.2 or higher is required.

###  1.1 Common Particle Shapes

Particle shapes that are intrinsically supported by Pomelo are mono- and
polydisperse spheres, tetrahedra, spherocylinders and ellipsoids. Pomelo comes
with a set of demos and tests.

One test case is a set of polydisperse spheres in a cubic cell. The input is a
xyzr file. Its first line is the number of particles. The second line is a
comment which contains some information (separated by comma) required by
Pomelo. This includes boundary conditions (periodic/non periodic in each
axis), box size, shrink and the number steps in $\phi$ and $\theta$ for
discretizing the sphere’s surface. Every following line describes one sphere
in the packing with its parameters (coordinates and radius). To run the demo,
call Pomelo with the command line argument -SPHEREPOLY, the path to the xyzr
file and the desired output folder.

Pomelo’s output are the vertices, faces and cells of the Set Voronoi Diagram.
The output can be set to different formats, like the POLY file format, off
(for geomview) or a gnuplot readable format for easy visualisation.

###  1.2 Generic Mode

Using Pomelo in generic mode allows to calculate the Set Voronoi Diagram for
generic particles. The input is a position file, which is a list of all
particles which are described by a set of parameters each. The parameters are
a complete description of the particles surface within the packing. The read
file is key for Pomelos versatility. It allows the user to create a surface
triangulation based on the parameters given in the position file with a script
written in the lua language RefLua . The script tells Pomelo how to create a
surface triangulation with the particles parameters as given in the position
file. It can be fully customized by the user to match any specific particle
shape. This allows the user to create even systems composed of a mixture of
diferent particles. The surface triangulation of all particles will then be
used to calculate the Set Voronoi Diagram as described above.

Pomelo’s demos include a variety of examples on how to use the read script to
handle different types of particles.

##  2 Applications of Set Voronoi Diagrams

###  2.1 Ellipsoid Packings

Figure 3: Ellipsoids ($\alpha=1.4$) in an isotopic system (left) and ther Set
Voronoi Cells (right).

To demonstrate the use of Set Voronoi Diagrams the local packing fraction
which is defined as $\Phi_{l}=\frac{v_{e}}{v_{i}}$ with $v_{i}$ being the
volume of the Set Voronoi Cell and $v_{e}=\frac{4}{3}\pi abc$ the volume of
the particle is calculated. It has been shown that the probability
distribution of local packing fractions is universal for sphere packings
RefAste  and jammed pro- and oblate ellipsoid packings RefSchaller2015 .

Experimental granular packings of triaxial ellipsoids (three distinct axis
lengths $a$, $b$, $c$) have been measured at various global packing fractions
using X-ray tomography, see figure 3. Particle positions, orientations and
size have been determined and the Set Voronoi Diagrams have been calculated
for each packing. Here the particles are shrinked by Pomelo to improve the
quality of the cells. The statistical distribution of local packing fractions
$\Phi_{l}$ for spheres ($\alpha=1.0$) and triaxial ellipsoids ($\alpha=1.1$,
$\alpha=1.4$) is shown in figure 4. To first order, the distribution can be
collapsed on a master curve, as shown in Ref. RefSchaller2015 , by subtracting
the global volume fraction $\Phi_{g}=\frac{<v_{e}>}{<v_{i}>}$ and scaling the
distribution by $\sigma(\Phi_{l})$, see figure 4. Thus the functional form of
the distribution is invariant to global packing fraction $\Phi_{g}$ and aspect
ratio $\alpha$, to first order.

Figure 4: Probability distribution of rescaled local packing fractions
$\Phi_{l}$ for packings of spheres (particle aspect ratio $\alpha=1.0$) and
triaxial ellipsoids ($\alpha=1.1$, $\alpha=1.4$). Global packing fractions
range from 0.62 to 0.65.

###  2.2 Tetrahedra

Packings of tetrahedra can also be treated with Pomelo without using the
generic mode. The data was obtained by X-ray tomography RefNeudecker . Each
tetrahedra is described by the positions of its 4 vertices. Pomelo is able to
process the data directly. This gives the advantage to treat not only
tetrahedra but pyramids in general.

The Set Voronoi tessellations of tetrahedra packings (figure 5) with different
packing fraction and contact numbers are calculated with Pomelo. The local
packing fraction is shown for those systems in figure 6.

Figure 5:  Tetrahedra of one of the experiments (left) and their Set Voronoi
cels (right).  Figure 6:  Probability distribution of rescaled local packing
fractions $\Phi_{l}$ for different systems of tetrahedra. The colorbar and the
point’s color show the global packing fraction of the system.

###  2.3 Pear Shaped Particles

Figure 7: Pear shaped by two BÃzier curves. Both are determined by the aspect
ratio $\alpha=\frac{\text{Height}}{\text{Width}}$ and the degree of tapering
$\alpha_{\theta}=\frac{\text{Height}_{\theta}}{\text{Width}}$.

The third example is the calculation of the Set Voronoi Diagrams of pear
shaped or tapered particles. The shape of these particles is described by the
aspect ratio $\alpha$ and the degree of tapering $\alpha_{\theta}$. By using
two BÃzier curves forming the bottom and top part of the pear and rotating
them around their symmetry axis the surface of the particle is generated
(figure 7) RefBarmes . A triangulation algorithm of this surface is
implemented within the read file, which allows Pomelo to process pear shaped
particles in the generic mode. The position file provides the position,
orientation, size, aspect ratio and degree of tapering for every individual
particle.

Figure 8: Pears ($\alpha=3.0,\alpha_{\theta}=3.8$) in an isotropic system
(left) and their Set Voronoi cells (right).

Systems of hard pear shaped particles ($\alpha=3.0$) with different degrees of
tapering are generated using Molecular Dynamics (see figure 8) RefBarmes . The
statistical distribution of the local packing fractions $\Phi_{l}$ for
different pear systems ($\alpha_{\theta}=3.0,\Phi_{g}=0.48$ ;
$\alpha_{\theta}=3.8,\Phi_{g}=0.50$ ; $\alpha_{\theta}=6.0,\Phi_{g}=0.54$) is
shown in figure 9. Similarly to the sphere and ellipsoid packings the
distributions of the different pear systems collapse on a master curve.
Accordingly the distribution is not only invariant to global packing fraction
$\Phi_{g}$ but also to degree of tappering $\alpha_{\theta}$ for pears.
However, the data of ellipsoids/spheres and pears do not collapse on the same
curve.

Figure 9:  Probability distribution of rescaled local packing fractions
$\Phi_{l}$ for pears ($\alpha=3.0$) with different degree of tapering and
packing fraction
$(\alpha_{\theta},\Phi_{g})\in\\{(3.0,0.48);(3.8,0.50);(6.0,0.54)\\}$. The
dashed line shows the master curve of the sphere distributions.

##  3 Outlook

We have shown that Pomelo is applicable for a variety of systems, including
spheres, ellipsoids, tetrahedra and pear shaped particles. We have illustrated
its use to calculate Set Voronoi volume distributions. However, it is
conceivable to analyse other interesting measures like Minkowski tensors, by
piping the Set Voronoi Cells into corresponding analysis tools RefSchaller2015
; RefSchroederTurk2011 ; RefSchroederTurk2013 . It is yet also unclear how
similar volume distributions of different systems are.

One way to improve the calculation of the Set Voronoi Diagram is to implement
an adaptive sampling of the particle’s surface. By comparing the distance
between two neighboring surface points to the distance to the face of the Set
Voronoi Cell it is possible to obtain an estimate on where the surface
triangulation needs a better resolution and where a coarse resolution is good
enough. This estimate will be used to change the density of the surface
triangulation of the particles.

## References

  * (1) T. Aste, T. Di Matteo, M. Saadatfar, T. J. Senden, M. SchrÃter, H. L. Swinney, EPL 79, 24003 (2007)

  * (2) M. van Hecke, J. Phys: Condens. Matter, 22, 0033101 (2010)

  * (3) P. Wang, C. Song, Y. Jin, K. Wang, H. Makse, J. Stat Mech 2010, P12005 (2010)

  * (4) F. M. Schaller, S. C. Kapfer, M. E. Evans et al., Philosophical Magazine 93, 3993-4017 (2013)

  * (5) V. Luchnikov, N. Medvedev, L. Oger, J. Troadec, Phys. Rev. E 59, 7205 (1999)

  * (6) E. Preteux, J. Math. Imaging Vision 1, 239 (1992)

  * (7) `http://theorie1.physik.uni-erlangen.de/`
`research/pomelo/index.html`

  * (8) `https://www.lua.org/`

  * (9) F. M. Schaller, S. C. Kapfer, J. E. Hilton et al., EPL 111, 24002 (2015)

  * (10) M. Neudecker, S. Ulrich, S. Herminghaus, M. SchrÃter PRL 111, 028001 (2013)

  * (11) F. Barmes, M. Ricci, C. Zannoni, D. J. Cleaver, Phys. Rev. E. 68, 021708 (2003)

  * (12) G. E. SchrÃder-Turk, W. Mickel, S. C. Kapfer et al., Adv. M. 23(22-23), 2535–2553 (2011)

  * (13) G. E. SchrÃder-Turk, W. Mickel, S. C. Kapfer et al., NJP 15, 083028 (2013)


Country | Infection Rate
---|---
Netherlands | 1.11e-6
Germany | 2.89e-6
France | 3.75e-6
Italy | 6.20e-6
Spain | 7.58e-6
Belgium | 8.03e-6
Denmark | 8.69e-6
Portugal | 3.49e-5
Switzerland | 3.49e-5
Austria | 4.05e-5

Table C.5: Objective 3: Travel rates output from the best solution from the RRHCA for model given in Figure C.6 (Infection rates in Table C.4). Columns: Connection from; Rows: Connection to. | AT | BE | CH | DE | DK | ES | FR | IT | NL | PT
---|---|---|---|---|---|---|---|---|---|---
AT |  |  |  3.553e-6 |  1.471e-6 |  |  |  |  1.086e-5 |  |
BE |  |  |  |  1.115e-6 |  |  |  5.525e-6 |  |  5.408e-7 |
CH |  3.786e-6 |  |  |  1.56e-6 |  |  |  3.389e-7 |  6.361e-6 |  |
DE |  3.390e-6 |  7.833e-6 |  3.443e-6 |  |  2.125e-6 |  |  4.267e-6 |  |  3.187e-6 |
DK |  |  |  |  1.960e-6 |  |  |  |  |  |
ES |  |  |  |  |  |  |  1.462e-6 |  |  |  1.565e-6
FR |  |  1.291e-6 |  1.425e-6 |  1.082e-5 |  |  2.630e-6 |  |  7.151e-6 |  |
IT |  2.535e-6 |  |  2.529e-6 |  |  |  |  1.75e-6 |  |  |
NL |  |  9.549e-6 |  |  2.342e-6 |  |  |  |  |  |
PT |  |  |  |  |  |  2.071e-5 |  |  |  |

###  C.5 Analysis

##### Visual analysis.

For an example see Figure C.8.

Figure C.8: P10SIR behaviour. Modelled infections for 10 Western European
countries connected by road. AT = Austria, BE = Belgium, CH = Switzerland, DE
= Germany, DK = Denmark, ES = Spain, FR = France, IT = Italy, NL = Netherlands
and PT = Portugal. Plots for countries under 100 Infected are not visible due
to the scale of the graph.

##### Model checking

For illustration we include one experiment involving model checking, see
Figure C.9.

Figure C.9: P${}_{\textrm{China}}$SIR model checking. (Top, left) Time series
traces of Infectious for the 34 provinces, with infection starting in Hubei.
The ratio between rate constants for Infect and Recover is $1:10^{4}$
($k_{infect}=1.0e^{-6}$, $k_{recover}=1.0e^{-2}$.) (Top, right) Model checking
reveals three provinces with more than one peak. (Bottom) Network of China
highlighted to indicate the geographical context of the three provinces: HB -
red, NX - green, QH - brown, XZ - blue.

##### Correlation analysis.

Table C.6 records the variables, correlation matrices are given in Figure
C.10, C.11, C.12. The correlation matrices for the unfitted models indicate
that the percent of population infected at peak is related to the population
size for P48SIR, but not for the smaller model P10SIR, perhaps because there
are more smaller countries in the smaller model. The number of connections is
related to the peak and peak pecentage in the larger model, but only to the
peak in the smaller model. Area is a partial proxy for the whole population
and is strongly related to peak and time (negatively - i.e. earlier) but not
to peak percentage in P10SIR. However in P48SIR area is related to peak and
peak percentage but not to time.

However, the correlation over the fitted data indicates that the percentage of
the population at the peak is less related to the population size, compared
with the unfitted model. Also the size of the peak was less related to length
of time from the start of infections in the fitted model compared to the
unfitted model. Finally, the time of the peak of infections was not correlated
with the number of inter-country connections in the fitted model, as opposed
to being negatively correlated in the unfitted model (more connections
implying earlier peaks), indicating that the intervention policies that
governments have taken to control the epidemic in their countries have
overriden the effect of geographical connections (i.e. closing down inter-
country travel).

Table C.6: Variables used to investigate road connected models further. Categories: (1) real world data external to model, (2) derived measure connecting external (black) to blue, (3) real world data incorporated in model, (4) data generated by model.  Category | Variable | Description
---|---|---
1 | Area.Km2 |  Area in km2 of a given country.
2 | Density | Population density of a given country (population/km2).
3 | Connections | Number of connections a country has in land borders.
| initSusc | Susceptible + Infected at time 0.
| Population | The real-world country populations.
| MapPop | The population used in the model (real world population / 1000).
4 | peak | The peak of Infected, this represents active cases.
| peakPerc | The percentage of initSusc that were infected at the peak.
| time | The number of steps into the simulation the peak occurred.
Figure C.10: Correlation matrix for P10SIR, not fitted. Green indicates a
strong positive correlation and red indicates a strong negative correlation.
Black crosses represent a non-significant correlation. Infection rates and
travel rates all identical at 1.0e-6; recovery rates 1.0e-2. Figure C.11:
Correlation matrix for P10SIR, fitted. Green indicates a strong positive
correlation and red indicates a strong negative correlation. Black crosses
represent a non-significant correlation. Figure C.12: Correlation matrix for
P48SIR, not fitted. Green indicates a strong positive correlation and red
indicates a strong negative correlation. Black crosses represent a non-
significant correlation.

##### Dendrograms.

Examples are given in Figures C.13, C.14.

Figure C.13: Dendrogram of P10SIR for West Europe, hierarchical clustering.
Height is on the x-axis. Data not normalised Figure C.14: Dendrogram of P48SIR
for Europe, hierarchical clustering. Height is on the x-axis. Data not
normalised.

###  C.6 Discussion.

Results of simulating lockdown and unlock by dynamic change of infection rates
in the SIR model; Figure C.15 \- lockdown only, Figure C.16 \- lockdown
followed by unlock.

Figure C.15: SIR model - Dynamic change of infection rates: lockdown. Output
from the SIR model in Figure 1 where infection rates were stepwise decremented
during simulation. Figure legend shows cumulative infections for different
number of decrements. Figure C.16: SIR model, dynamic change of infection rate
increased back to pre-intervention level at time step 20.

##  Appendix D List of Files Provided

##### All models given in figures (provided as Snoopy files):

  * ∙

SIR.spn \- standard SIR, see Figure 1

  * ∙

SIQR.spn \- SIR extended by quarantine, see Figure 2 (Top)

  * ∙

SIAR.spn \- SIR extend by symptomatic/asymptomatic compartments, see Figure 2
(Bottom)

  * ∙

SIR-S2.spn \- SIR-S${}_{2}^{age}$, uncoloured, see Figure 3 (Top)

  * ∙

SIR-S2_enum.colspn, SIR-S2_int.colspn \- SIR-S${}_{2}^{age}$, coloured, see
Figure 3 (Middle)

  * ∙

SIR-S20.spn \- SIR-S20, see Figure 3 (Bottom)

  * ∙

P2-SIR.spn \- P2SIR, uncoloured, see Figure 4 (Top)

  * ∙

P2-SIR.colspn \- P2SIR, coloured, see Figure 4 (Bottom)

  * ∙

network-Europe4.spn \- connectivity graph for Europe, 4 countries, uncoloured,
see Figure 5 (Top middle)

  * ∙

network-Europe4.colspn \- connectivity graph for Europe, 4 countries,
coloured, see Figure 5 (Top right)

  * ∙

P4-SIR.colspn \- P4SIR, coloured, see Figure 5 (Middle)

  * ∙

P4-SIR.spn \- P4SIR, uncoloured, see Figure C.2

  * ∙

P48-SIR-S10.colspn \- P48SIR-S10, see Figure 5 (Bottom)

  * ∙

P48-SIR-S10.andl.spn \- P48SIR-S10 unfolded, see Figure 6 (Bottom)

  * ∙

SIVR.hpn \- SIR with variant virus, see FigureC.1

  * ∙

P10-QI-SIR.colspn \- P10Q$\\_\textrm{I}$SIR, see Figure C.3 (Top)

  * ∙

P10-QSIR.colspn \- P10QSIR, see Figure C.3 (Bottom)

  * ∙

P10-SIR_Separate_Rates.cpn \- $\mathcal{CPN}$ used for parameter optimisation
- unfolded version of P10SIR, given in Figure C.6.

##### Connectivity networks (provided as CANDL files):

  * ∙

network-Europe04.candl
\- unfolding: |P|=4, |T|=8, |A|=16

  * ∙

network-Europe05.candl
\- unfolding: |P|=5, |T|=16, |A|=32

  * ∙

network-Europe10.candl
\- unfolding: |P|=10, |T|=30, |A|=60

  * ∙

network-Europe48.candl
\- unfolding: |P|=47, |T|=170, |A|=340
\- note: IS isolated, https://countrycode.org
\- unfolding time: 1.302sec

  * ∙

network-China.candl
\- unfolding: P|=33, |T|=142, |A|=286
\- note: 34 states; TW isolated
\- unfolding time: 0.283sec

  * ∙

network-USA.candl
\- unfolding: |P|=48, |T|=210, |A|=420
\- note: 50 states; AK, HI isolated
\- unfolding time: 1.643sec

##### PnSIR models (provided as CANDL files):

  * ∙

P-Europe02-SIR - Europe fragment, 2 countries

  * ∙

P-Europe04-SIR - Europe fragment, 4 countries

  * ∙

P-Europe10-SIR - Western Europe, 10 countries

  * ∙

P-Europe48-SIR - All Europe, 48 countries

  * ∙

P-China-SIR - China, 34 provinces

  * ∙

P-USA-SIR - USA, 50 states

##### Python programs:

  * ∙

RRHCA (Random Restart Hill Climbing)

  * ∙

Python Web Scraper, used to to collect COVID19 Data.

##### R code:

  * ∙

Europe.Rmd, Western_Europe.Rmd: Both these files take in csv outputs from
Snoopy, Spike to produce correlation matrices and regression analysis. The
Time and Infectious traces are extracted from the csv files and joined to real
world data (e.g. population, density, country size). Correlation matrices are
then produced. Regression analysis is also conducted along with checking the
assumptions of regression.

  * ∙

Clustering.Rmd: This code takes in one csv file (output from Snoopy, Spike)
and performs agglomerative hierarchical clustering using Euclidean distance
and complete-linkage. This code selects Time and Infected traces to perform
the clustering. The code plots a dendrogram of the clustering results for
visualisation.

#####  Spike:

  * ∙

config file + andl/candl file

  * ∙

all source files to produce Figure C.16.

##### MC2:

  * ∙

property library and Unix script to produce Figure C.9.


# Features of a nano-twist phase in the nanolayered Ti3AlC2 MAX phase

Julien Guénolé  Vincent Taupin  Maxime Vallet  Wenbo Yu  Antoine Guitton

###### Abstract

Complex intermetallic materials known as MAX phases exhibit exceptional
properties from both metals and ceramics, largely thanks to their nanolayered
structure. With high-resolution scanning transmission electron microscopy
supported by atomistic modelling, we reveal atomic features of a nano-twist
phase in the nanolayered Ti3AlC2. The rotated hexagonal single-crystal is
encompassed within basal symmetric twist interfaces similar to grain
boundaries. In particular, we show that air-oxidation at
$1000\text{\,}\mathrm{\SIUnitSymbolCelsius}$ can form a twisted phase that
leads to the formation of interfacial dislocation networks with screw
characters or to severe interfacial reconstructions. Additionally, we explore
the contribution of disclinations to the representation by continuum models of
the stress field generated by such nano-twist defect in the Ti3AlC2 bulk
phase. The occurrence of this unexpected defect is expected to impact the
physical response of this nanolayered-based material as such supports
property-by-design approaches.

††journal: Scripta Materialia

[lem3]organization=Université de Lorraine, CNRS, Arts et Métiers, LEM3,
city=Metz, postcode=57070, country=France

[labex]organization=Labex DAMAS, Université de Lorraine, city=Metz,
postcode=57070, country=France

[lms]organization=Laboratoire Mécanique des Sols, Structures et Matériaux,
CentraleSupélec, CNRS UMR 8579, Université Paris-Saclay, city=Gif-sur-Yvette,
postcode=91190, country=France [ls]organization=Laboratoire Structures,
Propriétés et Modélisation des Solides, CentraleSupélec, CNRS UMR 8580,
Université Paris-Saclay, city=Gif-sur-Yvette, postcode=91190, country=France

[cmse]organization=Center of Materials Science and Engineering, School of
Mechanical and Electronic Control Engineering, Beijing Jiaotong University,
city=Beijing, postcode=100044, country=China

Ti3AlC2 belongs to the wide family of MAX phases [1]. These ternary compounds
arise a keen interest, because they are stiff, lightweight, machinable, made
from relatively inexpensive raw materials, resistant to oxidation and thermal
shock, and capable of remaining strong up to temperatures above 1,300°C in air
[2]. More than fifty compounds are thermodynamically stable and all of them
exhibit the same range of promising properties. Because of their composition,
they were called Mn+1AXn phases ($n=$ 1 to 3, M is a transition metal, A is an
A-group element and X is nitrogen and/or carbon) [3].

As typical MAX phase, Ti3AlC2 has a nanolayered structure with an hexagonal
lattice, whose the space group is $P6_{3}/mmc$ [2, 4]. The primitive cell can
be described as a stacking of two Ti6C octahedron layers with one layer of Al
[1]. Furthermore, measurements of lattice parameters with numerous methods
reveal that Ti3AlC2 exhibits elevated high crystalline anisotropy. The $c/a$
ratio is slightly higher than 6 [5].

It is established, that MAX phases experience plastic deformation by the glide
on basal planes of a-dislocations confined between the M6X and Al layers, thus
forming pile-ups and walls [6, 5, 7]. The latter can form local disorientation
areas, known as kink bands [8, 9]. Furthermore, numerous dislocation
interactions forming networks (dipoles, reactions) and high lattice friction
have been observed [10, 11]. Note, that out-of-basal plane dislocations have
been observed as well, but they do not play a role at room temperature in
standard deformation conditions [12, 13]. At high temperature, it has been
revealed, that out-of-basal plane ¡a¿-dislocations are common events and hence
cross-slip (from basal planes to prismatic or pyramidal planes) plays a key
role in the deformation [14, 15]. This increase of available glide systems is
likely to originate the brittle-to-ductile transition of MAX phases [14]. Note
also that Frank partial ¡c¿-dislocations correlated with a diffusion mechanism
were recently reported [16]. Moreover, observations of stacking faults is
another major microstructural feature, but their role in deformation both at
room and high temperatures remains unclear [6, 13, 15].

Of all MAX phases, those containing Al, such as Ti3AlC2, are the most
resistant to oxidation [2]. During exposure of Ti–Al–C MAX phases to oxidizing
environment at high temperatures, the outward diffusion of the weakly bonded
Al atoms is much faster compared to the more covalently bonded Ti atoms [2,
17, 18, 19]. Therefore, it results in a regime, where a superficial protective
layer of Al2O3 is formed. However, TiO2 can be formed as well, leading, this
time, to a catastrophic regime. Both mechanisms are strongly depend on the
oxidation conditions and the initial microstructure [18, 19]. Despite the
description of these oxides, operating oxidation mechanisms of Ti3AlC2 and
more generally of MAX phases remain poorly documented, especially on how the
crystallographic structure is evolving. During decomposition of 312 MAX phases
(such as Ti3AlC2, Ti3SiC2) at high temperatures in an oxygen-containing
atmosphere or in species having a high affinity to A-element, A-element is
observed to de-intercalate, thus forming Ti3C2 platelets [20, 21, 22, 23].
More precisely, in case of Ti2AlC-Cu composite, this de-intercalation of Al is
associated with a Frank partial dislocation-based mechanism [16]. This
formation of Ti${}_{\text{x}}$C${}_{\text{y}}$ platelets is then viewed as
reinforcement of the composites [20].

Such is the background of the current letter. This work is motivated by the
analysis of original phases that were observed after oxidation in a Ti3AlC2
matrix by high-resolution scanning transmission electron microscopy (HR-STEM).
Such nanolamellar phases appear to be twisted with respect to the surrounding
matrix and also associated with diffusion mechanisms in its neighborhood.
Hence, by interpreting our HR-STEM observations using molecular dynamics
simulations of model nano-twist phases, we study their atomic scale features
including interfacial crystal defects, excess energy, and mechanical fields.
Concerning the latter, a continuous description of the twist phase using a
disclination based mechanical framework is conducted complementary to
atomistic simulations to further evidence the complexity of such phases.

The specimen was prepared by hot isotropic pressing (HIP). Briefly, powders of
Ti, Al and TiC were mixed in stoechiometric proportions 2Ti:1.05Al:0.85C and
then cold-compacted into cylindrical steel dies using an uniaxial pressure.
Powder compacts were encapsulated into pyrex containers under high vacuum for
reactive sintering in the HIP machine. Afterwards, the specimen were oxidized
at $1000\text{\,}\mathrm{\SIUnitSymbolCelsius}$ for $25\text{\,}\mathrm{h}$ in
a flowing air atmosphere. More details are given in [24, 19]. Cross-sectional
TEM samples were prepared by focused ion beam (FIB) on a FEI Helios dual beam
Nanolab 600i. Atomically-resolved high-angle annular dark field scanning TEM
(HAADF-STEM) was performed on an Titan3 G2 S/TEM fitted with a double
aberration corrector for probe and image corrections and operating at
$300\text{\,}\mathrm{kV}$.

HR-STEM in HAADF mode of the lamella is shown in Fig. 1 and at lower
magnification in Fig. SM1 in supplementary materials. As HAADF detector senses
a greater signal from atoms with a higher atomic number Z, Ti columns appear
brighter in the resulting micrograph [25]. The nano-layered structure of
Ti3AlC2 highlighted in the circular inset is consistent with the expected one
for a 312-MAX phase projected along [-1 -1 2 0] _i.e._ zigzag stacking of
three TiC planes followed by one Al plane, the Al plane being the axis of
symmetry of the zigzag. Several large Ti${}_{\text{x}}$C${}_{\text{y}}$ laths
are clearly visible as well. These sub-phases are probably similar to Ti3C2
phases as they are formed by the diffusion of the Al interlayers of the
Ti3AlC2 phase. However, the experiments we used are not able to characterize
them precisely, and the mechanisms responsible for their formation is out of
the scope of this letter. The region of interest shows blurred atomic layers
(ROI, white square in Fig. 1). The blur is along the [1 -1 0 0] direction,
localized within the basal plane. This indicates that this part of the lattice
has undergone a rotation around the [0 0 0 1] direction. The boundaries
between this rotated phase and the other phase are marked with dashed white
lines. Particular orientation relationships for MX platelets in MAX interfaces
have been reported in the literature [26, 27]. However, such orientations will
not produce the high-resolution contrast we observe in Fig.1, in particular,
the blurred atomic layer within the NTP. Note, that the upper extremity of the
ROI shows a large, blurred area and visible inter-diffusion between Al and Ti
planes. Such mechanism that induce severe local lattice strain might be
responsible for the difference in the contrast of the
Ti${}_{\text{x}}$C${}_{\text{y}}$ phase on both side of the NTP (Fig. 1).

Figure 1:  Experimental observation of the nano-twist phase (region of
interest, ROI) in the nanolayered Ti3AlC2. Filtered HR-STEM micrography in
HAADF mode with electron direction along [1 1 -2 0]. Boundaries of the nano-
twist phase are indicated by dashed lines. A model of the Ti3AlC2
crystallographic structure is shown in the inset. Some
Ti${}_{\text{x}}$C${}_{\text{y}}$ similar to Ti3C2 laths are also observed.
Figure 2:  Modeling of the nano-twist phase in the nano-layered Ti3AlC2 by
means of a cookie-shape rotation region. (a) atomistic simulation setup of the
cookie-shape region (in blue) parallel to the basal planes, with [0001]
rotation axis and $\psi$ rotation angle. The black square indicate the
location of the inset. : magnified projected view of a
$10\text{\,}\mathrm{nm}$ thin slice of the side of the cookie-shape nano-twist
phase for $\Psi=$$$. Ti, Al and C are colored in green, red and black,
respectively.

To gain more insights on the observed defect and to characterize precisely the
structure of the twist boundary, we modeled the nano-twist phase by atomistic
simulations. Molecular dynamics (MD) based on interatomic potentials is an
excellent modelling method to investigate the atomic-scale configuration of
defects, in particular interfaces [28, 29]. It has been widely used over the
past decade to investigate grain boundaries [30, 31, 32, 33] and phase
boundaries [34, 35, 36, 37]. The only published model suitable to describe
Ti3AlC2 atomic interactions is the bond-order potential (BOP) recently
adjusted by Plummer and Tucker [38]. It is based on the formalism initially
proposed by Tersoff [39] and described in detail by Albe _et al._ [40]. Fig. 2
shows the atomistic simulation setup to compute the energy and the structure
of the rotated phase boundaries. Within a fully periodic bulk Ti3AlC2 phase, a
cookie shape region in the center is rotated with an angle $\Psi$ around an
axis normal to the basal plane, as indicated in blue in Fig. 2. The inset
shows as an example a rotation of the atomic structure for $\Psi=$$$. Note
that structures with similar contrasts to HAADF-STEM ones can be obtained with
different values of $\Psi$ (Fig. 1).

It is important to mention that we do not intend here to model the exact
system observed in our experiments, for several reasons. (1) The limitation of
our experimental approaches do not reveal the exact composition of the
Ti${}_{\text{x}}$C${}_{\text{y}}$ phase. (2) The interaction potential used in
our MD simulations has been designed to model the Ti3AlC2 phase and, thus,
might be less reliable for any non-stochiometric
Ti${}_{\text{x}}$C${}_{\text{y}}$ phase. (3) The local atomic environment of
the interface, up to three TiC interlayers, is identical for both Ti3AlC2 and
Ti${}_{\text{x}}$C${}_{\text{y}}$ phases. In this context, changing the
material that forms the interfaces of the NTP in our MD simulations should
influence the value of the interfacial energy, but will have a negligible
impact on what is the focus of our work: the energetic profile and the
crystallographic structure of the interfaces. The consideration of idealized
systems or surrogate materials is common practice with atomistic simulations,
including for direct comparisons with experimental results [35, 41, 37, 16].

The as-formed nano-twist phase (NTP) exhibits prismatic and basal interfaces
with the Ti3AlC2 phase. The prismatic interface is not clearly defined in our
experimental observation, whereas the basal interface appears atomically sharp
(See Fig. 1). In this work, we thus focus on the characterisation of the NTP
basal boundary (NTB). Fig. 3 shows the interfacial energy of the NTB as
function of the twist angle $\Psi$. The energy is computed by considering a
cylindrical region in the center of the NTP, that does not encompass the
prismatic boundaries. For each twist angle, the system is statically relaxed
by using conjugate gradients and fire [42, 43] methods. The dimensions of the
box are adapted to ensure a globally stress-free system. The energy of such
as-relaxed NTB is shown with blue small dots in Fig. 3. Selected
configurations have been annealed at $600\text{\,}\mathrm{K}$ for
$100\text{\,}\mathrm{ps}$. A Nose-Hoover thermostat control the temperature
and a Nose-Hoover barostat ensure constant zero pressure. The systems are
subsequently quenched and their energies are shown in Fig. 3 with orange large
dots.

Figure 3:  Energy of the basal interface of the nano-twist phase as function
of the twist angle $\Psi$. Interfacial energy of the ground state and annealed
structure in blue and orange dots, respectively. The dotted line is a guide
for the eyes. Insets show magnified views of the stress field $\sigma_{zz}$
for $\Psi=3,13,14,15,30$.

The insets in Fig. 3 present the out-of-plane stress component $\sigma_{zz}$
of the as-relaxed NTB for different $\Psi$. By showing the distortion of the
crystal lattice at the basal boundary, it reveals the structures of the
interface. More details are presented in Fig. 4 for $\Psi=$$$, the smallest
twist angle considered in this work. The distortion of the crystal lattice at
the NTB is captured by the out-of-plane stress component $\sigma_{zz}$ (Fig.
4a,b). This evokes the signature of an interfacial dislocation network, with
the three-fold symmetry of the basal plane and the nodes of interacting
interfacial dislocations. The atomic positions for different basal interlayers
shown in Fig. 4(c)(d) confirm this observation, as the relative displacement
of the atoms is characteristic of screw dislocations.

Figure 4:  Interfacial structure of the Ti3AlC2 nano-twist phase for
$\Psi=$$$. (a) global and (b) magnified views of the local stress field
$\sigma_{zz}$ revealing the typical dislocation network of a low angle grain
boundary. (c) identical view as (b) but showing chemical species as in Fig. 2.
(d) magnified view of (c) revealing the atomic displacements induced by a
screw dislocation. (c) Atomic radius are proportional to the stacking position
in the [0001] direction: light green, red and dark green, from top to bottom.
The screw dislocation line is indicated as a blue dashed line.

The NTP we observed experimentally is not a grain on itself, but is located
within one crystallographically homogeneous area containing perfectly coherent
Ti3AlC2 and Ti${}_{\text{x}}$C${}_{\text{y}}$ phases. The phase
Ti${}_{\text{x}}$C${}_{\text{y}}$ being directly issued from the Ti3AlC2
phase, this area can be considered as a Ti3AlC2-base grain. This nano-phase
exhibits limited degree of freedom, with its twist axis and at least one
interface being precisely defined: the ¡c¿ axis and the basal planes,
respectively. It is clear that the basal interface, we called NTB, shares
similarities with pure twist grain boundaries, such as (1) an interfacial
energy evolution with a typical bell-shape, (2) an interfacial structure that
can be described by a network of screw dislocations lying on basal planes and
cross-sliping to 1st-order prismatic planes (consistent with previous
experimental observations [14]) for low twist angles, (3) a spacing between
interfacial dislocations inversely proportional to low twist angle values and
(4) a description of the interface in term of dislocations network that
vanished for high twist angles. The transition from a low angle configuration
to a high angle configuration is for $\Psi\approx$$$ (Fig. 3, insets).
However, the NTB also exhibits clear differences towards grain boundaries,
namely (1) only one macroscopic degree of freedom (DoF) instead of five for
grain boundaries and (2) interfacial reconstructions for high twist angle
leading to non-monotonous structural patterns.

Interestingly, such NTP were already reported by Drouelle _et al._ in Ti3AlC2
deformed by creep at high temperature. Indeed, they observed by conventional
TEM, highly disorganized lenticular defects with a very high density of screw
dislocations [15]. In addition, Zhang _et al._ reported low-angle twist grain
boundaries (disorientation around $0.5\text{\,}\mathrm{\SIUnitSymbolDegree}$)
in Ti3AlC2 compressed uniaxially at
$1200\text{\,}\mathrm{\SIUnitSymbolCelsius}$ [44]. Such sub-grain boundaries
are formed by screw dislocation networks originating from basal dislocation
reactions, like those predicted here in Fig. 4.a. Both studies conclude that
such features may play a role in the plasticity of Ti3AlC2.

It is worth noticing that, within the time frame considered by our atomistic
simulations, neither the interfacial energy nor the interfacial structure is
significantly altered by annealing. The configuration of the NTB predicted by
atomistic modelling appears as such stable.

Figure 5:  Shear stress fields $\sigma_{xz}$ and $\sigma_{yz}$ of a high-angle
nano-twist phase represented by (a) a dipole of twist disclination loops and
(b) an atomistic cookie-shape region. Only voxels (a) and atoms (b) with
stresses with magnitude larger than 600 MPa are shown. In purple (a) is shown
the norm of the disclination density tensor.

For high twist angles ($\Psi>$$$), the interfacial configuration exhibits
severe reconstructions. By comparison to configurations with low twist angles,
this reconstruction releases localized atomic stresses but does lead to lower
interfacial energies. Interfacial energies for tilt and twist boundaries in
crystals with cubic structures often show clear drops for particular $\Psi$
values that correspond to high symmetry misorientations. Yet, such stable high
angle configuration are not observed with the NTB investigated in this work.
This can be related to the limited DoF accessible to this boundary, as
described in the following. Grain boundaries are classically considered within
the coincidence site lattice theory (CSL), which predict periodic interfaces
for particular misorientation. As such, they are defined by five macroscopic
DoF (directions and rotation), among which some lead to particularly stable
configurations with low interfacial energy. Additionally, it is known that
considering the microscopic DoF (translation) accessible to any grain
boundaries is crucial to fully determine the most stable configurations.
However, the NTB we observed does not have so many DoF. In particular, the NTP
being confined within a well defined bulk Ti3AlC2 phase, it does not allows
the NTB for translation DoF. Allowing the NTB for microscopic DoF might thus
lead to lower energy configurations. This is out of the scope of the present
work, which is focused on the characterization of the NTP we observed
experimentally.

For large angles $\Psi$, the NTB cannot be conveniently described by
dislocations anymore as discussed above and illustrated in the insets of Fig.
3. We propose that for high-angles of rotation, the NTB can be appropriately
represented by disclination loops. As the rotational counterpart of
dislocations, disclinations are line defects which introduce a discontinuity
of the elastic rotation field, referred to as the Frank vector. The concept of
disclinations is appropriate for the description of elastic fields of tilt
boundaries with high misorientations [45] and for nanotwin microstructures
[46, 47]. Regarding the cookie shape NTP considered in this work, the stress
field due to the network of overlapping dislocation arrays can be equivalently
and more simply described by using two disclination loops delimiting the NTP.
Such representation was succesfully applied for nanolamellar twins with tilt
misorientations [47]. Here, we propose to apply this concept to twist
interfaces with high misorientations, which has never been attempted so far to
our knowledge. More details can be found as Supplementary Material.

As shown in Fig. 5(a), a NTP with misorientation $30^{\circ}$ is modelled by a
dipole of twist disclination loops with Frank vector magnitude $30^{\circ}$
about the Z ([0001]) axis. A recent field disclination mechanics framework is
used to built this loop dipole [45, 47]. The two main Cauchy stress tensor
components predicted by the model are shear stresses $\sigma_{xz}$ and
$\sigma_{yz}$. The latter are characteristic of screw dislocation out of plane
shear stresses and clearly arise from the twist rotation discontinuity applied
within the NTB.

Fig. 5(b) shows $\sigma_{xz}$ and $\sigma_{yz}$ as predicted by atomistic
simulations and evidenced a fair match with what is predicted by the
disclination model in Fig. 5(a). The similarities are (1) the morphology of
the surrounding stress field maxima with two inverse dipoles in a particular
crystallographic direction, and (2) the magnitude of these stress dipoles
around $\pm 1\text{\,}\mathrm{GPa}$. This match between atomistic and
continuum representation is however partial. The differences are (1) the
direction of stress maxima rotated by $None$, and (2) stresses within the NTB
predicted by MD but not by the disclinations. These differences clearly
originate from characteristics not considered by the disclination based model.
Typically, dislocation cores effects are only considered by atomistic
modelling and leads to $\epsilon_{zz}$ eigenstrains. Additionally, the
reconstruction of the NTB for high twist angles we observed by MD, appears to
play a crucial role in the distribution of the stress field. Disclination
based models will have to be enriched to consider such atomistic mechanisms,
but they are yet much promising to model general interfaces in complex
materials at the continuum scale.

Regarding the formation mechanism of the NTP, our preliminary observations
suggest a strong influence of the Ti-Al diffusion in the basal planes at the
onset of the NTB. The diffusion could be favoured by the interfacial
dislocation network that locally induce a variation of free volume within
basal planes.

The exact impact of a nano-twist phase on the properties of the Ti3AlC2
nanolayered structure is yet unclear. But it seems that they are not anecdotal
events, as reported in [18, 44]. Similarly to what has been observed with the
formation of Frank partial dislocations [16], the interfacial dislocation
network could limit the propagation of ¡a¿ dislocations and thus participate
to the strengthening of Ti3AlC2 phases. The prismatic interfaces of the NTP
would require more investigations as they should play an important role by
hindering the propagation of basal dislocation more efficiently than the basal
interfaces. From a more general perspective, nano-twist crystals have been
shown to exhibits peculiar optical properties [48].

Additional investigations are ongoing to gain a comprehensive understanding on
the formation mechanisms of this nano-twist phase, in order to opens up new
possibilities for tailored properties MAX phases.

## Acknowledgments

This project has received financial supports from the CNRS through the MITI
interdisciplinary programs and from the National Natural Science Foundation of
China (no. 52175284). This work was performed using HPC resources from GENCI-
TGCC (grant 2020-A0080911390 and 2021-A0100911390) and from the EXPLOR center
of the Université de Lorraine.

## References

  * [1] X. Wang, Y. Zhou, J Mater Chem 12 (3) (2002) 455–460.  doi:10.1039/b108685e.
  * [2] M. Barsoum, T. El-Raghy, Am. Sci. 89 (4) (2001) 334.  doi:10.1511/2001.4.334.
  * [3] M. W. Barsoum, Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, Germany, 2013.  doi:10.1002/9783527654581.
  * [4] M. W. Barsoum, Prog Solid State Ch 28 (1-4) (2000) 201–281.  doi:10.1016/S0079-6786(00)00006-6.
  * [5] M. W. Barsoum, M. Radovic, Ann Rev Mater Res 41 (1) (2011) 195–227.  doi:10.1146/annurev-matsci-062910-100448.
  * [6] L. Farber, M. W. Barsoum, A. Zavaliangos, T. El-Raghy, I. Levin, J. Am. Ceram. Soc. 81 (6) (2005) 1677–1681.  doi:10.1111/j.1151-2916.1998.tb02532.x.
  * [7] K. Gouriet, P. Carrez, P. Cordier, A. Guitton, A. Joulain, L. Thilly, C. Tromas, Philos. Mag. 95 (23) (2015) 2539–2552.  doi:10.1080/14786435.2015.1066938.
  * [8] M. W. Barsoum, L. Farber, T. El-Raghy, Metall. Mater. Trans. A 30 (7) (1999) 1727–1738.  doi:10.1007/s11661-999-0172-z.
  * [9] A. Guitton, S. Van Petegem, C. Tromas, A. Joulain, H. Van Swygenhoven, L. Thilly, Appl. Phys. Lett. 104 (24) (2014) 241910.  doi:10.1063/1.4884601.
  * [10] A. Guitton, A. Joulain, L. Thilly, C. Tromas, Philos. Mag. 92 (36) (2012) 4536–4546.  doi:10.1080/14786435.2012.715250.
  * [11] G.-P. Bei, A. Guitton, A. Joulain, V. Brunet, S. Dubois, L. Thilly, C. Tromas, Philos. Mag. 93 (15) (2013) 1784–1801.  doi:10.1080/14786435.2012.755272.
  * [12] C. Tromas, P. Villechaise, V. Gauthier-Brunet, S. Dubois, Philos. Mag. 91 (7-9) (2011) 1265–1275.  doi:10.1080/14786435.2010.494584.
  * [13] A. Joulain, L. Thilly, J. Rabier, Philos. Mag. 88 (9) (2008) 1307–1320.
  * [14] A. Guitton, A. Joulain, L. Thilly, C. Tromas, Sci. Rep. 4 (1) (2015) 6358.  doi:10.1038/srep06358.
  * [15] E. Drouelle, A. Joulain, J. Cormier, V. Gauthier-Brunet, P. Villechaise, S. Dubois, P. Sallot, J. Alloys 693 (2017) 622–630.  doi:10.1016/j.jallcom.2016.09.194.
  * [16] W. Yu, J. Guénolé, J. Ghanbaja, M. Vallet, A. Guitton, Scr. Mater. 191 (2021) 34–39.  doi:10.1016/j.scriptamat.2020.09.007.
  * [17] E. Drouelle, V. Brunet, J. Cormier, P. Villechaise, P. Sallot, F. Naimi, F. Bernard, S. Dubois, J. Am. Ceram. Soc. 103 (2) (2020) 1270–1280.  doi:10.1111/jace.16780.
  * [18] E. Drouelle, V. Gauthier-Brunet, J. Cormier, P. Villechaise, P. Sallot, F. Naimi, F. Bernard, S. Dubois, J. Alloys 826 (2020) 154062.  doi:10.1016/j.jallcom.2020.154062.
  * [19] W. Yu, M. Vallet, B. Levraut, V. Gauthier-Brunet, S. Dubois, J. Eur. Ceram. Soc. 40 (5) (2020) 1820–1828.  doi:10.1016/j.jeurceramsoc.2020.01.042.
  * [20] J. Zhang, J. Y. Wang, Y. C. Zhou, Acta Mater. 55 (13) (2007) 4381–4390.  doi:10.1016/j.actamat.2007.03.033.
  * [21] X. H. Wang, Y. C. Zhou, Chemistry of Materials 15 (19) (9 2003).  doi:10.1021/cm030022v.
  * [22] J. Emmerlich, D. Music, P. Eklund, O. Wilhelmsson, U. Jansson, J. M. Schneider, H. Högberg, L. Hultman, Acta Materialia 55 (4) (2 2007).  doi:10.1016/j.actamat.2006.10.010.
  * [23] M. W. Barsoum, T. El‐Raghy, L. Farber, M. Amer, R. Christini, A. Adams, Journal of The Electrochemical Society 146 (10) (10 1999).  doi:10.1149/1.1392573.
  * [24] W. Yu, V. Gauthier-Brunet, T. Cabioc’h, S. Dubois, J. Am. Ceram. Soc. 97 (7) (2014) 2308–2313.  doi:https://doi.org/10.1111/jace.12930.
  * [25] D. B. Williams, C. B. Carter, Springer US, Boston, MA, 2009.  doi:10.1007/978-0-387-76501-3.
  * [26] X. Ma, Y. Zhu, X. Wang, Y. Zhou, Philosophical Magazine 84 (28) (2004) 2969–2977.
  * [27] Z. Lin, M. Li, Y. Zhou, et al., J. Mater. Sci. Technol 23 (2) (2007) 145–65.
  * [28] L. Priester, Springer, 2013.
  * [29] C. R. Weinberger, G. J. Tucker (Eds.), Vol. 245 of Springer Series in Materials Science, Springer International Publishing, Cham, 2016.  doi:10.1007/978-3-319-33480-6.
  * [30] F. Sansoz, J. F. Molinari, Acta Mater. 53 (7) (2005) 1931–1944.  doi:10.1016/j.actamat.2005.01.007.
  * [31] T. Frolov, S. V. Divinski, M. Asta, Y. Mishin, Phys. Rev. Lett. 110 (25) (2013) 255502\.  doi:10.1103/PhysRevLett.110.255502.
  * [32] M. A. Tschopp, S. P. Coleman, D. L. McDowell, Integr. Mater. Manuf. Innov. 4 (1) (2015) 11.  doi:10.1186/s40192-015-0040-1.
  * [33] H. Zhao, L. Huber, W. Lu, N. J. Peter, D. An, F. De Geuser, G. Dehm, D. Ponge, J. Neugebauer, B. Gault, D. Raabe, Phys. Rev. Lett. 124 (10) (2020) 106102.  doi:10.1103/PhysRevLett.124.106102.
  * [34] Y. Mishin, Acta Mater. 52 (6) (2004) 1451–1467.  doi:10.1016/j.actamat.2003.11.026.
  * [35] A. Prakash, J. Guénolé, J. Wang, J. Müller, E. Spiecker, M. J. Mills, I. Povstugar, P. Choi, D. Raabe, E. Bitzek, Acta Mater. 92 (2015) 33–45.  doi:10.1016/j.actamat.2015.03.050.
  * [36] A. Vaid, J. Guénolé, A. Prakash, S. Korte-Kerzel, E. Bitzek, Materialia 7 (2019) 100355.  doi:10.1016/j.mtla.2019.100355.
  * [37] J. Guénolé, M. Zubair, S. Roy, Z. Xie, M. Lipińska-Chwałek, S. Sandlöbes-Haut, S. Korte-Kerzel, Materials & Design 202 (2021) 109572.  doi:10.1016/j.matdes.2021.109572.
  * [38] G. Plummer, G. J. Tucker, Phys. Rev. B 100 (21) (2019) 214114.  doi:10.1103/PhysRevB.100.214114.
  * [39] J. Tersoff, Phys. Rev. B 37 (12) (1988) 6991–7000.  doi:10.1103/PhysRevB.37.6991.
  * [40] K. Albe, K. Nordlund, J. Nord, A. Kuronen, Phys. Rev. B 66 (3) (2002) 352051–3520514.  doi:10.1103/PhysRevB.66.035205.
  * [41] Y. Chang, W. Lu, J. Guénolé, L. Stephenson, A. Szczpaniak, P. Kontis, A. Ackerman, F. Dear, I. Mouton, X. Zhong, S. Zhang, D. Dye, C. Liebscher, D. Ponge, S. Korte-Kerzel, D. Raabe, B. Gault, Nature Communications 10 (2019).  doi:10.1038/s41467-019-08752-7.
  * [42] E. Bitzek, P. Koskinen, F. Gähler, M. Moseler, P. Gumbsch, Phys. Rev. Lett. 97 (17) (2006) 1–4.  doi:10.1103/PhysRevLett.97.170201.
  * [43] J. Guénolé, W. G. Nöhring, A. Vaid, F. Houllé, Z. Xie, A. Prakash, E. Bitzek, Comput. Mater. Sci. 175 (2020) 109584.  doi:10.1016/j.commatsci.2020.109584.
  * [44] H. Zhang, C. Zhang, T. Hu, X. Zhan, X. Wang, Y. Zhou, Sci. Rep. 6 (1) (2016) 23943\.  doi:10.1038/srep23943.
  * [45] C. Fressengeas, V. Taupin, L. Capolungo, Int. J. Solids Struct. 51 (6) (2014) 1434–1441.  doi:https://doi.org/10.1016/j.ijsolstr.2013.12.031.
  * [46] B. Reinholz, S. Brinckmann, A. Hartmaier, B. Muntifering, W. B. Knowlton, P. Müllner, Acta Mater. 108 (2016) 197–206.  doi:https://doi.org/10.1016/j.actamat.2016.02.007.
  * [47] L. Capolungo, V. Taupin, Mater. theory 3 (1) (2019) 2\.  doi:10.1186/s41313-018-0013-9.
  * [48] H. Y. Lee, M. M. A. Ezzi, N. Raghuvanshi, J. Y. Chung, K. Watanabe, T. Taniguchi, S. Garaj, S. Adam, S. Gradečak, Nano Lett. 21 (7) (2021) 2832–2839.  doi:10.1021/ACS.NANOLETT.0C04924.