# Adaptation of the tapered double cantilever beam test for the measurement of
fracture energy and its variations with crack speed

Aditya Vasudevan***Equal contribution Thiago Melo Grabois†††Equal contribution
Guilherme Chagas Cordeiro  Stéphane Morel  Romildo Dias Toledo Filho  Laurent
Ponson<EMAIL_ADDRESS>Institut Jean Le Rond d’Alembert, (UMR 7190),
Sorbonne Université - CNRS, Paris, France.  Programa de Engenharia Civil,
COPPE, Universidade Federal do Rio de Janeiro, Rio de Janeiro, RJ, Brazil.
Laboratório de Engenharia Civil, Universidade Estadual do Norte Fluminense
Darcy Ribeiro, Campos dos Goytacazes, RJ, Brazil.  Institut de Mécanique et
d’Ingénierie (UMR 5295), Département Génie Civil et Environnement, Université
de Bordeaux, Bordeaux, France.

###### Abstract

In this work we present the design of a new test geometry inspired by the
Tapered Double Cantilever Beam (TDCB) specimen that is shown to provide an
improved characterization of the fracture properties of brittle solids. First,
we show that our new design results in an exponential increase of the specimen
compliance with crack length, leading to an extremely stable crack growth
during the test. We determine an analytical description of this behavior,
which provides a simple procedure to extract the fracture energy without
depending on finite element calculations. Validation tests are done on
polymethylmethacrylate (PMMA) specimens. We use both finite element
simulations and our analytical model to interpret the data. We find a very
good agreement between the toughness determined by both methods. The stable
nature of crack growth in our improved TDCB specimens results in a precise
control of the crack speed. This feature is employed to go one step further
and characterize the variations of toughness with crack speed. We propose an
original optimization procedure for the determination of the material
parameters characterizing the kinetic law describing the toughness rate
dependency. Overall, the approach proposed together with the newly designed
test geometry offer unprecedented possibilities for the full and accurate
characterization of the fracture behavior of brittle materials such as rocks,
sandstone, mortar etc.

###### keywords:

Fracture test geometry , toughness measurement , stable crack growth , kinetic
law.

††journal: International Journal of Solids and Structures

##  1 Introduction

Crack propagation is the main mode of failure of materials under traction. As
a result, characterizing the resistance of materials to crack growth is
crucial in engineering science for designing safe structures and equipments.
Unfortunately, measuring these properties is not straightforward at all. The
fracture energy $G_{c}$, that provides the amount of energy required to make a
crack propagate over a unit surface can be measured by mechanical tests (Lawn,
1993; Atkinson, 1987; ASTM, 2005), and more recently, from the statistical
analysis of fracture surfaces (Vernéde and Ponson, 2014; Auvray et al., 2019).
Contrary to statistical fractography that can achieve a direct measurement of
the dissipation rate $G_{c}$, mechanical test provides an indirect measurement
of the fracture energy, through the measurement of the elastic energy release
rate $G$, that corresponds to the amount of mechanical energy released by the
loaded specimen as the crack propagates over an infinitesimal length, and that
is equal to $G_{c}$ during crack growth. As a result, measuring $G_{c}$ by
mechanical test requires to track both the crack tip position and the elastic
energy stored (and so released) by the specimen during the entire test. For
that reason, toughness measurement by fracture tests remains complex and prone
to experimental errors. In addition, fracture is out-of-equilibrium by nature
so that the resistance to failure $G_{c}$ is rate-dependent, and varies with
crack speed. For certain materials, the fracture energy can vary from a factor
two or three with the crack propagation velocity ($v$) so this effect cannot
be neglected. Recent studies showed that the rate-dependency of fracture
energy plays a crucial role in setting the effective toughness of
heterogeneous materials (Chopin et al., 2018). These rate-dependent effects
can be embedded in the so-called fracture kinetic law $G_{c}(v)$ that provides
the variations of the fracture energy with the crack speed $v$. Measuring
accurately $G_{c}(v)$ requires stable fracture both in terms of crack
propagation (i.e., nearly constant crack speed) and trajectory (straight crack
growth). These can be particularly challenging in brittle solids such as
glass, rocks or mortar for which crack growth is hard to control for several
reasons, then the choice of the geometry of the fracture test that controls
the crack stability becomes essential.

Figure 1: The adapted TDCB fracture test geometry mounted on the loading
machine.

Different fracture tests have been proposed to handle these challenges. Tests
such as the Single Edge Notch Bending (SENB) and Compact Tension (CT) present
a short crack propagation length and a rather poor stability of the fracture
process (Sih, 1974), while the Double Torsion (DT) has a complex crack
geometry rendering the interpretation of mechanical data rather complex.
Alternatively, the Tapered Double Cantilever Beam (TDCB) test, which has a
tapered geometry, enables significant stability of the fracture process and is
also easily amenable to DIC measurements (Grabois et al., 2018). But so far,
it has been mainly used to study adhesively bonded joints (Marcus and Sih,
1971; Gallagher, 1971; Davalos et al., 1998) and quasi-brittle materials like
wood (Coureau et al., 2013) or mortars (Morel et al., 2003) that exhibit a
pronounced R-curve behavior providing naturally stability to the crack growth
process. Here, we take inspiration from these tests, but modify the original
TDCB geometry to make it amenable to the study of the failure behavior of
brittle solids over a wide range of crack speeds. Through the example of PMMA,
we show that the newly proposed fracture test constitutes a rather simple and
powerful mean to measure the toughness of brittle solids and its variations
with crack speed.

The classical TDCB geometry uses a straight portion in front of the taper to
facilitate the application of the load. Moreover, the ratio of the width to
the length of the specimen is generally less than 0.4. We remove these
features (see Fig. 1) and explore ratios of width to length close to unity
which results in an unexpected exponential variation of the compliance with
the crack length and hence an exponential decay of the energy release rate
with crack length under fixed imposed displacement. This behavior results in
an increased stability of the crack propagation process. As our data will
show, this allows us to measure accurately the fracture energy over a wide
range of velocities. We harness this high stability and develop a global
optimization procedure allowing to accurately determine the material
parameters characterizing the kinetic law, $G_{c}(v)$, relating toughness and
crack speed.

The paper is organized as follows: The first part presents the new design of
the TDCB fracture test. We demonstrate, via numerical simulations, how this
new design significantly enhances the stability of crack propagation. We also
present a simple analytical description of the mechanical response of the
newly designed TDCB specimen from which experimental data can be interpreted
simply for measuring toughness. The next section concerns the application of
this new fracture test to PMMA specimens allowing us to test the validity of
our analytical model and illustrate the practical implementation of our
method. The last section is dedicated to the measurement of the fracture
kinetic law. Mechanical data obtained from tests achieved at different driving
rates are treated using an original optimization method in order to measure
the material parameters, characterizing the kinetic law $G_{c}(v)$. The
application of our method to other materials in other contexts is finally
discussed.

##  2 Design of the improved TDCB fracture test

Figure 2: (a) Schematic of the sample geometry with the finite element mesh
superposed on the upper half (b) Finite element mesh at the crack tip
vicinity. Note the exponential decay of the mesh size as one gets closer to
the crack tip as shown in (b1). The mesh used for analysis contains about
$\sim 10^{5}$ nodes and $2\times 10^{5}$ degrees of freedom. (b2) and (b3)
show the experimental methodology employed to prepare the tip of the notch and
ensure a smooth crack initiation. In (b2) the crack is sharpened by a razor
blade while in (b3) a pre-crack is created to result in a sharp notch.

A schematic of the proposed test geometry is shown in Fig. 2(a). The narrow
arms of the original design of the TDCB specimen (Blackman et al., 2003) have
been removed and the height of the sample is now comparable to its length. To
measure the fracture properties, we use the compliance method that extracts
the crack length by comparing the experimentally measured compliance to a
theoretically predicted curve. For slender specimens – such as the Double
Cantilever Beam (DCB) – beam theory can be applied to obtain this theoretical
curve, but as our geometry has a large height/length ratio we use finite
element (FE) simulations to calculate it. Quadrilateral elements with linear
interpolation functions are used and the mesh size in both the $\mathrm{x}$
and $\mathrm{y}$ directions are denoted as $\mathrm{e_{x}=e_{y}}$. A typical
specimen has dimensions $\mathrm{h_{1}}=120~{}\mathrm{e_{x}}$,
$\mathrm{h_{2}}=180~{}\mathrm{e_{x}}$, $\mathrm{L}=400~{}\mathrm{e_{x}}$ ,
$\mathrm{l_{x}}=\mathrm{l_{y}}=60~{}\mathrm{e_{x}}$ and
$\mathrm{R}=12~{}\mathrm{e_{x}}$ that correspond to
$\mathrm{h_{1}}=3~{}\mathrm{cm}$, $\mathrm{h_{2}}=4.5~{}\mathrm{cm}$,
$\mathrm{L}=10~{}\mathrm{cm}$ ,
$\mathrm{l_{x}}=\mathrm{l_{y}}=1.5~{}\mathrm{cm}$ and
$\mathrm{R}=0.3~{}\mathrm{cm}$ with
$\mathrm{e_{x}=e_{y}}=250~{}\mu\mathrm{m}$. We run 2D plane stress calculation
assuming that the problem is invariant along the third direction z.

CASTEM, an open-source finite element package developed by CEA, France, is
used for our simulations. Exploiting the symmetry in the system, we simulate
only the upper half of the sample. A typical mesh configuration used in the
analysis is shown in Fig. 2. The mesh is coarse in the body with a mesh size
$\rm{e_{x}}=250\rm{\mu m}$, that reduces exponentially to $\rm{e_{x}}=1\times
10^{-11}\rm{m}$ as we reach the tip of the crack as shown in Figure 2(b). A
unit force is applied at the point P on the hole where the pins apply a
tensile force and the displacement is constrained along the Y-direction on the
ligament (OA in Fig.2(a)) to ensure pure Mode-I fracture. Crack length is
increased incrementally and the corresponding elasticity problem is solved
under plane stress and linear elastic conditions.

We extract the displacements $\mathrm{\delta_{F}}$ and $\mathrm{\delta}$ at
the point P of application of the force and at the crack mouth opening (i.e.,
location of the clip gauge in the experimental setup) (B’), respectively (see
Fig. 2(a)). We define the ratio $\mathrm{r=\delta_{F}/\delta}$ that turns out
to depend weakly on the crack length and remains close to $\mathrm{r\simeq
0.8}$. The compliance at these two locations defined as the ratio of the
displacement over the force are computed and noted $\mathrm{\lambda_{F}}$ and
$\mathrm{\lambda}$ respectively. Note that they are related by the same ratio
$r$. All the theoretical calculations are performed using the compliance at
the point of application of force. However, as we experimentally measure crack
mouth opening displacement (CMOD) at B’ using a clip gauge (see Fig. 2), the
ratio $r$ relates to both quantities. We assume a Poisson’s ratio $\nu=0.37$
that corresponds to PMMA and set the Young’s modulus $E$ to unity as it will
be subsequently corrected analytically in the expressions of the compliance
and elastic energy release rate. Note, however, that the computed value of the
compliance $\lambda_{\mathrm{F}}$ and stress intensity factor $K_{\mathrm{I}}$
in our 2D simulations depend only slightly on the Poisson’s ratio.

###  2.1 Numerical investigation of the mechanical response of the TDCB
specimen

A semi-log plot of the non-dimensional compliance obtained from the finite
element simulations is shown in Fig. 3 as a function of the normalized crack
length. We see that for a certain range of crack lengths, the dependence of
the compliance on the crack length can be described by an exponential function
such as (solid line in Fig. 3):

$\lambda_{F}=\frac{\lambda_{0}}{Eb}e^{c/c_{0}}$ (1)

where $E$ is the Young’s modulus, $b$ the thickness of the specimen while
$\mathrm{\lambda_{0}}$ and $\mathrm{c_{0}}$ are parameters of the exponential
fit that depend on the geometry of the specimen. The compliance follows an
exponential increase in a mid-range of crack lengths from
$\mathrm{L/3<c<2L/3}$, where $L$ is the specimen length along the propagation
direction. Once the crack reaches the vicinity of the free surface, the
ligament becomes small and the compliance increases faster than an exponential
and it deviates away from the exponential fit. This exponential variation of
the compliance provides immense stability to crack propagation.

Figure 3: Normalized compliance as a function of normalized crack length. The
black dots are finite element calculations while the red dots are experimental
measurements of the compliance that is obtained by 14 cycles of loading-
unloading-reloading as shown in the $F-\delta$ curve in the inset.

We then perform a systematic analysis on the scaling of the parameters of the
fit ($c_{0}$ and $\lambda_{0}$) by performing simulations for different
geometries. We choose two different geometries that have the same taper angle
$\alpha\sim 8.5^{\mathrm{o}}$, but different lengths $l_{x}$ (i.e., that
corresponds to the distance from the hole, where a pin apply the force, to the
boundary) and $h_{1}$ (i.e., that represents the height of the specimen). We
also vary the length of the specimen $L$. The first geometry has
$l_{x}=90e_{x}$ and $h_{1}=160e_{y}$, while the second geometry has
$l_{x}=90e_{x}$ and $h_{1}=120e_{y}$. The inset of Fig. 4 shows the parameter
$c_{0}$ extracted from the fit for all the simulations. It turns out that
$c_{0}$ depends only on $L$ and varies as $\mathrm{c_{0}=0.389L-0.0229}$,
where $\mathrm{c_{0}}$ and $L$ are in meter. We find the other fitting
parameter, the amplitude $\lambda_{0}$ of the exponential law, scales as
$\mathrm{c_{0}^{2}/(h_{1}l_{x})}$. This scaling is shown from the collapse of
the curves in Fig. 4 corresponding to different values of $c_{0}$, $h_{1}$ and
$l_{x}$. This behavior remains valid in the range $\mathrm{60e_{x}\leqslant
l_{x}\leqslant 90~{}e_{x}}$, $\mathrm{120~{}e_{x}\leqslant h_{1}\leqslant
160~{}e_{x}}$ and $400~{}e_{x}\leqslant L\leqslant 1200~{}e_{x}$. Using this
scaling to normalize the axis in Fig. 4, different geometries lead to the same
exponential regime in some limited range of crack lengths. In this range, the
compliance can then be approximated by

$\lambda_{F}\simeq\frac{1}{Eb}\beta\frac{c_{0}^{2}}{h_{1}l_{x}}e^{c/c_{0}}$
(2)

where $\beta~{}\simeq~{}3.0$ is a constant obtained from the numerical fit.
This expression is represented by the solid line in Fig. 4. This results in a
simple formula for the variation of compliance for different geometries,
provided the $8.5^{\mathrm{o}}$ tapered angle that remains the same. We can
now predict the mechanical energy release rate $G$ that will be subsequently
required to interpret the experimental data. $G$ can be deduced from the
compliance formula also referred to as the Irwin-Keis equation (Lawn, 1993) by
the relation:

Figure 4: Variation of the non-dimensional specimen compliance with the non-
dimensional crack length for different specimen length and geometry. The
simulations are carried out for two different geometries with the same
tapering angle $\mathrm{\alpha}$ but with different $\mathrm{l_{x}}$ and
$\mathrm{h_{1}}$. The geometry number 1 is represented by circles and has
dimensions $\mathrm{l_{x}=90e_{x}}$ and $\mathrm{h_{1}=160e_{y}}$ while
geometry number 2 is represented by triangles with dimensions
$\mathrm{l_{x}=60e_{x}}$ and $\mathrm{h_{1}=120e_{y}}$. The solid line shows
the exponential fit of the collapsed curve and corresponds to Eq. 2.

$G=\frac{F^{2}}{2b}\frac{d\lambda_{F}}{dc}=F^{2}g$ (3)

where $g$ that corresponds to the geometry dependent part of the elastic
energy release rate is given by

$g=\frac{1}{2b}\frac{d\lambda_{F}}{dc}$ (4)

Substituting Eq. 2 in Eq. 4 one obtains its expression as a function of the
geometry of the specimen:

$g=\frac{1}{2Eb^{2}}\frac{\beta c_{0}}{h_{1}l_{x}}e^{c/c_{0}}$ (5)

Our motivation to define $g$ is that it depends only on the specimen geometry,
so $G$ can be decomposed into the product of the square of the force applied
with a geometry dependent function. As a result, for a unit force, $G=g$ as
measured in our simulations where a unit force is imposed. Taking inspiration
from Eq. 5, we plot $g$ on a semi-log scale as a function of $c/c_{0}$ after
normalization by $c_{0}/(2Eb^{2}h_{1}l_{x})$. The value of $g$ is here
obtained by three independent methods: (i) The first method consists in using
Eq. 4 with the derivative of the compliance computed numerically using the
finite element solution (see e.g., Fig. 4); (ii) the second method uses the
J-Integral that computes energy flow to the tip of the crack through a closed
contour integral around the crack tip (Rice, 1968); (iii) and the final and
third method relies on the crack opening displacement. From linear elastic
fracture mechanics, we know that the crack opening profile follows $u(x)\sim
K_{I}\sqrt{x}$, so the stress intensity factor, $K_{I}$, is obtained by
fitting the computed opening profile $u(x)$. This is done in a logarithmic
representation using a linear fit of slope 1/2. The energy release rate is
then obtained from Irwin’s relation $G=K_{I}^{2}/E$. These methods are
employed for different crack lengths and they give very similar results as
shown in Fig. 5.

Figure 5: Variation of non-dimensional elastic energy release rate $g$ as a
function of normalized crack length computed from finite element simulations
using three different methods. The simulations correspond to a specimen with
$\mathrm{l_{x}=90~{}e_{x}}$, $\mathrm{L=400~{}e_{x}}$ and
$\mathrm{h_{1}=160~{}e_{x}}$.

As the compliance is known, we can now go from force constant (dead-weight)
loading used in the simulation to displacement constant (fixed-grips) loading
used in the experiments. To convert the previous expression into fixed grip
loading, we must replace $F$ in Eq. 3 by $\delta_{F}/\lambda_{F}$. Using
$r=\delta_{F}/\delta$ and the expression of the compliance $\lambda_{F}$ given
by Eq. 1, one obtains

$G=\frac{\delta_{F}^{2}}{\lambda_{F}^{2}}g=\delta^{2}\frac{Er^{2}e^{\frac{-c}{c_{0}}}}{2\lambda_{0}c_{0}}.$
(6)

We observe that G decreases exponentially fast with crack length under fixed
imposed displacement $\delta$. This is a key property of the modified TDCB
test that allows us to explore different velocities and fracture energy. For
sake of comparison, let us now consider classical fracture tests like DCB, DT
and CT. In DCB and classic TDCB, the compliance increases as the cube of crack
length so from Eq. 6, the energy release rate decays as one over fourth power
of crack length (Sih, 1971). In DT test, the compliance varies linearly with
crack length and the energy release rate varies inversely to the square of
crack length. To have a stable and controlled crack growth, it is decisive
that the energy release rate decreases as fast as possible with crack length.
The exponential decrease achieved by the TDCB specimens is faster than the one
obtained by these classical tests that display only a polynomial decrease,
thus ensuring an extremely stable crack growth as illustrated in the
following.

##  3 Experimental Results

We use a transparent thermoplastic, the so-called Polymethylmethacrylate
(PMMA) or Plexiglas, as a model brittle solid. Samples are obtained by
machining $8~{}\mathrm{mm}$ thick PMMA sheets through laser cutting. Fracture
tests are then carried out under uniaxial tension using a Shimadzu (model AG-
Xplus) universal testing machine of $10~{}kN$ maximum load capacity. In this
setup, a force gauge of $1~{}kN$ load capacity measures the tensile forces
applied to the specimen by the two steel pins located in the sample holes and
a clip gauge measures the displacements between the lips of the crack (see
Fig. 1). Experiments are controlled by the clip gauge, imposing a constant
crack opening speed. In order to ensure a stable crack propagation and avoid
overloading the sample at initiation, we employed two different methodologies
to sharpen the V-shape notch machined during the cutting process that provided
satisfying results (see Fig. 2(b2) and (b3)): (i) we use a razor blade of 0.2
mm thickness to produce propagation over approximately $1~{}\mathrm{mm}$ and;
(ii) after laser cutting a slit of 3.5 mm we pre-load the sample at a crack
mouth opening rate of 0.5 $\mathrm{\mu m/s}$ until crack initiates. The sample
is then unloaded and the new initial crack length is recorded.

Fig. 6 shows a typical Force vs. displacement curve for a test that is driven
with a clip gauge opening rate $\dot{\delta}=d\delta/dt=2.5~{}\mu m.s^{-1}$.

Figure 6: Typical load-crack opening displacement ($F-\delta$) curve for a
single test under mode I tensile fracture in a TDCB specimen of PMMA. The blue
curve is the experimentally measured curve, while the green and the red curve
are theoretically predicted by assuming a kinetic law $G_{c}(v)$ (see Section
5) that is obtained by fitting for three or one experiment respectively. The
dotted lines of respective color are the error-bars in the theoretical
prediction.

###  3.1 Experimental validation of the compliance method

Before measuring the fracture properties of the material, we first test our
approach by measuring experimentally the compliance as a function of crack
length and comparing it to to our simulation results (see Fig. 3). The
experimental compliance is obtained by unloading and reloading the sample at
some fixed crack length that is measured independently through optical means.
As the sample is fully unloaded, we do not observe any residual displacement,
thus validating the linear elastic approach (see inset of Fig. 3). The slope
of the force-displacement curve provides the compliance value after each
measurement, we propagate the crack over a small distance to measure the new
compliance value for a longer crack. Figure 3 shows that the experimental and
theoretical compliance agree very well.

###  3.2 Measurement of crack length and crack speed

The crack tip position $c(t)$ is inferred from the comparison of the
experimentally measured compliance $\lambda(t)=\delta/F$ with the one
$\lambda^{\mathrm{FE}}(c)$ computed from finite element simulations. The crack
length normalized by the full length L of the sample is plotted in Fig. 7(a)
as a function of the imposed crack mouth opening displacement $\delta$. Along
with the use of the full FE results, we also plot the crack length predicted
from the analytical formula of Eq. (6) and observe a good agreement until
$c/L\simeq 0.7$, as expected from the limited range of validity of the
analytical formula, $L/3\lesssim c\lesssim 2L/3$. The crack growth velocity
$\displaystyle{v(t)=\frac{dc}{dt}}$ is then computed and its variations is
shown in Fig. 7(b). We observe a fairly constant crack speed in agreement with
the finite element study that predicts very stable crack growth.

---

(a) (b)

(c)

Figure 7: Experimental measurement of the crack length $c$, its velocity $v$
and energy release rate $G$ computed. (a) The crack length is extracted from
the experimental compliance value and comparison either with values computed
from FE or with the analytical formula of Eq. 2; (b) measured velocity
variations as a function of normalized crack length using the compliance
method based on finite element calculations; while (c) Elastic energy release
rate vs. the normalized crack length obtained from the compliance based on
finite element calculations (in blue) and the area method (in red) described
in Sec. 3.3. Finally, the numerical predictions derived from the kinetic law
$G_{c}(v)$ with the parameters $G_{c0}=2000\pm 100$ and $\gamma=0.2\pm 0.01$
determined from the optimization method presented in Sec. 5 is represented in
green, with the dotted green line providing the uncertainty propagated from
the error bar on the parameters values $G_{c0}$ and $\gamma$.

###  3.3 Measurement of the elastic energy release rate

Once crack length is known, we then have two methods that can be used to
measure the elastic energy release rate, $G$. The first method relies on the
estimation of the energy dissipated for an incremental crack advance using the
force-displacement response. Assume that in Fig. 6 in B the crack length is
equal to $c$ and that after the time $\Delta t$, it is equal to $c+\Delta c$.
The mechanical energy release rate is then given by the energy dissipated
during this time divided by the area of the newly created fracture surface.
The energy dissipated during this time interval $\Delta t$ is equal to the
area OBB’O (see Fig. 6) and the newly created fracture surface area is $\Delta
c\,b$ where $b$ is the width of the specimen (Morel et al., 2005). Thus,
$G(t)$ is given by

$G(t)=\frac{E^{d}(t)}{\Delta c(t)b}$ (7)

where $E^{d}(t)$ is the energy dissipated that corresponds to the area OBB’O
in Fig. 6 after replacing $\delta$ by $\delta_{F}$ using the ratio $r$. In the
second method, we use directly the finite element predictions of Eq. 3 derived
in section 2.1. The energy release rate obtained by both methods is shown in
Fig. 7(c) presenting very similar results. So from now on, Eq. 3 is used to
measure the energy release rate $G$ in our experiments.

###  3.4 Prediction of the average crack speed

Our theoretical description of crack growth in TDCB specimens is used to
predict the characteristic crack speed $v_{\mathrm{m}}$ as a function of the
imposed crack opening rate $\dot{\delta}$. Assuming that the fracture energy
values remains uniform during a test, the equation of motion
$G=G_{\mathrm{c}}$ together with the expression 6 of the elastic energy
release rate provides

$v_{m}=r\sqrt{\frac{2Ec_{0}}{\lambda_{0}G_{c}}}e^{\frac{-c}{2c_{0}}}\dot{\delta}$
(8)

after derivation with respect to time. First, we notice that the crack speed
decreases exponentially with crack length. However, as the characteristic
decay length $2c_{0}\simeq 0.8L$ is relatively large, the range of crack
speeds explored during a single test is relatively modest, a feature that can
be noticed on Fig. 7(b). Furthermore, $c$ can be taken equal to a the
characteristic crack length, typically $c\simeq 60$ mm, thus providing the
characteristic crack speed explored during the test proportional to the
imposed crack opening rate $\dot{\delta}$.

Figure 8: Log-log representation of the average crack speed during the TDCB
fracture test as a function of the imposed crack opening rate. The straight
line corresponds to the prediction of Eq. 8 where the parameters
$c_{0}=21~{}\mathrm{mm}$, $\lambda_{0}=4.81$, and $r=0.8$ are measured from
finite element simulations. The value $E=1.82~{}\mathrm{GPa}$ is the PMMA
Young’s modulus. The adopted characteristic value of the PMMA fracture energy
is $G_{\mathrm{c}}=380~{}\mathrm{J.m^{-2}}$.

Figure 8 shows the variations of the average experimental crack speed for
different loading rates in the range
$0.5~{}\mu\mathrm{m}.s^{-1}\leq\dot{\delta}\leq 100~{}\mu\mathrm{m}.s^{-1}$
imposed during the experiments. Its comparison with the analytical prediction
of Eq. 8 – displayed in a logarithmic representation – shows an excellent
agreement, reflecting that the crack speed explored experimentally can be
finely tuned by changing the loading rate imposed to the specimen. In
particular, the crack velocity is about fifty times larger than the loading
rate for PMMA material parameters in this new TDCB fracture test
configuration. Thus, the application of the analytical formula 8 enables to
select some desired crack velocity by choosing the value of the crack opening
rate.

##  4 Measurement of the kinetic law $G_{\mathrm{c}}(v)$

In polymeric materials like PMMA, slow crack propagation proceeds in two
steps: the polymeric chains are first gradually aligned until full elongation
and then failure takes place (Berry, 1962). This elongation phase is a rate-
dependent process and so depends on the velocity at which the crack
propagates. As a result, the fracture energy is also a function of the cracks
velocity (Green and Pratt, 1974). The so-called $G_{c}-v$ curve, also called
the fracture kinetic law, is a material property that can be described by a
power law in polymeric materials (Maugis, 1985). The kinetic law obtained from
compiling all our experimental data is provided in Fig. 9. Interestingly just
few experiments are sufficient to capture the kinetic law over a rather
extended range of crack speeds. This can be achieved since one single test
(for example see the test shown in green in Fig. 9 that corresponds to
$\dot{\delta}=2.5~{}\mu\mathrm{m.s^{-1}}$) explores a range of crack speed
varying over about a factor three.

Here, around six TDCB fracture tests at each imposed crack opening rate
$\dot{\delta}$ were required to characterize the rate-dependency of the
fracture energy. In the following section, we describe a methodology that can
achieve a similar accuracy, but from much less fracture tests.

Figure 9: Diagram of the variations of the fracture energy with crack speed.
All the experimental data are in black circles while the three colored symbols
highlight the experimental results used for the optimization procedure,
fitting to the kinetic law $G_{c}(v)=G_{c0}v^{\gamma}$. The red and green
lines show the fit whether it has considered one or three experiments,
respectively. The dashed lines represent the error bars of the fit.

##  5 Quantitative characterization of the kinetic law $G_{c}(v)$ from a
limited number of fracture tests

###  5.1 Prediction of the mechanical response of the modified TDCB specimens
from the kinetic law

We first show here how the rate dependency of the fracture energy can be
included into our description of the failure behavior of TDCB specimens. For
polymers, the fracture energy is generally found to depend rather strongly on
crack speed (Maugis, 1985) which is consistent with our experimental
observations. To describe this variation, a two parameter kinetic law of the
form

$G_{c}(v)=G_{c0}\,v^{\gamma}$ (9)

can be suitably used. The determination of the kinetic law of a material has
many interests as it allows to make predictions on the lifetime of a structure
or it can be used to infer the microscopic failure mechanisms taking place at
the crack tip vicinity during crack growth (like stress corrosion, etc.).
However, the experimental determination of the material parameters (like here
$G_{c0}$ and $\gamma$) that characterize it is generally very challenging. We
here propose a methodology that allows for an improved determination of these
material parameters from a limited number of tests. For illustrating the
method, we consider the kinetic law of Eq. 9, and seek to determine that
material parameters $G_{c0}$ and $\gamma$.

Let’s assume first some values of these two material parameters and predict
the corresponding force-displacement response of the modified TDCB specimen
under some prescribed opening rate $\dot{\delta}$. We start with the
Griffith’s criterion that reads as

$G=G_{c}(v)\Longrightarrow\frac{1}{2b}\left(\frac{\delta(t)}{\lambda^{FE}[c(t)]}\right)\left.\frac{d\lambda^{FE}}{dc}\right|_{c(t)}=G_{\mathrm{c}}[\dot{c}(t)]$
(10)

where $\delta(t)=\dot{\delta}\,t$ and $G_{\mathrm{c}}[\dot{c}(t)]$ is given by
Eq. 9. Equation 10 is a nonlinear first order differential equation where the
crack length variation $c(t)$ is unknown. This equation can be solved
numerically using an explicit scheme with the initial condition
$c(t=0)=c_{\mathrm{ini}}$, where $c_{\mathrm{ini}}$ is the crack length at
initiation. The force-displacement response is then reconstructed using the
parametric functions $\delta(t)=\dot{\delta}\,t$ and
$F=\dot{\delta}\,t/\lambda^{FE}[c(t)]$, which can be compared with the
experimental force-displacement curve as shown in Fig. 6. Comparing the
predicted and the experimental mechanical response of the TDCB specimen, one
can thus determine the material parameters involved in the kinetic law of Eq.
9 as described in the following section.

###  5.2 Global optimization procedure for the accurate determination of the
kinetic law

To carry out the determination of the kinetic law, we choose three experiments
performed at three different loading rates that are represented in color in
Fig. 9. The three experiments are chosen so that a large range of velocities
are covered, resulting in an improved determination of the kinetic law. For
each experiment, we assume some values of the parameters $G_{\mathrm{c}}^{0}$
and $\gamma$ and compute the error

$\epsilon_{r}=\sqrt{\sum_{\delta_{i}}\left[\frac{F_{\mathrm{test}}(\delta_{i})-F_{\mathrm{exp}}(\delta_{i})}{F_{\mathrm{exp}}}\right]^{2}}$
(11)

between the numerically predicted force-displacement curve
$F_{\mathrm{test}}(\delta)$ (following the procedure described in the previous
section) and the experimentally measured response $F_{\mathrm{exp}}(\delta)$.
The procedure is then applied for multiple experiments and the cumulative
error for all the experiments is subsequently defined as

$E_{r}=\sqrt[]{\sum_{k=1}^{n}\frac{\epsilon_{k}^{2}}{n}}$ (12)

where $n$ is the total number of experiments. Figure 10 shows a 2D map of the
error $E_{r}$ in the parametric plane $(\gamma,G_{c0})$ for $n=3$. The minimum
of the error, that corresponds to the best match between the theoretically
predicted mechanical response and the experimentally measured one, is found
for $G_{c0}=2000$ and $\gamma=0.20$. The uncertainty of both parameters is
determined by one percent deviation of the forces predicted from the optimum
fit. We find $\delta G_{c0}\sim\pm 100$ and $\delta\gamma\sim\pm 0.01$ and the
error bars are plotted in dotted lines in Fig. 9. In addition, we push the
optimization procedure by considering only one experiment (the one that
corresponds to $\dot{\delta}=2.5~{}\mu\mathrm{m}.{s}^{-1}$ in Fig. 9) and we
measure $G_{c0}=1750$ and $\gamma=0.192$ and is plotted in red in Fig. 9 while
the forces are plotted in Fig. 6. We find that with this geometry a single
experiment could be sufficient to measure the kinetic law, though using at
least three experiments that drives cracks at different velocities yields
better results as it spans over a broader range.

Once the kinetic law has been determined, one can predict all the other
features of the fracture test like $c(t)$, $v(c)$ and $G(c)$. The prediction
of each of these curves is shown in solid line in Fig. 7. They are in very
good agreement with the experimentally measured behavior.

Figure 10: Error on the prediction of the force-displacement response of the
modified TDCB specimen as a function of the parameters $G_{\mathrm{c}}^{0}$
and $\gamma$ involved in the kinetic law of Eq. 9.

##  6 Conclusion

This work demonstrated how an improved TDCB test geometry enhanced the
determination of the fracture properties of brittle materials. We proposed a
new fracture test inspired from the tapered double cantilever beam specimen
that ensures steady crack growth and allows for the precise measurement of the
fracture energy. In particular, the numerical results show that the elastic
energy release rate decays exponentially fast with crack length, providing
extreme stability to the fracture process. Moreover, the experimentally
measured compliance and crack length data were in agreement to the
simulations, which validate the designed approach.

The improved TDCB geometry was used in experiments at a macroscopic scale to
determine fracture properties (e.g., crack tip positions, crack velocities,
fracture energy) of PMMA samples. We showed the ability of the experimental
configuration to calculate the sought parameters by correlating the
experimental force and displacement data to FE results. In addition, data from
different driven load velocity tests enabled to build a relationship between
the average crack velocity and the imposed opening rate (i.e.,
$v_{\mathrm{m}}\sim\dot{\delta}$), allowing to explore any desired crack
speed.

The proposed fracture test combined with a numerical resolution of the linear
elastic fracture mechanics problem was harnessed to accurately determine the
full kinetic law $G_{c}(v)$, characterizing the variations of toughness with a
large range of crack speeds. As a result, this law permits us to predict the
sought fracture test features such as $c(t)$, $v(c)$, and $G(c)$. Overall, we
endorse that this improved TDCB methodology overcomes the limitation of other
fracture test geometries by allowing a simple and accurate determination of
the material toughness and its variations with crack speed. Another route,
which is worth further investigating to validate the approach, is to consider
different brittle solids such as rocks, glasses, and cement-based materials.

## Acknowledgement

LP and AV acknowledge the support of the City of Paris through the Emergence
Program. This study was financed in part by the Coordenação de Aperfeiçoamento
de Pessoal de Nível Superior - Brasil (CAPES) - Finance Code 001. TG’s stay at
UPMC was supported by the Doctoral Sandwich Program Abroad (PDSE fellowship)
of CAPES, grant $\\#99999.009821/2014-07$. TG, GC, and RTF also acknowledge
the financial support of the Brazilian agencies FAPERJ and CNPq to this
investigation.

## References

  * ASTM (2005) ASTM.  _ASTM E 1820: Standard Test Method for Measurement of Fracture Toughness_.  ASTM International, 2005.
  * Atkinson (1987) B. K. Atkinson.  _Fracture mechanics of rock_.  Academic Press, 1987.
  * Auvray et al. (2019) N. Auvray, N. Negi, S. Trancart, and L. Ponson.  L’analyse statistique des faciès de rupture: la science de la donnée au service de l’analyse de défaillance.  _Traitements & Matériaux_, 459:47–50, 2019.
  * Berry (1962) J. P. Berry.  Fracture processes in polymeric materials III - topography of fracture surfaces of polymethylmethacrylate.  _J. Appl. Phys._ , 33:1741–1744, 1962.
  * Blackman et al. (2003) B.R.K Blackman, H Hadavinia, A.J Kinloch, M Paraschi, and J.G Williams.  The calculation of adhesive fracture energies in mode I: revisiting the tapered double cantilever beam test.  _Engineering Fracture Mechanics_ , 70(2):233 – 248, 2003.  ISSN 0013-7944.
  * Chopin et al. (2018) Julien Chopin, Aditya Bhaskar, Atharv Jog, and Laurent Ponson.  Depinning dynamics of crack fronts.  _Phys. Rev. Lett._ , 121:235501, Dec 2018.
  * Coureau et al. (2013) Jean-Luc Coureau, Stephane Morel, and Nuno Dourado.  Cohesive zone model and quasibrittle failure of wood: A new light on the adapted specimen geometries for fracture tests.  _Engineering Fracture Mechanics_ , 109(0):328 – 340, 2013.  ISSN 0013-7944.
  * Davalos et al. (1998) J.F. Davalos, P. Madabhusi-Raman, P.Z. Qiao, and M.P. Wolcott.  Compliance rate change of tapered double cantilever beam specimen with hybrid interface bonds.  _Theoretical and Applied Fracture Mechanics_ , 29(2):125 – 139, 1998.  ISSN 0167-8442.
  * Gallagher (1971) J.P Gallagher.  Experimentally determined stress intensity factors for several contoured double cantilever beam specimens.  _Engineering Fracture Mechanics_ , 3(1):27 – 43, 1971.  ISSN 0013-7944.
  * Grabois et al. (2018) Thiago Melo Grabois, Jan Neggers, Laurent Ponson, François Hild, and Romildo Dias Toledo Filho.  On the validation of integrated dic with tapered double cantilever beam tests.  _Engineering Fracture Mechanics_ , 191:311 – 323, 2018\.  ISSN 0013-7944.
  * Green and Pratt (1974) A. K. Green and P. L. Pratt.  Measurement of the dynamic fracture toughness of polymethylmethacrylate by high-speed photography.  _Eng. frac. mech._ , 6:71–80, 1974.
  * Lawn (1993) B. Lawn.  _Fracture of brittle solids_.  Cambridge University Press, 1993.
  * Marcus and Sih (1971) H.L. Marcus and George C. Sih.  A crackline-loaded edge-crack stress corrosion specimen.  _Engineering Fracture Mechanics_ , 3(4):453 – 461, 1971.  ISSN 0013-7944.
  * Maugis (1985) D. Maugis.  Subcritical crack growth, surface energy, fracture toughness, stick-slip and embrittlement.  _J. Mater. Sci._ , 20:3041–3073, 1985.
  * Morel et al. (2003) S. Morel, G. Mourot, and J. Schmittbuhl.  Influence of the specimen geometry on r-curve behavior and roughening of fracture surfaces.  _Int. J. Frac._ , 121:23–42, 2003.
  * Morel et al. (2005) S. Morel, N. Dourado, G. Valentin, and J. Morais.  Wood: A quasi-brittle material. r-curve behavior and peak load evaluation.  _Int. J. Frac._ , 131:385–400, 2005.
  * Rice (1968) J. Rice.  A path-independant integral and the approximate analysis of strain concentration by notches and cracks.  _J. Appl. Mech._ , 35:379–386, 1968.
  * Sih (1971) G. C. Sih.  A review of the threedimensional stress problem for a cracked plate.  _Int. J. Frac. Mech._ , 7(1):39–61, 1971.
  * Sih (1974) G. C. Sih.  _Handbook of stress intensity factors_.  Institute of fracture and solid mechanics, 1974.
  * Vernéde and Ponson (2014) S. Vernéde and L. Ponson.  Method for characterizing the cracking mechanism of a material from the fracture surface thereof, 2014.  French patent 1459525.


# Understanding the role of FFNs in driving multilingual behavior in LLMs

Sunit Bhattacharya and Ondřej Bojar
Institute of Formal and Applied Linguistics
Faculty of Mathematics and Physics
Charles University
<EMAIL_ADDRESS>

###### Abstract

Multilingualism in Large Language Models (LLMs) is an yet under-explored area.
In this paper, we conduct an in-depth analysis of the multilingual
capabilities of a family of a Large Language Model, examining its
architecture, activation patterns, and processing mechanisms across languages.
We introduce novel metrics to probe the model’s multilingual behaviour at
different layers and shed light on the impact of architectural choices on
multilingual processing.

Our findings reveal different patterns of multilinugal processing in the
sublayers of Feed-Forward Networks of the models. Furthermore, we uncover the
phenomenon of "over-layerization" in certain model configurations, where
increasing layer depth without corresponding adjustments to other parameters
may degrade model performance. Through comparisons within and across
languages, we demonstrate the interplay between model architecture, layer
depth, and multilingual processing capabilities of LLMs trained on multiple
languages.

##  1 Introduction

Large Language Models (LLMs) are consistently getting better at multilingual
NLP, e.g. at machine translation (MT) or performing cross-lingual tasks. But
is this enough to conclude that these models are learning representations or
patterns that generalise across multiple languages? Such an expectation is an
old one, dating back to the proposal of “universals of language” Greenberg
(1966). Previous work (see Section 2) with word embeddings and intermediate
representations of layers learnt by auto-encoding models such as BERT Rogers
et al. (2021) has shown some degree of alignment in internal representation
among languages. However, similar analysis of decoder-only Transformer models
and especially their Feed-Forward Network (FFN) layers is so far insufficient.

In this work, we focus on the FFNs of large autoregressive language models.
FFNs are characterised by their layer-by-layer processing of input
representations and make up the bulk of the total parameters of all LLMs. They
process individual tokens in parallel, with all cross-token information
transfer happening in the self-attention layers. Still, despite their
seemingly straightforward structure, their function in Transformer models has
not been fully understood yet. Although recent work has explored how the FFNs
process layer-wise intermediate representations leading to the next word
prediction in LLMs, there are still many unanswered questions. Multilingualism
(how different languages are represented inside the language models) in FFNs
is one such topic.

Figure 1: Transformer block and the structure of FFN

In this work, we build upon the view of the FFNs introduced by Bhattacharya
and Bojar (2023). We consider the lower sublayer of the FFN as _detectors_ and
the higher sublayer as _combinators_ (Figure 1) (Figure 1). Over the duration
of training, the detectors at the different layers start getting triggered by
specific patterns in the input data (Geva et al., 2021). The activation
function “selects” the important aspects (selection) and the combinators
“combine” them to emit an output which can be interpreted as a prediction of
the planned next token for that layer Belrose et al. (2023); Geva et al.
(2022); Dar et al. (2023). In fact, the final prediction in the Transformer
models might result from a process of “incremental prediction” (based on the
idea of iterative inference in Jastrzębski et al., 2017) through the layers.
Hence, during training, the model is trained to slowly convert the input
representations close to the expected output, one layer (processing step) at a
time. Thus, the intermediate activations during inference can reveal how input
representations evolve within the model, mapping input prefixes to predicted
output.

In this paper:

  * •

Refine the analysis method of Bhattacharya and Bojar (2023) and propose a
novel measure called activation flatness to identify and study activity
patterns corresponding to multiple languages in LLMs.

  * •

Describe how this activity varies across model layers and across model sizes.

  * •

Identify distinct zones of multilingual and language-specific sub-layers in
the FFNs.

We study the working of detectors and combinators across layers of four XGLM
models Lin et al. (2022) (568M, 1.7B, 2.9B, 7.5B) which are freely available,
share the same architecture (decoder-only Transformer) and are trained on the
same data comprising of 500B tokens from 30 diverse languages. We select this
particular model family because the training data curation involved up-
sampling of under-resourced languages to achieve a balanced language
representation. This makes the models perfect candidates to study multilingual
representations in LLMs.

##  2 Related Work

#### Multilingual alignment in LLMs.

There is no consensus yet on multilinugal LLMs learning universal patterns
across languages. However there is clear evidence that they learn embeddings
which have high overlap across languages, primarily between those of the same
family Doddapaneni et al. (2021). Muller et al. (2021) show that a
multilinugal BERT model can be seen as a stack of a language-specific encoder
at the lower layers and a language-agnsotic predictor in the upper layers. The
language agnosticity of the upper layers was also confirmed by Liu et al.
(2020); Pires et al. (2019). Stańczak et al. (2022) have demonstrated
significant cross-lingual overlap between neurons in multilingual auto-encoder
language models. A probing based analysis of multilinguality of autoregressive
models was done by Mueller et al. (2022). But on the question of language
models learning universal patters, recently Gurnee et al. (2024) have
demonstrated the existance of universal neurons in GPT2 models. Xie et al.
(2021) attempted to identify language-specific neurons in an encoder-decoder
based machine translation systems using Transformers.

#### Model Sparsity.

In recent literature Liu et al. (2023); Mirzadeh et al. (2023), sparsity has
been thought of as a way to increase the efficiency of models during
inference. To do that, analyses of models employing the ReLU activation
functions have been done. Some works have also explored the “ReLU-fication”
Song et al. (2024) of pretrained networks. Sparse models also should be easier
to interpret Tamkin et al. (2023).

##  3 Methodology

###  3.1 Model snapshots

In order to better understand how language models arrive at their prediction
during inference, we collect “model snapshots”, i.e. intermediate
representations of the layers of different models of the XGLM family. To do
that, we feed the language models with prefixes of sentences one word at a
time and ask them to predict the next word.

For an English sentence like “Elementary, my dear Watson.”, we first remove
all punctuations except the last i.e. period in this case. For this sentence,
we obtain 8 subwords: [‘Elementar’, ‘y’, ‘my’, ‘de’, ‘ar’, ‘Watson’,
‘.’].111Note that the subword model used in XGLM does not distinguish between
word-internal and word-final subwords, increasing the conflation of tokens
across languages (e.g. “a” being the determiner in English and a conjunction
in Czech) The prefixes (at the subword level) would thus be: [Elementar’,
‘y’], [‘Elementar’, ‘y’, ‘my’], [‘Elementar’, ‘y’, ‘my’, ‘de’, ‘ar’],
[‘Elementar’, ‘y’, ‘my’, ‘de’, ‘ar’, ‘Watson’] and [‘Elementar’, ‘y’, ‘my’,
‘de’, ‘ar’, ‘Watson’, ‘.’].

###  3.2 Snapshots from parallel test-sets

Corresponding to each prefix, we save the model snapshots of all the layers,
focusing only on the FFN layers (i.e. the detector and combinator sublayers)
and discarding everything else. Since we are explicitly interested in the
multilingual aspects of the model(s), we feed them with prefixes obtained from
sentences across 4 different WMT222https://machinetranslate.org/wmt test-sets.
Each test set contains parallel sentences in English and another language
(Czech, French, German or Hindi). Hence, for each language pair, we end up
with model snapshots for prefixes of both languages.

Thus, each prefix yields a _detector_ and _combinator_ representation for each
layer of the model in the form of a tensor of shape $(l,d)$ where $l$
indicates the number of subwords in the prefix and $d$ indicates the dimension
of the layer. So, for a layer of the model with 1024 neurons, a prefix with 4
subwords would yield a vector of shape $(4,1024)$. This notation will be
relevant later in the paper. For the most of our analyses, we extract the
representation of the last subword of the prefix assuming that it contributes
most closely to the prediction of the next word, leading us to a vector of
shape (1,1024). Later, in Section 4.2, we will use of all the prefix
representations.

###  3.3 Focus of Interest

Our primary objective is to reveal and understand “concentrated activity”
across the layers in LLMs in the context of multilingual language processing.
We do it in three levels of granularity. We start by observing the patterns of
sparsity in the models across languages. Then we dive deeper and focus on the
nature of distribution of the activations. And finally, we determine the
degree of multilinguality of the neurons through the layers by looking at
patterns of shared activation.

##  4 Observations

###  4.1 Sparsity patterns across languages

We start with a very basic analysis of the activations of detectors and
combinators of the model. The XGLM models utilize the GELU activation function
which, similarly to ReLU, allows the training to reach a sparse
representation. Thus, we start with examining if the processing of certain
languages correspond to sparser representations. In order to assess the level
of sparsity across model sizes and layers, we use activation frequency.

Given a set of model snapshots over a set of input prefixes, we define the
activation frequency for a particular neuron as the ratio of activation counts
(number of instances where the activation value was non-zero) to the total
number of prefixes. This gives an idea about the overall importance of this
neuron for the given test inputs. Activation frequency close to one indicates
that the neuron has fired for almost all prefixes. For each layer, we collect
activation frequencies of all the neurons and consider its average and, more
importantly, standard deviation. Highly varying activation frequency across
neurons of a layer indicates that some neurons are very important and some are
very unimportant for the given set of prefixes, i.e. that the overall
representation is sparse.

We present the average and standard deviation of activation frequency of
detectors across the studied model sizes in Figure 2

Figure 2: Activation frequency for detectors along with standard deviation
plotted for English, German and Hindi.

We observe that for the detectors, the average activation frequency drops near
the output layers for all languages. But it remains stable (and close to 1)
for all other layers. Also, the standard deviation of activation frequencies
increase near the input and output layers. In other words, input and output
layers have far sparser detector representations than the middle layers. We
also observe that the representations in middle layers become increasingly
dense (smaller standard deviations) with increase in model size. Finally, the
greatest drop in the activation frequency occurs for languages like English,
French and German. In Figure 2, we show results from only three languages
(i.e. English, German and Hindi) for readability.

Figure 3: Activation frequency for combinators along with standard deviation
plotted for English, German and Hindi.

For the combinators (Figure 3), we see a similar pattern, i.e. representations
are sparser close to the input and output than the middle layers.

The authors of the XGLM model claim that the training data was curated to
reflect a balanced representation of languages. Yet, it appears that lesser-
represented languages (in terms of total tokens in the training data) like
Czech and Hindi exhibit more combinator sparsity while well-represented
languages like English, French and German exhibit greater detector sparsity in
the later layers.

One possible interpretation for this observation is that owing to their
relative under-representation in the training data, languages like Czech and
Hindi learn to utilise only a subset of the total neurons in the combinators
for generation (greater sparsity). We therefore posit that over-sampling
under-resourced languages does not improve the model’s capabilities to
generate tokens in those languages.

###  4.2 Activation flatness across layers

Do neurons in FFNs show similar activation patterns across languages? How does
the distribution of the values of activation change across layers? Are they
more _peaked_ for certain languages in some layers? To answer questions like
these, we make use of a novel metric called “activation flatness”. We describe
the metric and the observations from using it, in the section below.

#### Activation flatness

We define the activation flatness for a particular layer as :

$A=\sum_{i=1}^{n}flatness(X_{i})$ (1)

where

$flatness(X)=\frac{-\sum_{i=1}^{m}S(x_{i})\log_{2}(S(x_{i}))}{m}$ (2)

Here $m$ corresponds to the number of neurons in the layer and $S(X)$ is the
neuron activation value scaled linearly to the range [0,1] within the layer,
formally defined as:

$S(x_{i})=\frac{x_{i}-min(x_{1},x_{2},...,x_{n})}{max(x_{1},x_{2},...,x_{n})-min(x_{1},...,x_{n})}$
(3)

At a high level, activation flatness measures the entropy of normalized neuron
activations in a layer. If all neurons return similar values, the entropy will
be high, making our measure of flatness high. If only a handful of neurons
fire, the entropy of these activations will be low. Thus, the activation
flatness measures if the activations for a particular layer are more _peaked_
or _uniform_. Lower flatness indicates that the activations were more peaked
at a few neurons than layers with higher flatness.

Figure 4: Activation flatness for detectors

Figure 5: Normalized activations for detectors: XGLM 2.9B for layers
1,27,45,48

The activation flatness pattern for detectors (Figure 4) shows a consistent
pattern of decreasing flatness through the layers across all models. In other
words, the activations in the detectors get more and more peaked with
increasing layer depth. We see this clearly in Figure 5 where we look at the
activations across all prefixes for detectors of some layers. The activation
patterns for layers 1 and 27 (uniform), are for instance very different from
those of layers 45 and 48 (peaked at certain neurons).

Figure 6: Activation flatness for combinators

Figure 7: Normalized activations for combinators: XGLM 2.9B for layers 1, 21,
22, and 47

The activation flatness pattern for combinators (Figure 6) shows that all
models follow a common pattern. Regardless of the number of parameters and the
total number of layers, there is a distinct decrease in the activation
flatness near layer 20 for all models. A visual inspection of Figure 7 shows
that the activations for the prefixes indeed show a pattern of certain
‘peakedness’ in layers of low activation flatness (i.e. layers 1, 21, 22 vs.
47). Also the drop in activation flatness occurs across all languages for
every model. From the detector-combinator perspective, the detector does
pattern matching from the input representation. The combinator then uses the
output of the detector (after being passed through an activation function) to
make an intermediate prediction for that layer. Any drop in the activation
flatness for combinators thus indicates more ‘focused’ firing along specific
dimensions.

Figure 8: Combinator activation for all prefixes for layer 1: 2.9B (top,
activation flatness value of $\sim$1000) vs. 7.5B (bottom, activation flatness
of $\sim$2000) models.

To better illustrate how activation flatness works, we also compare the
activation snapshots for layer 1 of the 2.9B and 7.5B parameter model in
Figure 8. We observe that for the first layer of the 2.9B parameter model,
where the combinator activation flatness value was around 1000 for all
languages, there is a single distinct peak in the activation distribution
compared to the first layer of the 7.5B parameter model (with activation
flatness of around 2000 across languages). Upon closer inspection we find that
for the first layer of the 2.9B model, the maximum activation is always
recorded for neuron 1849 for every prefix across all languages. Also across
all languages, the maximum standard deviation across prefixes was recorded for
neuron 218 for all languages. The single peak that we observe for the 2.9B
model in Figure 8 corresponds to the neuron 218. Thus, the aggregate value of
our activation flatness reflects the intuitive differences in representations.

For both detectors and combinators, we observe that the low values of
activation flatness occurs similarly across all languages. We therefore
investigate the extent to which the prefixes across languages exhibit
similarity in their layer-wise representations. Representational similarity
would point to the presence of multilingual neurons in those layers.

Input: Layer snapshot $L_{1}$ with shape $(m,d)$ and snapshot $L_{2}$ with
shape $(n,d)$

Output: Aggregate minimum distance across all rows of $L_{1}$, i.e. all
prefixes of the source sentence and their most similar counterpart prefixes in
the target language

$total\\_dist\leftarrow 0$;

for  _$i\leftarrow 1$ to $m$_ do

$min\\_dist\leftarrow+\infty$;

for  _$j\leftarrow 1$ to $n$_ do

$score\leftarrow dist({L_{1}[i],L_{2}[j]}$;

if  _$score <min\\_dist$_ then

$min\\_dist\leftarrow score$

end if

end for

$total\\_dist\leftarrow total\\_dist+min\\_dist$

end for

return _$total\\_dist$_

Algorithm 1 Representation Distance Score

We posit that if a layer has more multilingual neurons, then the
representation of two sentence prefixes with the same meaning should be
similar in that layer. In other words, the representation of “I am a computer”
should be very similar to “Ja jsem počítač” (in Czech). Thus for every prefix
in one language, we find the representational difference with all prefixes in
the second language. We then select the minimum distance among them. The idea
is to identify ‘similar’ prefixes across both languages. For simplicity, we
assume monotonic translation, which more or less holds for our studied
language pairs. We then aggregate the distances across all prefixes. The layer
with the minimal aggregate distance should exhibit the greatest
representational similarity for the two languages, as formally captured in
Algorithm 1.

Figure 9: Distance between representations: Detectors (568M, 1.7B. 2.9B, 7.5B)

Figure 10: Distance between representations: Combinators (568M, 1.7B. 2.9B,
7.5B)

From Figure 9 we see that the representational distance for detectors
decreases with layer depth for all models except the 7.5B parameter model.
Figure 10 shows that for the 2.9B parameter model, the representational
distance for combinators decreases in the middle, before increasing again.
Interestingly, the pattern of detector representational distance for 7.5B
seems similar to the pattern of combinator representational distance for the
2.9B model. We also observe that around layer 30 in the 7.5B parameter model,
the representational distance drops significantly for the combinators before
increasing abruptly through layers 31 and 32. We label the layers with the
minimal distance between prefixes of different language pairs as relatively
more _multilingual_.

Figure 11: Representational Similarity between models (combinators): 1.7B vs
2.9B, 2.9B vs 7.5B, 1.7B vs 7.5B

Finally, we also compare the “representational similarity” Kriegeskorte et al.
(2008) for German between the 1.7B, 2.9B and 7.5B parameter models in Figure
11. We calculate the representational similarity on the basis of the
normalized values of activation corresponding to each prefix for the
detectors/combinators of all layers. We see that the similarity is the maximum
for 2.9B and 7.5B models. Particularly, some of the middle layers of the 7.5B
parameter model shows high correlation with early and middle layers of the
2.9B parameter model.

The representational distance for all language pairs at the last layers of the
combinator is higher for the models with larger parameters. In other words, it
appears that the final layer of the larger models can make distinction between
different languages. We posit that this might be an explanation why the 7.5B
model outperforms models with smaller parameters of the same family as
reported by Lin et al. (2022).

Thus, based on the observations of activation flatness and representational
distances/similarity, we have a comprehensive view of the nature of
multilinguality in detectors and combinators. For all models except 7.9B, the
detectors keep getting more multilingual with layer depth. For the 7.5B model,
the middle layers are more multilingual than the layers near the input/output.
For the combinators, final layers are the most multilingual except the 2.9B
parameter model where the middle layers are the most multilingual.

###  4.3 The strange behaviour of the 2.9B model

We find that the combinators of the 2.9B billion parameter behave very
differently through the layers than the combinators of the other models. In
terms of the architectural micro-details, model dimension size and the hidden
dimension size of the 2.9B model is the same as that of the 1.7B model i.e.
2048 and 8192 respectively. The 2.9B model differs from the 1.7B model only in
terms of layer depth i.e. 24 vs 48. There is also representational similarity
across layers of the two models. We see that for both models, there is a
minimum in activation flatness (as well as minimum representational distance)
near layer 20. For the 2.9B billion parameter model, we also observe that the
representational distance gradually increases in the later layers. For the
final layer, we also find that for both the 7.5B and 2.9B models, the
representational distance among language pairs is almost comparable
($\sim$20000). One major point of difference between the 2.9B and 7.5B
parameter model (except the architectural micro-details) is the number of
layers in the model.

We suspect that this might make the 2.9B model bad at generation (in
comparison to the 7.5B model). We thus posit that the increase in number of
layers harms the generation quality of the 2.9B parameter model (over-
layerization). We suspect that longer training of the 2.9B model could have
fixed this but we have no further details about the training of XGLM.

The authors of the XGLM report that the 7.5B parameter model outperforms other
smaller models of the same family in mutilingual tasks including few-shot
translation. However, recent work by Zhang et al. (2023) involving QLoRA fine
tuning of LLMs for translation shows that after finetuning, the XGLM 2.9B
parameter model outperforms larger XGLM models in translation tasks. Is this
boost in performance somehow related to the activation flatness patterns
observed by us? We plan to investigate this aspect in future work,
hypothesizing that flatter activation patters could be easier to adapt in
QLoRA.

###  4.4 Activation similarity through languages

We have already seen some evidence that certain layers are more multilingual
than others based on our metric of representational distance/similarity. In
this section, we validate the observations and confirm if the same neurons
process different languages in the layers.

For the model snapshots of detectors and combinators from each layer, we take
the sum of activations across all prefixes and rank the neurons based on that.
The motivation here is that an aggregate picture might reveal which neurons
‘light up’ the most across prefixes.

We then calculate the Spearman rank correlation for a layer between language
pairs. Our hypothesis posits that high rank correlations between two languages
within a layer suggest similarities in their processing patterns. We present
the results below.

Figure 12: Rank correlation of neurons in detectors

Figure 13: Rank correlation of neurons in combinators

We observe that layers of detectors (Figure 12) exhibit decreasing rank
correlation across languages with increasing model size. However, the rank
correlation is consistently high for German-English-French language pairs
across models. From Figure 12 and Figure 13, we also find that the combinators
show greater rank correlation than detectors across languages and across
layers, aligning again with the observations made previously with
representational similarity.

For both detectors and combiantors, we observe high degrees of correlation in
the early and middle layers for all models. For the combinator we observe an
‘emerging’ pattern where the final layers initially exhibit greater values of
rank correlation (568M<1.7B) before decreasing again (1.7B>2.9B>7.5B). We also
observe that the middle combinator layers of the 7.5B model exhibited the
greatest rank correlation for the English-German-French pairs.

Thus, the early detectors are more multilingual for all models. We also
observe that the early and middle combinator layers become multilingual with
increasing model size.

Figure 14: Similarity of ranks (English) : detectors

Figure 15: Similarity of ranks (English) : combinators

#### Intra-language comparison

Next, we compare English sentences across different test sets. Figure 16 shows
the extent of overlap between the sentences across the 4 different English
test-sets used for the analysis. The hypothesis is that the analysis should
show if there were dedicated layers for language specific shallow processing
in the detectors or combinators.

Figure 16: Overlap of English sentences

We see from Figure 14 that detectors in the early and middle layers of all
models exhibit high degrees of rank correlation. We had already posited that
the early detectors are multilingual. We use these results to extend that
argument and state that this indicates the presence of “shallow” detectors
i.e. they do not deal with the ‘semantic’ content of the sentences. On the
other hand, the combinators (Figure 15) exhibit the maximal rank correlation
in the later layers (close to the output). We posit that this indicates the
presence of greater language specific neurons close to the output.

Figure 17: Similarity of ranks (parallel data) : detectors

Figure 18: Similarity of ranks (parallel data) : combinators

#### Parallel Data Comparison

Finally, we compare languages from parallel test-sets. The idea is to observe
if certain layers processed “semantic” content similarly for certain models.
We observe that for the English-German and English-French data, the detectors
exhibit high degree of correlation in the early layers reaffirming earlier
observations about their multilinugal and shallow processing abilities. For
the combinators we find high values of rank correlation across layers for
smaller models. The correlation is in fact highest when comparing German,
French and English sentences. The combinators of the 7.5B parameter model
however exhibits rank correlation close to zero across all layers for all
language pairs except German-English and French-English. Infact, for 7.5B, the
similarity is the most at layer 15, which also recorded the maximal
representational similarity in Figure 10. We posit that like detectors, the
combinators are multilinugal at early layers. The fact that the last layers of
the XGLM models are language specific is pretty interesting as it has already
been observed that for BERT, the upper layers (close to the ouptut) are more
language specific. We note that this is an interesting parallel between the
two kinds of LLMs (i.e. autoregressive versus auto-encoder).

##  5 Conclusion

In this work, we present our observations from analysing the FFN layers of
various models from the XGLM language family. We first identify that across
models, FFN representations are more sparse near the input and the output
layers. Adopting the detectors-combinators view of FFNs and looking deeper
into the nature of distribution of the activations, we find that for
detectors, the activations get peaked through the layers. For the combinators,
we observe that the activations get peaked across languages around layer 20,
irrespective of the model depth, model dimension size or the size of the
hidden representation in the FFNs. For combinators in fact, we discover that
the representations across languages are very similar until very deep layers.
Finally, we find that the detectors in the early layers process shallow
features and are multilingual. They become language specific close to the
output layers. We also find that the combinators are mostly multi-lingual
around the middle layers.

##  6 Limitations

While our study offers insights into the multilingual capabilities of the XGLM
model, several limitations must be acknowledged.

Firstly, our analysis relies only snapshots of the model’s activations, which
does not fully capture the complexities of LLMs. In addition to model
snapshots, probing tasks or even looking at the attention matrices could
provide a more comprehensive understanding of how the model processes
different languages.

Secondly, our investigation is limited to the XGLM model, and our findings may
not generalize to other LLM architectures or implementations.

Furthermore, our analysis is limited to a set of predefined languages and test
sets, which may not fully represent the diversity of languages and linguistic
phenomena captured by these models. Extending our analysis to more languages
and exploring cross-linguistic variation in more detail could yield further
insights into the generalizability of our findings.

## References

  * Belrose et al. (2023) Nora Belrose, Zach Furman, Logan Smith, Danny Halawi, Igor Ostrovsky, Lev McKinney, Stella Biderman, and Jacob Steinhardt. 2023.  Eliciting latent predictions from transformers with the tuned lens.  _arXiv preprint arXiv:2303.08112_.
  * Bhattacharya and Bojar (2023) Sunit Bhattacharya and Ondrej Bojar. 2023.  Unveiling multilinguality in transformer models: Exploring language specificity in feed-forward networks.  _EMNLP 2023_ , page 120.
  * Dar et al. (2023) Guy Dar, Mor Geva, Ankit Gupta, and Jonathan Berant. 2023.  Analyzing transformers in embedding space.  In _Proceedings of the 61st Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , pages 16124–16170.
  * Doddapaneni et al. (2021) Sumanth Doddapaneni, Gowtham Ramesh, Mitesh M Khapra, Anoop Kunchukuttan, and Pratyush Kumar. 2021.  A primer on pretrained multilingual language models.  _arXiv preprint arXiv:2107.00676_.
  * Geva et al. (2022) Mor Geva, Avi Caciularu, Kevin Wang, and Yoav Goldberg. 2022.  Transformer feed-forward layers build predictions by promoting concepts in the vocabulary space.  In _Proceedings of the 2022 Conference on Empirical Methods in Natural Language Processing_ , pages 30–45.
  * Geva et al. (2021) Mor Geva, Roei Schuster, Jonathan Berant, and Omer Levy. 2021.  Transformer feed-forward layers are key-value memories.  In _Proceedings of the 2021 Conference on Empirical Methods in Natural Language Processing_ , pages 5484–5495.
  * Greenberg (1966) Joseph H Greenberg. 1966.  _Language universals_.  M.I.T. Press.
  * Gurnee et al. (2024) Wes Gurnee, Theo Horsley, Zifan Carl Guo, Tara Rezaei Kheirkhah, Qinyi Sun, Will Hathaway, Neel Nanda, and Dimitris Bertsimas. 2024.  Universal neurons in gpt2 language models.  _arXiv preprint arXiv:2401.12181_.
  * Jastrzębski et al. (2017) Stanisław Jastrzębski, Devansh Arpit, Nicolas Ballas, Vikas Verma, Tong Che, and Yoshua Bengio. 2017.  Residual connections encourage iterative inference.  _arXiv preprint arXiv:1710.04773_.
  * Kriegeskorte et al. (2008) Nikolaus Kriegeskorte, Marieke Mur, and Peter A Bandettini. 2008.  Representational similarity analysis-connecting the branches of systems neuroscience.  _Frontiers in systems neuroscience_ , 2:249.
  * Lin et al. (2022) Xi Victoria Lin, Todor Mihaylov, Mikel Artetxe, Tianlu Wang, Shuohui Chen, Daniel Simig, Myle Ott, Naman Goyal, Shruti Bhosale, Jingfei Du, et al. 2022.  Few-shot learning with multilingual generative language models.  In _Proceedings of the 2022 Conference on Empirical Methods in Natural Language Processing_ , pages 9019–9052.
  * Liu et al. (2020) Chi-Liang Liu, Tsung-Yuan Hsu, Yung-Sung Chuang, and Hung-Yi Lee. 2020.  A study of cross-lingual ability and language-specific information in multilingual bert.  _arXiv preprint arXiv:2004.09205_.
  * Liu et al. (2023) Zichang Liu, Jue Wang, Tri Dao, Tianyi Zhou, Binhang Yuan, Zhao Song, Anshumali Shrivastava, Ce Zhang, Yuandong Tian, Christopher Re, et al. 2023.  Deja vu: Contextual sparsity for efficient llms at inference time.  In _International Conference on Machine Learning_ , pages 22137–22176. PMLR.
  * Mirzadeh et al. (2023) Iman Mirzadeh, Keivan Alizadeh, Sachin Mehta, Carlo C Del Mundo, Oncel Tuzel, Golnoosh Samei, Mohammad Rastegari, and Mehrdad Farajtabar. 2023.  Relu strikes back: Exploiting activation sparsity in large language models.  _arXiv preprint arXiv:2310.04564_.
  * Mueller et al. (2022) Aaron Mueller, Yu Xia, and Tal Linzen. 2022.  Causal analysis of syntactic agreement neurons in multilingual language models.  _arXiv preprint arXiv:2210.14328_.
  * Muller et al. (2021) Benjamin Muller, Yanai Elazar, Benoît Sagot, and Djamé Seddah. 2021.  First align, then predict: Understanding the cross-lingual ability of multilingual bert.  _arXiv preprint arXiv:2101.11109_.
  * Pires et al. (2019) Telmo Pires, Eva Schlinger, and Dan Garrette. 2019.  How multilingual is multilingual bert?  _arXiv preprint arXiv:1906.01502_.
  * Rogers et al. (2021) Anna Rogers, Olga Kovaleva, and Anna Rumshisky. 2021.  A primer in bertology: What we know about how bert works.  _Transactions of the Association for Computational Linguistics_ , 8:842–866.
  * Song et al. (2024) Chenyang Song, Xu Han, Zhengyan Zhang, Shengding Hu, Xiyu Shi, Kuai Li, Chen Chen, Zhiyuan Liu, Guangli Li, Tao Yang, et al. 2024.  Prosparse: Introducing and enhancing intrinsic activation sparsity within large language models.  _arXiv preprint arXiv:2402.13516_.
  * Stańczak et al. (2022) Karolina Stańczak, Edoardo Ponti, Lucas Torroba Hennigen, Ryan Cotterell, and Isabelle Augenstein. 2022.  Same neurons, different languages: Probing morphosyntax in multilingual pre-trained models.  _arXiv preprint arXiv:2205.02023_.
  * Tamkin et al. (2023) Alex Tamkin, Mohammad Taufeeque, and Noah D Goodman. 2023.  Codebook features: Sparse and discrete interpretability for neural networks.  _arXiv preprint arXiv:2310.17230_.
  * Xie et al. (2021) Wanying Xie, Yang Feng, Shuhao Gu, and Dong Yu. 2021.  Importance-based neuron allocation for multilingual neural machine translation.  In _Proceedings of the 59th Annual Meeting of the Association for Computational Linguistics and the 11th International Joint Conference on Natural Language Processing (Volume 1: Long Papers)_ , pages 5725–5737.
  * Zhang et al. (2023) Xuan Zhang, Navid Rajabi, Kevin Duh, and Philipp Koehn. 2023.  Machine translation with large language models: Prompting, few-shot learning, and fine-tuning with qlora.  In _Proceedings of the Eighth Conference on Machine Translation_ , pages 468–481.

## Appendix