# Quenched disorder and instability control dynamic fracture in three
dimensions

Yuri Lubomirsky  Eran Bouchbinder  Chemical and Biological Physics Department,
Weizmann Institute of Science, Rehovot 7610001, Israel

###### Abstract

Materials failure in 3D still poses basic challenges. We study 3D brittle
crack dynamics using a phase-field approach, where Gaussian quenched disorder
in the fracture energy is incorporated. Disorder is characterized by a
correlation length $R$ and strength $\sigma$. We find that the mean crack
velocity $v$ is bounded by a limiting velocity, which is smaller than the
homogeneous material’s prediction and decreases with $\sigma$. It emerges from
a dynamic renormalization of the fracture energy with increasing crack driving
force $G$, resembling a critical point, due to an interplay between a 2D
branching instability and disorder. At small $G$, the probability of localized
branching on a scale $R$ is super-exponentially small. With increasing $G$
this probability quickly increases, leading to misty fracture surfaces, yet
the associated extra dissipation remains small. As $G$ is further increased,
branching-related lengthscales become dynamic and persistently increase,
leading to hackle-like structures and to a macroscopic contribution to the
fracture surface. The latter dynamically renormalizes the actual fracture
energy until eventually any increase in $G$ is balanced by extra fracture
surface, with no accompanying increase in $v$. Finally, branching width
reaches the system’s thickness such that 2D symmetry is statistically
restored. Our findings are consistent with a broad range of experimental
observations.

Materials failure is a physical phenomenon of prime scientific and
technological importance, which continues to pose fundamental challenges [1,
2, 3, 4, 5]. It is mediated by the propagation of cracks, which are dynamic
defects that feature extreme conditions near their edges — mimicking a
mathematical singularity — and accelerate to relativistic velocities (relative
to elastic wave-speeds) in brittle materials [6, 7]. Moreover, crack dynamics
involve strongly non-equilibrium physics, long-range elastodynamic
interactions and a wide range of spatiotemporal scales, which still resist a
complete theoretical understanding.

Things somewhat simplify in 2D, i.e., when experiments in thin samples are
considered [8, 9], where a crack is the trajectory left behind a propagating
tip (a point). Major progress in understanding rapid crack dynamics in 2D has
been recently made [9, 10, 11, 12, 13]. In particular, high-velocity symmetry-
breaking instabilities have been predicted, in quantitative agreement with
experiments, including the 2D oscillatory instability [8] and the 2D branching
instability (sometimes also termed macro-branching/tip-splitting) [9]. These
two instabilities are linear instabilities, which are spontaneously triggered
by infinitesimal noise, once the critical conditions are met.

Our corresponding understanding of dynamic fracture in 3D lags far behind. In
3D, the sample thickness is no longer a negligible dimension and a crack is
the surface left behind a propagating line, i.e., here the crack edge is a 2D
front, see Fig. 1a. While available experiments on brittle amorphous materials
shed important light on the phenomenology of dynamic fracture in 3D [14, 15,
16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28], our fundamental
understanding of the underlying physics remains limited. These experiments
revealed a variety of 3D out-of-plane crack structures, including micro-
cracking damage [14, 15, 16, 17], surface steps [18, 19], cross-hatching
patterns [20], micro-branches [21, 22, 23, 24] and mirror-mist-hackle patterns
[25, 26, 27, 28], all emerging at mean crack front velocities smaller than the
2D instability thresholds.

Figure 1: (a) A sketch of a 3D tensile crack experiment in a long bar of
height $L_{y}$ and thickness $L_{z}$. The long bar feature a length
$L_{x}\\!\gg\\!L_{y},L_{z}$, and only a small section of it is shown. The
crack is located in the middle (symmetry) plane and will subsequently
propagate predominantly in the $x$ direction in response to sufficiently large
crack driving force $G$ (related to the tensile loading, applied in the $y$
direction). We employ periodic boundary conditions along the $z$ direction
[29]. The coordinates are normalized by the dissipation length $\xi$, see text
for definition. (b) A linear crack branching instability in 2D, i.e.,
$L_{z}\\!=\\!0$, occurring at a very high propagation velocity $v_{{}_{\rm
B}}\\!=\\!0.95c_{{}_{\rm R}}$ ($c_{{}_{\rm R}}$ is the Rayleigh wave-speed),
once a threshold tensile driving force $G_{\rm B}$ is surpassed. See text and
Fig. 2 for additional discussion.

Here, we study 3D dynamic fracture using a flexible computational framework in
which material quenched disorder is incorporated. It is based on a phase-field
fracture approach [9, 12, 13, 30, 31, 32, 33], which is particularly suitable
for studying 3D dynamic fracture as it allows cracks to self-consistently
select complex 3D trajectories (including topological changes) without
imposing any external path selection criteria. Moreover, it allows to track
the in silico real-time 3D spatiotemporal dynamics of cracks in a way that
goes well beyond current experiments. This framework quantitatively predicted
the 2D high-velocity oscillatory and branching instabilities in homogeneous
materials [9, 12, 13], as well as the dynamics of crack front waves in 3D in
homogeneous materials with isolated heterogeneities [33].

Evidence regarding 3D dynamic fracture is largely obtained in experiments on
amorphous materials, such as various glassy polymers [14, 15, 16, 17, 21, 22,
23, 28], silica glasses [22, 25, 26, 27] and brittle elastomers [24]. These
materials are intrinsically disordered, featuring fluctuations in material
properties over various lengthscales. Material disorder and heterogeneity have
been extensively discussed mainly in the context of in-plane fracture
roughness in 3D and mostly in the quasi-static regime (e.g., [3, 34]). Some
recent studies considered the effect of material heterogeneity on dynamic in-
plane cracks [35] and on quasi-static out-of-plane cracks [36]. Yet, the
interaction of a 3D elastodynamic crack with continuous quenched disorder has
not been studied theoretically. We show that the interplay between quenched
disorder and the 2D branching instability controls dynamic fracture in 3D, in
agreement with a wide range of experimental observations.

Dimensionality, stability and disorder-induced limiting velocity

A canonical experiment for probing failure dynamics involves an initial planar
crack located in the middle plane of a long 3D rectangular bar, see Fig. 1a.
The bar features height $L_{y}$, thickness $L_{z}$ and length
$L_{x}\\!\gg\\!L_{y},L_{z}$, where the $(x,y,z)$ coordinate system is defined
in Fig. 1. It is loaded symmetrically by small tensile displacements
$u_{y}(x,y\\!=\\!\pm L_{y}/2,z,t)\\!=\\!\delta/2$, where $\bm{u}(x,y,z,t)$ is
the 3D displacement field and $t$ is time. The tensile loading corresponds to
a crack driving force $G\\!\sim\\!E\delta^{2}/L_{y}$ (the linear elastic
energy per unit area stored far ahead of the initial crack) [6, 7], where $E$
is Young’s modulus. Once $G$ surpasses the minimal quasi-static fracture
energy $\Gamma_{0}$, the tensile (mode-I) crack front starts propagating
predominantly in the $x$ direction. Understanding the subsequent
spatiotemporal dynamics as a function of $G\\!>\\!\Gamma_{0}$ remains a major
open challenge.

Experiments are limited in probing the real-time spatiotemporal dynamics of
cracks in 3D, mainly because continuous imaging of extended dynamic defects
that propagate at high velocities inside a material is currently not
available. Consequently, large-scale computer simulations offer a powerful
complementary approach. The phase-field approach, mentioned above, allows to
numerically solve the 3D problem delineated in Fig. 1a, while accurately
probing the in silico real-time 3D spatiotemporal dynamics of the crack. This
is achieved by solving the linear elastodynamic field equations for
$\bm{u}(x,y,z,t)$ inside the bar, coupled to an auxiliary field — the scalar
phase-field $\phi(x,y,z,t)$ [9, 12, 13, 33, 29]. The latter satisfies its own
dissipative field equation, which features a characteristic dissipation length
$\xi$ and a dissipation time $\tau$ [9, 12, 13, 33, 29].

In the presence of the intense, nearly singular deformation fields in the
vicinity of the moving crack front, $\phi(x,y,z,t)$ both spontaneously
generates the traction-free surfaces that define the crack [9, 12, 13, 33, 29]
— and hence self-consistently selects the crack trajectory — and gives rise to
a rate-dependent fracture energy $\Gamma(v)$ controlled by the dissipation
time $\tau$, where $v$ is the crack velocity (note that
$\Gamma(v\\!\to\\!0)\\!=\\!\Gamma_{0}$). Another major advantage of this
computational approach, essential for the present work, is the ability to tune
material properties in a controlled manner in ways that are difficult to
obtain in the laboratory. The phase-field fracture approach played a central
role in recent progress in understanding rapid crack dynamics in 2D [9, 12,
13]. The latter corresponds to taking the $L_{z}/\xi\\!\to\\!0$ limit in the
3D setup of Fig. 1a, leading to a long strip configuration [6, 7].

In Fig. 2, we plot (green circles) the imposed crack driving force
$G/\Gamma_{0}$ against the resulting steady-state crack velocity $v$
(normalized by the shear wave-speed $c_{\rm s}$) obtained in 2D phase-field
calculations in a homogeneous material [13]. It is observed that $G$ is a
monotonically increasing function of $v$, where each point along the curve
corresponds to a straight crack that respects the global tensile (mode-I)
symmetry of the system. Moreover, the curve is linear up to a velocity close
to the Rayleigh wave-speed, here $c_{{}_{\rm R}}\\!=\\!0.93c_{\rm s}$, which
is predicted to be the upper bound on the velocity of cracks based on Linear
Elastic Fracture Mechanics (LEFM) theory, in the absence of symmetry-breaking
instabilities [6, 7]. Beyond this regime, $G$ strongly increases with $v$ due
to crack tip blunting, leading to increased dissipation [9, 13].

Figure 2: The crack driving force $G/\Gamma_{0}$ vs. steady-state crack
velocity $v/c_{\rm s}$ for a 2D homogeneous material (green circle) and two 3D
heterogeneous materials with $\sigma\\!>\\!0$ (yellow diamonds for
$\sigma\\!=\\!0.25$ and brown diamonds for $\sigma\\!=\\!0.50$), obtained in
phase-field simulations [29]. The thick-boundary diamonds correspond to the 5
snapshots shown in Fig. 3. Steady-state solutions in 2D homogeneous materials
do not exist for $G\\!>\\!G_{\rm B}$ (horizontal dashed-dotted line), upon
which a linear branching instability sets in (see Fig. 1b). In 3D materials,
the average steady-state crack velocity is bounded by a $\sigma$-dependent
limiting velocity $v_{{}_{\rm lim}}(\sigma)$ (vertical dashed lines), where
$G$ increases in a critical-like manner. (inset) $G$ (normalized by the first
data point $G_{1}$) vs. steady-state crack velocity $v/c_{\rm R}$ measured in
experiments on a glassy polymer (Polymethyl Methacrylate, PMMA), data
extracted from Fig. 5 of [21]. The vertical dashed line, highlighting the
similarity to the simulational data in the main panel, is added as a guide to
the eye.

The 2D $G$ vs. $v$ curve terminates at $G_{\rm B}$ (horizontal dashed-dotted
line), at a velocity very close to — yet slightly smaller than — $c_{{}_{\rm
R}}$, upon which a branching instability (sometimes termed macro-
branching/tip-splitting) takes place. That is, a straight tensile crack
branches/splits into two coexisting cracks [9, 13], shown in Fig. 1b,
providing a mechanism for excess dissipation due to additional fracture
surfaces [37, 38]. For $G\\!>\\!G_{\rm B}$, steady-state solutions no longer
exist in 2D. The high-velocity 2D branching instability, previously
demonstrated in [9, 13], has been experimentally observed in quasi-2D brittle
materials [15, 9, 28, 23, 39, 40, 41]. Finally, since in steady-state $G$ is
balanced by the fracture energy, we identify the $G$ vs. $v$ curve as
$\Gamma_{\\!{}_{\rm 2D}}(v)$ and hence the 2D branching velocity $v_{{}_{\rm
B}}$ corresponds to $\Gamma_{\\!{}_{\rm 2D}}(v_{{}_{\rm B}})\\!=\\!G_{\rm B}$.
We note in passing that the 2D oscillatory instability, which emerges at a
velocity slightly smaller than $v_{{}_{\rm B}}$ [12], does not emerge here
[29].

Branching in 2D exists in purely homogeneous materials and hence is a linear
instability. We next considered the effect of dimensionality on the stability
of cracks in homogeneous materials by studying 3D systems with a finite
thickness $L_{z}$, as in Fig. 1a, driven by $G\\!<\\!G_{\rm B}$. We found that
tensile cracks in our 3D homogeneous system are remarkably stable, both
linearly and nonlinearly. In particular, for a wide variety of crack shape
perturbations (of different spatial extent and amplitudes) as well as for a
wide range of isolated heterogeneities (of different sizes, shapes and
locations), cracks transiently go out of their initial plane, but then decay
back to the tensile symmetry plane. That is, no persistent 3D out-of-plane
crack instabilities exist in homogeneous materials for $G\\!<\\!G_{\rm B}$,
hence for $v\\!<\\!v_{{}_{\rm B}}$ (persistent in-plane crack front waves are
triggered, as discussed very recently in [33]). This is in sharp and
qualitative disagreement with a wide range of experimental observations that
show that tensile cracks in 3D brittle amorphous materials reveal persistent
out-of-plane structures at velocities smaller than $v_{{}_{\rm B}}$, as
already noted above.

What is the origin of this understanding gap? What is the crucial physical
ingredient missing? In order to gain insight into these fundamental questions,
we stress that the experimental evidence at hand was obtained for amorphous
materials, which are homogeneous at macroscopic scales, but feature disorder
at smaller scales. The quantitative characterization and understanding of
disordered materials is still incomplete; in particular, while it is known
that various material properties (e.g., elastic moduli [42, 43]) feature
significant fluctuations over supermolecular lengthscales, we do not yet have
robust and quantitative knowledge of the statistics of fundamental physical
quantities such as the fracture energy. Consequently, in order to test the
possible importance of material quenched disorder, our strategy would be to
adopt minimal assumptions about material quenched disorder and explore its
implications for bridging the qualitative gap in our understanding of 3D
dynamic fracture.

Specifically, we introduce quenched disorder in the form of a Gaussian
distribution in the dimensionless quasi-static ($v\\!\to\\!0$) fracture energy
$\bar{\Gamma}/\Gamma_{0}$, with unity mean, standard deviation determined by
$\sigma$ and spatial correlation length $R$ [29]. We hereafter set
$R\\!=\\!10\xi$ and study the effect of $\sigma\\!>\\!0$ on 3D dynamic
fracture. In Fig. 2, we plot (diamonds) $G/\Gamma_{0}$ against the average
steady-state crack velocity $v/c_{\rm s}$ obtained in 3D phase-field
calculations as in Fig. 1a, for two values of $\sigma$. Recall that in steady
state, energy balance implies $\Gamma(v)\\!=\\!G$ [6, 7], where $\Gamma(v)$ is
the actual fracture energy. It is observed that $\Gamma(v)\\!=\\!G$ closely
follows the 2D homogeneous material curve at relatively small $v$, then it
deviates from it at larger $v$’s until it reaches a limiting (terminal)
velocity $v_{{}_{\rm lim}}(\sigma)\\!<\\!v_{{}_{\rm B}}\\!<\\!c_{{}_{\rm R}}$,
consistently with numerous experiments (e.g., [2, 21] and the compilation of
experimental works in Table. 1 of [15]). Moreover, $d\Gamma/dv$ appears to
diverge for $v\\!\to\\!v_{{}_{\rm lim}}$, also consistently with experiments
[21], cf. inset of Fig. 2.

The results in Fig. 2 can be described by the steady-state energy balance

$\Gamma(v)=\Gamma_{\\!{}_{\rm 2D}}(v)+\delta\Gamma_{\\!{}_{\rm
3D}}(v,\sigma)=G\ ,$ (1)

where $\delta\Gamma_{\\!{}_{\rm 3D}}(v,\sigma)$ corresponds to 3D disorder-
induced corrections to the 2D homogeneous-material $\Gamma_{\\!{}_{\rm
2D}}(v)$. We note that we do not explicitly incorporate the dependence of
$\delta\Gamma_{\\!{}_{\rm 3D}}(v,\sigma)$ on the correlation length $R$ since
it is kept constant, as stated above. For relatively small velocities, we have
$\delta\Gamma_{\\!{}_{\rm 3D}}(v,\sigma)\\!\ll\\!\Gamma_{\\!{}_{\rm
2D}}(v)\\!\simeq\\!\Gamma(v)$ and the behavior essentially identifies with
that of 2D homogeneous materials. At some $\sigma$-dependent characteristic
$v$, $\delta\Gamma_{\\!{}_{\rm 3D}}(v,\sigma)$ is no longer negligible
compared to $\Gamma_{\\!{}_{\rm 2D}}(v)$, i.e., 3D disorder-induced effects
become important. Finally, as $v\\!\to\\!v_{{}_{\rm lim}}$,
$\delta\Gamma_{\\!{}_{\rm 3D}}(v,\sigma)$ becomes sizable/dominant and appears
to feature a critical behavior, i.e., $d\Gamma(v)/dv$ appears to diverge. Our
main challenge in the remaining of the paper is to understand the physics
behind $\delta\Gamma_{\\!{}_{\rm 3D}}(v,\sigma)$.

Figure 3: (left column, panels a-e) Uppermost out-of-plane crack structures —
corresponding to the phase-field $\phi({\bm{x}})\\!=\\!1/2$ iso-surface [29] —
inside the bulk from a top view, for the five $G$ values indicated
(corresponding to the thick-boundary diamonds in Fig. 2), see text for
extensive discussion. In each panel, the long axis is the propagation ($x$)
direction and the perpendicular one in the crack front direction ($z$ axis).
Both axes are expressed in units of $\xi$, and the basic length along $z$ is
twice that of along $x$, for visual clarity. The initial crack front is marked
by the bright line in the left part of each panel. (middle column, panels f-g)
Tilted zoom-in 3D view of the regions marked by the rectangles in panels (b)
and (e), respectively. See text for discussion. Note that in panel (g), the
crack symmetry plane corresponds to $y\\!=\\!0$ in order to highlight the out-
of-plane structures emerging both above and below the symmetry plane. (right
column, panels i-h) Fracture surface morphologies (fractography) measured in
dynamic fracture experiments in a glassy polymer (PMMA) involving accelerating
3D cracks (propagating from left to right), obtained through SEM imaging [28].
Panel (h) shows the mirror-to-mist transition, where the mist region is
characterized by rather isolated, small-amplitude out-of-plane surface
structures (adapted from Fig. 4b of [28]). Initially, the crack does not
capture the entire sample thickness (vertical axis in the image) and it
propagates predominantly radially. Panel (h) bears qualitative similarities to
panels (b) and (f). Panel (i) shows the hackle region, appearing at larger
crack propagation distances in the same experiment (adapted from Fig. 4d of
[28]). The observed wider and elongated out-of-place surface structures bear
clear similarities to panels (c) and (d).

Out-of-plane crack structures and the localized branching instability

To start addressing this challenge, we focus on out-of-plane crack structures
that correspond to the 3D results in Fig. 2 with $\sigma\\!=\\!0.5$ (to be
considered hereafter). Since the crack front experiences both in- and out-of-
plane fluctuations, the quantity $v$ hereafter stands for the averaged front
velocity [29]. We take advantage of our computational platform and consider
out-of-plane crack structures inside the bulk of the material, not just their
fractographic signature (i.e., the patterns left on the postmortem fracture
surface). In Fig. 3 (left column), we present a sequence of out-of-plane crack
structures as a function of the crack driving force $G$, covering the various
regimes of interest (the data points used are marked by thick boundaries in
Fig. 2).

Figure 3a corresponds to $G/\Gamma_{0}\\!=\\!2$, for which
$\delta\Gamma_{\\!{}_{\rm 3D}}(v,\sigma)\\!\ll\\!\Gamma_{\\!{}_{\rm 2D}}(v)$.
It is observed that out-of-plane protrusions are very rare and small, i.e.,
the crack is essentially planar (mirror-like), which is consistent with the
agreement between $\Gamma(v)$ and $\Gamma_{\\!{}_{\rm 2D}}(v)$ in this regime.
Figure 3b corresponds to $G/\Gamma_{0}\\!=\\!2.5$, where the deviation of
$\Gamma(v)$ from $\Gamma_{\\!{}_{\rm 2D}}(v)$ is still small. It is observed
that out-of-plane structures, of width $\Delta{z}$ and length $\Delta{x}$,
emerge at a significantly higher probability compared to panel (a). The out-
of-plane structures, leading to a misty fracture surface, are mainly localized
branching events of characteristic width comparable to the disorder
correlation length $R\\!=\\!10\xi$, as highlighted in the tilted zoom-in view
in Fig. 3f. In addition, each localized branching event is accompanied by
v-shaped tracks (clearly observed in Fig. 3f), which are reminiscent of
typical tracks left by crack front waves [44, 45]. Indeed, the propagation
velocity of these tracks is consistent with that of crack front waves [29].

The character of the out-of-plane structures appearing in Fig. 3c,
corresponding to $G/\Gamma_{0}\\!=\\!3.75$, qualitatively changes. The
localized branching events therein become significantly wider and longer,
i.e., leading to hackle-like structures featuring larger $\Delta{z}$ and
$\Delta{x}$, hence larger areas. With increasing $G$, as in Figs. 3d-e,
$\Delta{z}$ and $\Delta{x}$ further increase. In Fig. 3e, the width
$\Delta{z}$ becomes (statistically) comparable to the system thickness
$L_{z}$, which implies that 2D symmetry is statistically restored. This
3D-to-2D transition is manifested in branching events that are no longer
localized, i.e., can be regarded as macro-branches [15, 23, 39] since they
capture the entire system thickness and propagate sizable distances, as
highlighted in the tilted zoom-in view in Fig. 3g. Figures 3c-e correspond to
$G$ values for which the crack reached the limiting velocity,
$v\\!\to\\!v_{{}_{\rm lim}}$, where $\delta\Gamma_{\\!{}_{\rm 3D}}(v,\sigma)$
makes a sizable/dominant contribution to $\Gamma(v)$ and $d\Gamma(v)/dv$
appears to diverge.

Figure 4:  (a) The number of branches per unit planar area as a function of
$G/G_{\rm B}\\!\leq\\!1$ (squares). The solid line corresponds to the
prediction in Eq. (2), where $\sigma_{\\!{}_{\rm R}}\\!=\\!0.233$ is
independently predicted [29] and the overall amplitude is set to agree with
the first data point. (b) The normalized weighted branch length
$\langle\Delta{x}\rangle_{\rm w}/R$ (see text for definition) vs. $G/G_{\rm
B}$. The vertical dashed line marks $G\\!=\\!G_{\rm B}$. (c) The same as panel
(a), but vs. $v$. The vertical dashed line corresponds to $v_{{}_{\rm lim}}$,
as in Fig. 2.

The probabilistic crossover between Fig. 3a and Fig. 3b closely resembles the
widely-observed mirror-to-mist transition in brittle amorphous materials [25,
26, 27, 28], see Fig. 3h. The transition in out-of-plane crack structures
observed between Fig. 3b and Fig. 3c bears strong similarity to the widely
documented mist-to-hackle transition, se Fig. 3i, observed in many brittle
amorphous materials, including glassy polymers (e.g., [28]) and silica glasses
(e.g. [27]). Moreover, the 3D-to-2D transition, which is accompanied by
macroscopic crack branching and that follows the hackle regime with increasing
crack driving force $G$, is also well documented in the very same materials
(e.g., [15, 23, 27, 28]). The correspondence between our findings and a wide
range of experimental observations will be further discussed once we gain
deeper understanding of the former.

A physical picture and supporting numerics

Our first goal is to understand the emergence of the localized branching
instability in 3D, along with its statistical and dynamical properties. The 2D
branching instability in homogeneous materials, cf. Figs. 1a and 2, is
triggered at a critical driving force $G_{\rm B}/\Gamma_{0}$ (see horizontal
dashed-dotted line in Fig. 2). What is the relevance of this instability in 3D
in the presence of disorder?

On the scale of the correlation scale $R$, the disordered material can be
viewed as locally 2D. Assuming that $G$ in this regime of the 3D dynamics is
rather homogeneously distributed along the thickness ($z$ axis), we expect
every fluctuation $\bar{\Gamma}$ in the quasi-static fracture energy that
satisfies $G/\bar{\Gamma}\\!>\\!G_{\rm B}/\Gamma_{0}$ to give rise to a
localized branching event on a scale $R$. Rearranging the inequality to read
$\bar{\Gamma}/\Gamma_{0}\\!<\\!G/G_{\rm B}$ and assuming localized branching
events to be independent of one another in this regime such that the Gaussian
(normal) distribution of $\bar{\Gamma}/\Gamma_{0}$ (with unity mean and
standard deviation $\sigma$) can be invoked, we obtain [29]

$p(G;\sigma)\sim 1+\text{erf}\left(\frac{G/G_{\rm B}-1}{\sigma_{\\!{}_{\rm
R}}}\right)\ .$ (2)

Here $p(G;\sigma)$ is the localized branching probability, $\text{erf}(\cdot)$
is the error function, which varies super-exponentially with its argument, and
$\sigma_{\\!{}_{\rm R}}\\!\sim\\!\sigma$ is the renormalized strength of
disorder on a scale $R$ [29].

Equation (2) predicts an extremely small localized branching probability at
small $G$ values and a strong increase of the probability with increasing $G$
below $G_{\rm B}$, in qualitative agreement with the results Fig. 3a-b. To
more directly test the prediction in Eq. (2), we plot in Fig. 4a the number of
localized branches per unit planar area, which is proportional to
$p(G;\sigma)$, as a function of $G\\!<\\!G_{\rm B}$. It is observed that the
simple model agrees with the data reasonably well when $G$ is not too close to
$G_{\rm B}$.

While 2D branching is a linear instability triggered in purely homogeneous
materials (that can be realized in silico) by infinitesimal noise, the
corresponding 3D localized branching events with
$\bar{\Gamma}/\Gamma_{0}\\!<\\!G/G_{\rm B}$ are finite-disorder instabilities,
and hence are non-perturbative with respect to a homogenized material.
Moreover, 3D localized branching events with
$\bar{\Gamma}/\Gamma_{0}\\!<\\!G/G_{\rm B}\\!<\\!1$ are expected to be rather
compact, roughly of linear scale $R$, and hence short-lived and contributing
small excess fracture surface. The reason for this is that while localized
branching can be triggered with $G/G_{\rm B}\\!<\\!1$, it cannot be sustained
over scales significantly larger than $R$ (since large local fracture energy
fluctuations are likely to be encountered). On the other hand, we expect a
change in the dynamics of 3D localized branching events for $G\\!>\\!G_{\rm
B}$, since under these conditions the branching solution has a sizable
probability to be sustained over larger distances — i.e., to feature
significantly larger length $\Delta{x}$ — while co-existing with the planar
crack solution.

To test this expectation, we plot in Fig. 4b $\langle\Delta{x}\rangle_{\rm
w}/R$, the average of $\Delta{x}$ weighted by the branch area (normalized by
$R$), against $G/G_{\rm B}$. It is observed that $\langle\Delta{x}\rangle_{\rm
w}$ indeed significantly grows with $G$ above $G/G_{\rm B}\\!\simeq\\!1$
(marked by the vertical dashed line), reaching values much larger than $R$, as
expected and in qualitative agreement with Figs. 3c-e. We note that we
consider the averaged branch length weighted by its area because while at
smaller $G$ values branching events are independent of one another (cf. Fig.
3b), for $G\\!>\\!G_{\rm B}$ branches effectively “screen” the nucleation of
other branches and in fact the number of primary branches decreases with
increasing $G$ (though secondary structures emerge as well [29]). Moreover,
for $G\\!>\\!G_{\rm B}$ the branch width $\Delta{z}$ is no longer expected to
inherit its scale from $R$. Indeed, the results presented in the inset of Fig.
5 indicate a quasi-linear relation between the averaged branch width
$\langle\Delta{z}\rangle$ and length $\langle\Delta{x}\rangle$. Fluctuations,
to be discussed below, also grow significantly with $G$, eventually
corresponding to the statistical 3D-to-2D transition observed in Fig. 3g.
Overall, the results indicate that for $G\\!>\\!G_{\rm B}$ branching
lengthscales become dynamic and a significant extra fracture surface emerges.

The emergence of a limiting velocity and dynamic renormalization of the
fracture energy

We are now in a position to rationalize the steady-state dynamics of 3D cracks
as encapsulated in Eq. (1) and in particular the physics underlying
$\delta\Gamma_{\\!{}_{\rm 3D}}(v,\sigma)$. The physical picture discussed
above and its numerical support indicate that 3D out-of-plane fracture
dynamics are related to a branching instability controlled by the driving
force $G$ and to quenched disorder, where both the branching probability
$p(G;\sigma)$ and the excess fracture surface associated with branching play a
role. In particular, at small $G$ the evolution of $p(G;\sigma)$ is dominant,
while for $G\\!>\\!G_{\rm B}$ the extra surface of the branches dominates,
with a transition in between.

To relate this picture to the emergence of a limiting velocity $v_{{}_{\rm
lim}}$ and to the properties of $\delta\Gamma_{\\!{}_{\rm 3D}}(v,\sigma)$, we
replot in Fig. 4c the number of localized branches per unit planar area of
Fig. 4a as a function of $v$, i.e., a quantity proportional to $p(v;\sigma)$
instead of $p(G;\sigma)$. It is observed that $p(v;\sigma)$ appears to feature
a diverging derivative as $v_{{}_{\rm lim}}$ is approached (vertical dashed
line). That is, the branching probability $p(v;\sigma)$ alone reveals a clear
signature of the emergence of a limiting velocity as $G$ approaches $G_{\rm
B}$. For $G\\!>\\!G_{\rm B}$, the branching probability does not further grow
(as noted above, it actually decreases), but rather the area of the branches
increases such that any increase in $G$ is expected to be balanced by excess
fracture surface, without increasing $v$, leading to an increase in
$\delta\Gamma_{\\!{}_{\rm 3D}}(v,\sigma)$ at a fixed $v\\!=\\!v_{{}_{\rm
lim}}$.

To test this expectation, we define the apparent areal ratio as the total area
of the uppermost out-of-plane crack structures inside the bulk (as seen in the
top view in the left column of Fig. 3) divided by the nominal (planar) area.
The true areal ratio may in fact be larger as secondary structures can develop
beneath the uppermost out-of-plane crack structures [29]. The apparent areal
ratio is plotted (diamonds, left $y$-axis) in Fig. 5 vs. $v$ and is observed
to be indistinguishable from unity until $v_{{}_{\rm lim}}$ is approached. At
$v\\!=\\!v_{{}_{\rm lim}}$, the apparent areal ratio increases in a critical-
like manner, as expected. It thus supports the idea that the 3D disorder-
induced excess dissipation $\delta\Gamma_{\\!{}_{\rm 3D}}(v,\sigma)$ balances
any increase in $G$ at $v\\!=\\!v_{{}_{\rm lim}}$ by creating extra fracture
surface. That is, the fracture energy $\Gamma(v)$ is dynamically renormalized
by the out-of-plane spatiotemporal dynamics of the 3D crack.

Figure 5: The apparent area ratio (see text for definition) vs. $v/c_{\rm s}$
(green diamonds, left $y$-axis). Out-of-plane fluctuations (see text for
definition) in units of $\xi$ vs. $v/c_{\rm s}$ (brown circles, right
$y$-axis, and note the logarithmic scale). The vertical dashed line marks
$v_{{}_{\rm lim}}$. (inset) The averaged branch length
$\langle\Delta{x}\rangle/R$ vs. width $\langle\Delta{z}\rangle/R$.

Finally, we expect fluctuations in out-of-plane crack structures to
significantly increase as $v\\!\to\\!v_{{}_{\rm lim}}$. To test this
expectation, we extracted for each $(x,y\\!=\\!0,z)$ point along the original
symmetry plane the $y$ difference (in units of $\xi$) between the uppermost
out-of-plane crack location (shown on the left column of Fig. 3) and its
lowermost counterpart. We then averaged over all $(x,y\\!=\\!0,z)$ points and
plotted the results in Fig. 5 (circles, right $y$-axis) vs. $v$ on a log-
linear scale. It is observed that the out-of-plane fluctuations are very small
for small $v$ (note the logarithmic right $y$-axis), consistently with the
branching probability $p(v;\sigma)$ being small in this regime. Around
$v/c_{\rm s}\\!=\\!0.7$, the fluctuations start to grow, yet they still make a
negligible contribution to the extra fracture area (see diamonds and left
$y$-axis). As $v\\!\to\\!v_{{}_{\rm lim}}$, out-of-plane fluctuations grow
significantly, featuring the same critical-like behavior as the apparent areal
ratio, as expected from the emerging physical picture.

Discussion

Our results indicate that 3D dynamic fracture in brittle amorphous materials
is controlled by two physical ingredients. The first is a 2D linear
instability that is characterized by a crack driving force threshold $G_{\rm
B}$, providing a mechanism for generating extra fracture surface, i.e., for
increasing the effective fracture energy. The 2D homogeneous material
instability involves a topological change, going beyond a single crack, in the
form of crack branching. The second ingredient is finite-strength quenched
disorder that features a spatial correlation length $R$.

In the presence of finite disorder, the 2D branching instability is
transiently excited in 3D at crack driving forces $G$ well below $G_{\rm B}$.
It features localization over a scale $R$ that breaks translational invariance
along the crack front. Upon increasing $G$ in 3D, the localized branching
instability becomes longer-lived and features larger lengthscales for
$G\\!\simeq\\!G_{\rm B}$. The extra fracture surface associated with these
increasingly larger branches dynamically renormalizes the effective fracture
energy, giving rise to a disorder-dependent limiting velocity $v_{{}_{\rm
lim}}$ that the crack cannot surpass. The emergence of $v_{{}_{\rm lim}}$ is
accompanied by a critical-like behavior in which the fracture energy features
a divergent variation with the crack velocity. In this regime, for
sufficiently large $G$, branches become so wide that the system statistically
recovers translational invariance along the crack front, and branching becomes
macroscopic.

This physical picture and its manifestations are consistent with a broad range
of experimental observations that were not fully understood previously. At the
macroscopic scale, it is widely observed (e.g., in [2, 15, 21]) that cracks in
brittle materials attain limiting velocities well below $c_{{}_{\rm R}}$, the
ideal theoretical limit [6, 7]. Moreover, the experiments on steady-state
cracks in a glassy polymer — reproduced in the inset of Fig. 2 — are
consistent with a critical-like behavior of the fracture energy upon
approaching $v_{{}_{\rm lim}}$. At smaller scales, localized instabilities in
3D, mainly in the form of micro-cracking damage [14, 15, 16, 17] and micro-
branching [21, 22, 23, 24] that involve a topological change and the
generation of extra fracture surface, are widely observed.

The localized branching instability extensively discussed above bears close
similarity to the micro-branching instability [21, 22, 23, 24]. The latter
features characteristic lengthscales at initiation, most notably minimal
micro-branch width and length [21, 22, 23, 24], consistently with our findings
in which the latter are related to the disorder correlation length $R$. Our
results show that localized branching is a probabilistic phenomenon, yet it
appears to be rather sharp due to the super-exponential increase of the
branching probability at small $G$, which depends on finite-strength disorder.
The probabilistic nature of the micro-branching instability and its dependence
on finite disorder are indeed experimentally demonstrated [24, 46], in
agreement with our findings.

With increasing driving force, micro-branches have been shown to significantly
increase in size, to lead to increased effective fracture energy associated
with extra fracture surface, to be accompanied by increased out-of-plane
fluctuations, to feature a 3D-to-2D transition [21] and eventually to
transform into macro-branches [23], all in agreement with our results. The
approximate linear scaling
$\langle\Delta{x}\rangle\\!\sim\\!\langle\Delta{z}\rangle$, cf. inset of Fig.
5, is also in agreement with experimental micro-branching observations [44].

At intermediate scales, commonly probed by fractography in experiments, our
findings are consistent with and offer novel understanding of the widely-
observed sequence of mirror-mist-hackle-macro-branching morphological
transitions [25, 26, 27, 28]. This sequence of transitions is commonly
observed as a function of crack propagation distance in experiments that
involve accelerating cracks. We take advantage of our computational framework,
which allows to simulate very long systems [29], to study this sequence of
transitions in a controlled manner at the steady-state level, as a function of
$G$. Not only our observations, cf. Fig. 3, unprecedentedly reproduce this
sequence of transitions as $G$ is varied, but we also provided a theoretical
understanding of the physics underlying each transition.

The mirror-mist transition is a probabilistic crossover, which is rather sharp
due to the super-exponential variation of the localized branching probability
for $G\\!<\\!G_{\rm B}$. Transiently-excited localized branching events in
this regime are spatially uncorrelated and small, featuring a scale comparable
to the disorder correlation length $R$, responsible for the misty appearance
of the fracture surface. The mist-hackle transition occurs at
$G\\!\simeq\\!G_{\rm B}$, upon which localized branching events become longer-
lived and feature dynamic lengthscales that increase with increasing $G$. The
wider and longer branching events are responsible for the hackle-like
appearance of the fracture surface. Finally, at sufficiently large $G$,
branches can capture the entire system thickness $L_{z}$, statistically
recovering translational symmetry in the thickness direction. This 3D-to-2D
transition results in macro-branching.

Our findings also raise some pressing questions for future investigation. A
central element in the emerging physical picture is finite-amplitude quenched
disorder. In the absence of well-established data about spatial disorder in
the fracture energy of amorphous materials, we invoked in our calculations a
minimal assumption in the form of a Gaussian disorder characterized by a
dimensionless strength $\sigma$ and correlation length $R$. It would be most
desirable to quantitatively characterize the statistical properties of various
important physical quantities in amorphous materials, including the amplitude
of fluctuations and their spatial correlation length. A possible indirect way
to address this point might be to study surface roughness at low propagation
velocities [3, 34].

Furthermore, we have shown that various quantities, e.g., the onset of
localized branching velocity and the limiting velocity, are disorder
dependent. Consequently, extracting the properties of the disorder may lead to
a quantitative comparison with experiments. Finally, while our large-scale
calculations have been carried out using cutting-edge GPU-based computational
platforms [29], we cannot yet rule some finite-size effects in our results,
especially with respect to $L_{y}$, but possibly also to $L_{z}$. Future
research, potentially employing enhanced computational power, should shed
light on this issue as well.

Acknowledgements. This work has been supported by the United States-Israel
Binational Science Foundation (BSF, grant no. 2018603). E.B. acknowledges
support from the Ben May Center for Chemical Theory and Computation, and the
Harold Perlman Family.

## References

  * [1] Ravi-Chandar, K.  Dynamic fracture of nominally brittle materials.  _International Journal of Fracture_ 90, 83–102 (1998).
  * [2] Fineberg, J. & Marder, M.  Instability in dynamic fracture.  _Physics Reports_ 313, 1–108 (1999).
  * [3] Bonamy, D. & Bouchaud, E.  Failure of heterogeneous materials: A dynamic phase transition?  _Physics Reports_ 498, 1–44 (2011).
  * [4] Bouchbinder, E., Goldman, T. & Fineberg, J.  The dynamics of rapid fracture: instabilities, nonlinearities and length scales.  _Reports on Progress in Physics_ 77, 046501 (2014).
  * [5] Fineberg, J. & Bouchbinder, E.  Recent developments in dynamic fracture: Some perspectives.  _International Journal of Fracture_ 196, 33–57 (2015).
  * [6] Freund, L. B.  _Dynamic Fracture Mechanics_ (Cambridge University Press, Cambridge, 1990).
  * [7] Broberg, K. R.  _Cracks and Fracture_ (Academic Press, New York, 1999).
  * [8] Livne, A., Ben-David, O. & Fineberg, J.  Oscillations in rapid fracture.  _Physical Review Letters_ 98, 124301 (2007).
  * [9] Lubomirsky, Y., Chen, C.-H., Karma, A. & Bouchbinder, E.  Universality and stability phase diagram of two-dimensional brittle fracture.  _Physical Review Letters_ 121, 134301 (2018).
  * [10] Bouchbinder, E., Livne, A. & Fineberg, J.  Weakly nonlinear theory of dynamic fracture.  _Physical Review Letters_ 101, 264302 (2008).
  * [11] Bouchbinder, E.  Dynamic crack tip equation of motion: High-speed oscillatory instability.  _Physical Review Letters_ 103, 164301 (2009).
  * [12] Chen, C.-H., Bouchbinder, E. & Karma, A.  Instability in dynamic fracture and the failure of the classical theory of cracks.  _Nature Physics_ 13, 1186–1190 (2017).
  * [13] Vasudevan, A., Lubomirsky, Y., Chen, C.-H., Bouchbinder, E. & Karma, A.  Oscillatory and tip-splitting instabilities in 2D dynamic fracture: The roles of intrinsic material length and time scales.  _Journal of the Mechanics and Physics of Solids_ 151, 104372 (2021).
  * [14] Ravi-Chandar, K. & Knauss, W.  An experimental investigation into dynamic fracture: II. Microstructural aspects.  _International Journal of fracture_ 26, 65–80 (1984).
  * [15] Ravi-Chandar, K. & Knauss, W. G.  An experimental investigation into dynamic fracture: III. On steady-state crack propagation and crack branching.  _International Journal of fracture_ 26, 141–154 (1984).
  * [16] Scheibert, J., Guerra, C., Célarié, F., Dalmas, D. & Bonamy, D.  Brittle-quasibrittle transition in dynamic fracture: An energetic signature.  _Physical Review Letters_ 104, 045501 (2010).
  * [17] Guerra, C., Scheibert, J., Bonamy, D. & Dalmas, D.  Understanding fast macroscale fracture from microcrack post mortem patterns.  _Proceedings of the National Academy of Sciences_ 109, 390–394 (2012).
  * [18] Tanaka, Y., Fukao, K., Miyamoto, Y. & Sekimoto, K.  Discontinuous crack fronts of three-dimensional fractures.  _EPL (Europhysics Letters)_ 43, 664 (1998).
  * [19] Kolvin, I., Cohen, G. & Fineberg, J.  Topological defects govern crack front motion and facet formation on broken surfaces.  _Nature Materials_ (2017).
  * [20] T. Baumberger, D. M. O., C. Caroli & Ronsin.  Magic angles and cross-hatching instability in hydrogel fracture.  _Physical Review Letters_ 100, 178303 (2008).
  * [21] Sharon, E. & Fineberg, J.  Microbranching instability and the dynamic fracture of brittle materials.  _Physical Review B_ 54, 7128 (1996).
  * [22] Sharon, E. & Fineberg, J.  Universal features of the microbranching instability in dynamic fracture.  _Philosophical Magazine B_ 78, 243–251 (1998).
  * [23] Sharon, E. & Fineberg, J.  The dynamics of fast fracture.  _Advanced Engineering Materials_ 1, 119–122 (1999).
  * [24] Livne, A., Cohen, G. & Fineberg, J.  Universality and hysteretic dynamics in rapid fracture.  _Physical review letters_ 94, 224301 (2005).
  * [25] Lawn, B.  _Fracture of brittle solids_ (Cambridge University Press, 1993).
  * [26] Johnson, J. & Holloway, D.  Microstructure of the mist zone on glass fracture surfaces.  _The Philosophical Magazine: A Journal of Theoretical Experimental and Applied Physics_ 17, 899–910 (1968).
  * [27] Rabinovitch, A., Belizovsky, G. & Bahat, D.  Origin of mist and hackle patterns in brittle fracture.  _Physical Review B_ 61, 14968 (2000).
  * [28] Jiao, D., Qu, R. T. & Zhang, Z. F.  Macroscopic bifurcation and fracture mechanism of polymethyl methacrylate.  _Advanced Engineering Materials_ 17, 1454–1464 (2015).
  * [29] _Summplemental Materials._
  * [30] Karma, A., Kessler, D. & Levine, H.  Phase-field model of mode III dynamic fracture.  _Physical Review Letters_ 87, 45501 (2001).
  * [31] Karma, A. & Lobkovsky, A. E.  Unsteady crack motion and branching in a phase-field model of brittle fracture.  _Physical Review Letters_ 92, 245510 (2004).
  * [32] Hakim, V. & Karma, A.  Laws of crack motion and phase-field models of fracture.  _Journal of the Mechanics and Physics of Solids_ 57, 342–368 (2009).
  * [33] Das, S., Lubomirsky, Y. & Bouchbinder, E.  Dynamics of crack front waves in three-dimensional material failure.  _Physical Review E_ 108, L043002 (2023).
  * [34] Ponson, L.  _Fracture Mechanics of Heterogeneous Materials: Effective Toughness and Fluctuations_ , pp. 207–254 in Mechanics and Physics of Fracture: Multiscale Modeling of the Failure Behavior of Solids (Springer International Publishing, Cham, 2023).
  * [35] Roch, T., Lebihain, M. & Molinari, J.-F.  Dynamic crack-front deformations in cohesive materials.  _Physical Review Letters_ 131, 096101 (2023).
  * [36] Lebihain, M., Leblond, J.-B. & Ponson, L.  Effective toughness of periodic heterogeneous materials: the effect of out-of-plane excursions of cracks.  _Journal of the Mechanics and Physics of Solids_ 137, 103876 (2020).
  * [37] Eshelby, J.  _Energy relations and the energy-momentum tensor in continuum mechanics_ , pp. 77–115 in Inelastic behaviour of solids (McGraw-Hill, New York, 1970).
  * [38] Katzav, E., Adda-Bedia, M. & Arias, R.  Theory of dynamic crack branching in brittle materials.  _International Journal of Fracture_ 143, 245–271 (2007).
  * [39] Kobayashi, A., Wade, B., Bradley, W. & Chiu, S.  Crack branching in homalite-100 sheets.  _Engineering Fracture Mechanics_ 6, 81–92 (1974).
  * [40] Sun, Y., Edwards, M. G., Chen, B. & Li, C.  A state-of-the-art review of crack branching.  _Engineering Fracture Mechanics_ 257, 108036 (2021).
  * [41] Yoon, J. _et al._ In situ tensile and fracture behavior of monolithic ultra-thin amorphous carbon in tem.  _Carbon_ 196, 236–242 (2022).
  * [42] Wagner, H. _et al._ Local elastic properties of a metallic glass.  _Nature Materials_ 10, 439–442 (2011).
  * [43] Kapteijns, G., Richard, D., Bouchbinder, E. & Lerner, E.  Elastic moduli fluctuations predict wave attenuation rates in glasses.  _The Journal of Chemical Physics_ 154, 081101 (2021).
  * [44] Sharon, E., Cohen, G. & Fineberg, J.  Crack front waves and the dynamics of a rapidly moving crack.  _Physical Review Letters_ 88, 085503 (2002).
  * [45] Fineberg, J., Sharon, E. & Cohen, G.  Crack front waves in dynamic fracture.  _International Journal of Fracture_ 121, 55–69 (2003).
  * [46] Boué, T. G., Cohen, G. & Fineberg, J.  Origin of the microbranching instability in rapid cracks.  _Physical Review Letters_ 114, 054301 (2015).


# Performance of a plastic scintillator developed using styrene monomer
polymerization

A. Sadigov 11footnotemark: 122footnotemark: 2 F. Ahmadov 33footnotemark:
344footnotemark: 4<EMAIL_ADDRESS>G. Ahmadov 55footnotemark:
566footnotemark: 677footnotemark: 788footnotemark: 8 E. Aksu 99footnotemark: 9
D. Berikov 1010footnotemark: 101111footnotemark: 11 S. Nuruyev
1212footnotemark: 121313footnotemark: 13 R. Akbarov 1414footnotemark:
141515footnotemark: 151616footnotemark: 16 M. Holik 1717footnotemark:
171818footnotemark: 18 J. Nagiyev 1919footnotemark: 19 S. Gürbüz Güner
2020footnotemark: 20 A. Mammadli 2121footnotemark: 21 N. Suleymanova
2222footnotemark: 22 C. Abbasova 2323footnotemark: 23 S. Melikova
2424footnotemark: 24 E. Yilmaz 2525footnotemark: 25 O. Tagiyev
2626footnotemark: 26 S. Lyubchyk 2727footnotemark: 27 Z. Sadygov
2828footnotemark: 282929footnotemark: 293030footnotemark: 30

###### Abstract

This paper presents a newly developed plastic scintillator produced in
collaboration with Türkiye Energy, Nuclear and Mineral Research Agency
(TENMAK). The scintillator is manufactured using thermal polymerization of
commercially available styrene monomer. The absorption spectrum of the
scintillator exhibited two absorption bands at 225 nm and 340 nm, with an
absorption edge observed at 410 nm. The wavelength of the emitted light was
measured in the range of 400-800 nm, with a maximum intensity at 427 nm.
Monoenergetic electrons from the 137Cs source were used to evaluate the
characteristics of the new scintillator, particularly its light yield. As the
light readout the MAPD-3NM type silicon photomultiplier array (4 $\times$ 4)
with an active area of 15 $\times$ 15 mm2, assembled using single MAPDs with
an active area of 3.7 $\times$ 3.7 mm2, was used. The light yield of the
scintillator was determined to be 6134 photons/MeV. In addition, the
efficiency of the scintillator for gamma rays with an energy of 662 keV was
found out to be approximately 1.8 $\%$. A CmBe neutron source was employed to
evaluate its fast neutron detection performance. However, neutron/gamma
discrimination using pulse shape discrimination (charge integration) method
was not observed. The results demonstrate the potential of a newly produced
plastic scintillator for various applications, particularly in radiation
monitoring and detection systems.

###### keywords:

Micropixel avalanche photodiode , SiPM , plastic scintillator , gamma source ,
CmBe neutron source , styrene monomer , styrene monomer polymerization

††journal: Radiation Measurements

[fn1]organization=Institute of Radiation Problems under Ministry of Science
and Education,addressline=B.Vahabzade str. 9, city=Baku, postcode=AZ1143,
country=Azerbaijan

[fn2]organization=Innovation and Digital Development Agency Nuclear Research
Department,addressline=Gobu str. 20th km of Baku-Shamakhi Highway, city=Baku,
postcode=AZ1021, country=Azerbaijan

[fn3]organization=Innovative Electronics and Detectors
LLC,addressline=Badamdard STQ-1, city=Baku, postcode=AZ1021,
country=Azerbaijan

[fn4]organization=Joint Institute for Nuclear Research,addressline=Joliot-
Curie 6, city=Dubna, postcode=141980, country=Russia

[fn5]organization=Türkiye Energy, Nuclear and Mineral Research Agency, Nuclear
Energy Research Institution,addressline=Saray, Atom Cd. 27, city=Ankara,
postcode=06980, country=Türkiye

[fn6]organization=Institute of Nuclear Physics of the Ministry of Energy of
the Republic of Kazakhstan,addressline=Ibragimova 1, city=Almaty,
postcode=050032, country=Kazakhstan

[fn7]organization=Faculty of Electrical Engineering, UWB in
Pilsen,addressline=Univerzitni 2795/26, city=Pilsen, postcode=306 14 ,
country=Czech Republic

[fn8]organization=Institute of Experimental and Applied Physics, CTU in
Prague,addressline=Husova 240/5, city=Prague, postcode=110 00 , country=Czech
Republic

[fn9]organization=Department of Physics, Bolu Abant Izzet Baysal
University,addressline=Gölköy Kampüsü, city=Bolu, postcode=14280,
country=Türkiye

[fn10]organization=Institute of Physics – Ministry of Science and
Education,addressline=G.Javid ave. 131, city=Baku, postcode=AZ1143,
country=Azerbaijan

[fn11]organization=Universidade Lusofona,addressline=Campo Grande 376,
city=Lisbon, postcode=1749-024, country=Portugal

##  1 Introduction

Plastic scintillators with photo sensors have found wide application in high-
energy physics, space exploration, medical diagnostics, and security systems
[1, 2, 3, 4, 5]. The widespread use of plastic scintillators can be explained
by advances in their manufacturing technology as well as a number of valuable
properties of the scintillator itself, such as short decay time, high
radiation resistance, operating temperature, resistance to mechanical stress,
etc. Moreover, the ability to produce plastic scintillators in various shapes
and sizes makes them suitable converters of ionizing energy into visible
light, rendering them suitable for numerous experiments [4, 5, 6]. In light of
these factors, the production of plastic scintillators continues to be of
utmost significance, driving ongoing research and development efforts aimed at
further improving their performance, cost-effectiveness, and compatibility
with emerging detection systems.

One notable advantage of plastic scintillators is their exceptionally short
decay time, which makes them indispensable for experiments that require
precise time-of-flight measurements. Additionally, due to the low atomic
number, they are effective in detecting charged particles and neutrons [6].
The interaction of gamma rays with plastic scintillators is primarily related
to their atomic number (Z) and energy of gamma ray ($E_{\gamma}$), resulting
in a minimal probability of the photoelectric effect
($\sigma_{pe.}\sim\frac{Z^{n}}{E_{\gamma}^{3.5}}$). Consequently, the energy
of gamma rays in plastic scintillators is typically determined by the Compton
edge ($\sigma_{comp.}\sim const.\cdot Z$) and the maximum energy of a Compton
electron can be calculated using the following formula:

$E_{e}=\frac{2\cdot E_{\gamma}^{2}}{2\cdot E_{\gamma}+m_{e}\cdot c^{2}}$ (1)

where $E_{e}$ \- the maximum energy of a Compton electron, $m_{e}$ \- the mass
of an electron and $c$ \- the speed of light [7, 8]. Table 1 lists the maximum
energy of Compton electrons produced by various gamma rays.

Table 1: The maximum energy of Compton electrons produced by various gamma rays. Radionuclide | Energy of gamma-rays (keV) | Maximum energy of Compton electrons (keV)
---|---|---
241Am | 59.6 | 11.27
133Ba | 81, 276, 303, 356, 384 | 19.5, 143.3, 164.4, 207.2, 230.6
137Cs | 661.6 | 477.65
60Co | 1173, 1332 | 963.2, 1117.6

When ionizing radiation interacts with the scintillator material, the
scintillator molecules get excited, and then they return to their ground state
with emitting scintillation photons. All absorbed energy of ionizing radiation
in the scintillator is not converted to visible scintillation light. Ideally,
the number of photons generated should increase linearly with the incident
radiation energy, although deviations from linearity can occur depending on
the specific characteristics of the scintillation material. The light yield of
the scintillator is calculated as follows [4]:

$Y=\frac{E_{dep.}\cdot\eta}{E(\lambda)}$ (2)

where $E_{dep.}$ \- deposited energy of ionizing radiation, $\eta$ \-
efficiency of scintillator and $E(\lambda)$ \- the energy of an emitted
photon. The light yield of the scintillator depends on the scintillation
material, the type of incident particles, the energy of particles and the
temperature [5].

In scintillation detectors, the phenomenon of quenching plays a significant
role in the detection and measurement of ionizing radiation. Quenching occurs
due to various mechanisms associated with the interaction of radiation with
the scintillator material. Scintillator quenching is generally described by
different models (Birks, Yoshida, Voltz) [9, 10]. In these models are taken
into account affect electronic stopping, nuclear stopping, the singlet
quenching, the triplet quenching and others effects. The luminescence yield
per unit length in an organic scintillator ($\frac{dY}{dx}$) can be calculated
given the ionization density ($\frac{dE}{dx}$), the energy transfer
probability (k), a normalization factor (s), the constant of proportionality
(B) to determine the number of damaged molecules, the energy emitted as light
(Y), and particle energy dissipated in the scintillator (E)

$\frac{dY}{dx}=\frac{s\cdot\frac{dE}{dx}}{1+kB\cdot\frac{dE}{dx}}.$ (3)

Different types of particles, such as electrons, alpha particles, heavy ions,
can exhibit distinct quenching behaviors. Quenching of the light yield is
related to the energy transfer to damaged molecules, which do not convert the
received energy into scintillation photons. These features, which leads to
variations in the shape of the detected pulses, can be employed to distinguish
different types of ionizing radiation [5, 9]. It is possible to discriminate
particles type by analyzing pulse shapes, thereby enhancing the capabilities
of radiation detection systems. This performance of plastic scintillation
detectors allows use them in PSD experiments [4].

The light yield of the using plastic scintillators mainly changes in the range
of 3000 - 20000 photons/MeV e.g : EJ-232-0.5$\%$ (2900 photons/MeV),
EJ-256-5$\%$ (5200 photons/MeV), EJ-240G (6300 photons/MeV), EJ-254-5$\%$
(7500 photons/MeV), EJ-276G (8000 photons/MeV), EJ-290 (9000 photons/MeV),
EJ-212 (10000 photons/MeV) and others [11, 12].

The main objective of this study is to investigate the parameters and
performance of the plastic scintillator produced by TENMAK. The absorption and
photoluminescence performance of this scintillator were measured, and the
light yield was determined using the light output of the commercial standard
$LaBr_{3}(Ce)$ scintillator by Epic Crystal. To gain insight into the
scintillation mechanism, the optical performance, scintillation responses, and
pulse shape discrimination (PSD) performance on $\gamma$-rays, neutrons, and
beta particles were investigated using MAPD-3NM type silicon photomultiplier
as a light readout and various sources including 133Ba, 137Cs, and 60Co, 90Sr,
and 244CmBe neutron source.

##  2 Experimental

In this study, we used a polystyrene plastic scintillator with a rectangular
shape and dimensions of 15 $\times$ 15 $\times$ 50 mm3 to convert the energy
of ionizing radiation (gamma ray and beta particle) to visible light pulse.
The plastic scintillator was produced through thermal polymerization of
commercially available styrene monomer. The production process involved
several steps. Firstly, purified styrene monomer was heated to 120 ∘C over a
period of 5 - 10 hours. It was maintained at this temperature for 3 days to
ensure complete polymerization. Afterward, the material was slowly cooled to a
temperature below the glass transition temperature of polystyrene (Tg = 100
∘C) and underwent annealing to relieve internal stresses and avoid cracking.
Finally, the scintillator was cut into the desired shape from the fabricated
piece, and the surfaces were polished to achieve a mirror-like finish. To
assess the properties of the scintillator, including scintillation properties
and density, measurements were conducted.

Silicon photomultiplier (SiPM) array (type MAPD-3NM) was used as a detector to
read out scintillation photons produced by ionizing radiation in the
scintillator (figure 1 left). MAPD-type SiPMs have been designed and developed
by our research group, employing a deeply buried pixel design. This design
incorporates a double n-p-n-p junction with micro-well structures located
below the surface. The structure of these SiPMs does not include conventional
quenching resistors, in contrast to SiPMs produced using standard surface-
pixel technology. Triggered avalanche is quenched through the utilization of
specially designed potential barriers. For a comprehensive understanding of
the structure and operational principles of MAPD-type SiPMs, refer to [13, 14,
15, 16, 17, 18]. Additional insights into the performance of MAPD as both an
individual device and an array, with various scintillation materials, can be
found in references [19, 20, 21, 22, 23].

The plastic scintillator was covered with teflon layers and connected to the
photodiode via optical glue (figure 1 right). The SiPM array used consisted of
16 elements, with each element being connected in parallel. A voltage was
applied to SiPM array via red wire, which was linked to the cathode and the
signal was taken from the black wire, connected to the anode (figure 1 left).
The SiPM array has the following parameters: size - 15 $\times$ 15 mm2,
operation voltage - 55.2 V, capacitance - 155 pF, breakdown voltage - 51 V,
average gain - 2 $\times$ $10^{5}$, photon detection efficiency (PDE)- 30-35
$\%$.

Figure 1: Photo of the electronic system and the detector based on MAPD and
plastic scintillator.

The ionizing radiation measurements were carried out with the Spectrig MAPD
device. Detailed information about the Spectrig MAPD device can be found in
[22, 25, 26]. During measurement with the Spectrig MAPD, the following
parameters were selected: gate width - 180 ns, variable gain - 1 - 15 dB,
threshold - 43 mV, bias voltage - 55.2 V and measurement time - 200 sec and
500 sec. All measurements were carried out at a temperature of 22 ∘C. Figure 2
presents a photo of the experimental setup.

Figure 2: Photo of the experimental setup.

##  3 Result and discussion

Figure 3 shows the absorption and photoluminescence spectra of the plastic
scintillators. The absorption and photoluminescence spectra were recorded by
LS 55 Fluorescence Spectrometer (Perkin Elmer) [23] and Cary 50Scan UV-Vis
Spectrophotometer (Varian) [24], respectively. The absorption spectrum was
measured in the wavelength range of 200-800 nm.

Figure 3: The absorption and photoluminescence spectra of the plastic
scintillator.

The absorption spectrum exhibits two absorption bands at 225 and 340 nm, with
an absorption edge observed at 410 nm (figure 3). As shown in figure, the
transmission of light through the scintillators decreases in the region of
400-800 nm, with the transmission loss increasing as the radiation wavelength
increases up to 550 nm. Beyond 550 nm, no significant changes are observed in
the spectra. The photoluminescence spectrum of the scintillator was measured
using a xenon lamp with a continuous spectrum ranging from 230 to 50 nm at
room temperature. The spectrum exhibits a wide band covering the wavelength
range of 350-550 nm, with a maximum corresponding to 427 nm (figure 3).

Energy spectra of 137Cs measured with the SiPM array + $LaBr_{3}(Ce)$ [25, 26]
and the SiPM array + plastic scintillator were shown in figure 4. The purpose
of using the $LaBr_{3}(Ce)$ scintillator was to calculate the light yield of
the plastic scintillator. It is well known that 137Cs decays by beta emission
to a metastable state (137mBa), and directly the ground state of 137Ba, with
both cases resulting in the emission of gamma rays ($\sim$ 662 keV). In the
first case, there is approximately a 9.6 $\%$ probability that the emitted
gamma rays are captured by K- shell electrons of 137Ba, generating a
monoenergetic beta particle with an energy of 626 keV. The gamma ray losses
about 36 keV of energy due to the binding energy of 137Ba’s K-shell electrons
[3].

Figure 4: Energy spectra of 137Cs measured with the SiPM array +
$LaBr_{3}(Ce)$ and the SiPM array+plastic scintillator (the variable gain was
selected 1 dB).

The amplitude of the signal can be calculated as $A=PDE\times M\times N$,
where N (related to scintillator) is the number of scintillation photons
produced by 1 MeV of ionization radiation (or light yield of scintillator), M
is the gain of SiPM, and PDE (related to SiPM) is the photon detection
efficiency, which depends on the over-voltage and wavelength of the photon
[5]. Considering that the photon detection efficiency and gain of the SiPM are
about the same for both emission wavelengths of scintillators, the amplitude
of the signal will depend on the light yield of the scintillators. The light
yield of the $LaBr_{3}(Ce)$ scintillator is determined to be 68000 photons per
MeV of incident energy, corresponding to approximately 45016 photons for the
incident gamma ray with an energy of 662 keV [25, 26]. The detected gamma ray
signal with an energy of 662 keV, using the $LaBr_{3}(Ce)$ scintillator, has
an amplitude of 2122 ADC channel. In this case, 626 keV energy will
approximately correspond to the 2006.6 ADC channel, and the number of
scintillation photons will be approximately 42568 photons. Conversely, when
using the plastic scintillator, the amplitude of the detected mono-energetic
beta particle signal with an energy of 626 keV is measured to be 181 ADC
channel. Considering this information one can determine the light yield of the
plastic scintillator using the following ratio:

$\frac{A_{LaBr_{3}(Ce)}(626\,keV)}{N_{LaBr_{3}(Ce)}(626\,keV)}=\frac{A_{pl}(626\,keV)}{N_{pl}(626\,keV)}$
(4)

$N_{pl}\,(626keV)=N_{LaBr_{3}(Ce)}\,(626\,keV)\times\frac{A_{pl}\,(626\,keV)}{A_{LaBr_{3}(Ce)}}=42568\,ph.\times\frac{181\,ADC}{2006.6\,ADC}=3840\,ph.$
(5)

In the case of the plastic scintillator, the number of scintillation photons
corresponding to the incident gamma ray with an energy of 626 keV was found to
be 3840. The light yield of the plastic scintillator for 1 MeV of incident
energy was determined to be 6134 photons. It is important to note that in the
given calculation, factors such as light loss due to reflector materials,
linearity of scintillators and the dependence of PDE on wavelength were not
taken into account. Although in the works of the authors [27, 28] on the study
of the non-proportionality of electron response in scintillators, it can be
seen that for some plastic scintillators in the region of 626 keV, the
electron response and relative light yield are approximately 100 $\%$. The
efficiency ($\eta$) of the scintillator can be calculated as the following
formula [4]:

$\eta=\frac{E_{ph}(\lambda)\times N_{ph}}{E_{\gamma}}$ (6)

where $E_{ph}$ \- corresponding energy to a scintillation photon, $N_{ph}$ \-
the number of scintillation photons produced by 1 MeV of ionization radiation,
and E - the energy of the absorbed gamma ray. In this experiment, the energy
corresponding to 427 nm was measured to be 2.9 eV, the light yield was
determined to be 6134 photons/MeV. The obtained efficiency of the scintillator
was approximately 1.8 $\%$. The energy resolution of the monoenergetic
electron (626 keV) measured with the plastic scintillator was 33.9 $\%$. On
the other hand, in the case of the measurement with $LaBr_{3}(Ce)$, the energy
resolution for 662 keV gamma rays was 3.3 $\%$. In these measurements, a
variable gain of MAPD Spectrig was selected 1 dB.

Figure 5 shows the energy spectra of 133Ba, 137Cs and 60Co measured using the
Spectrig MAPD with the variable gain of 15 dB. To distinguish beta particles
from gamma rays emitted from 137Cs, a copper foil with a thickness of 1 mm was
placed between the scintillator and the source. Consequently, the
monoenergetic (770-1200 ADC channel) and low-energy parts (200-730 ADC
channel) of the beta particles were effectively absorbed by the foil,
resulting in the observation of the Compton continuum and edge in the
spectrum. The measurement time was selected based on the activity of the
sources. The measurement time was selected 200 seconds for 137Cs source, while
500 seconds for 133Ba, 137Cs with the copper foil and 60Co.

Figure 5: Energy spectra of 133Ba, 137Cs and 60Co measured with the SiPM array
+ plastic scintillators (the variable gain was selected 15 dB).

Energy resolution of the monoenergetic electrons with an energy of 626 keV was
measured to be 25$\%$ at 15 dB variable gain. Notably, the Compton edge of
gamma rays was observed to increase with the energy of the incident gamma
rays. Figure 6 shows the energy spectrum obtained from 90Sr. 90Sr undergoes
decay, emitting an electron with a maximum energy of 546 keV, transforming
into 90Y, which further decays, emitting an electron with a maximum energy of
2274 keV, and finally resulting in 90Zr [29]. When a copper foil was placed
between the scintillator and the source, the beta particles were completely
absorbed by the foil, and only gamma and x-ray events were observed in the
spectrum. The obtained results were found to be in agreement with the data
reported in other works [29].

Figure 6: Energy spectra of 90Sr measured with the SiPM array + plastic
scintillators (the variable gain was selected 15 dB).

CmBe neutron sources were utilized to test the neutron detection performance
of the plastic scintillator. The energy spectrum of CmBe is presented in
figure 7. The CmBe source emits both neutrons and gamma rays, both of which
have energies within the range of up to 10000 keV. The detection of fast
neutrons is achieved through elastic collisions between neutrons and hydrogen
atoms, leading to the generation of protons. The newly generated protons,
resulting from the nuclear interaction with hydrogen, can be effectively
detected using the scintillator. In order to minimize the effects of gamma
rays, Pb shielding material was employed, which effectively attenuates gamma
rays while allowing neutrons to penetrate and interact with the scintillator.

Figure 7: Energy spectra of CmBe measured with the SiPM array + plastic
scintillators (the variable gain was selected 15 dB).

The source was placed inside Pb shielding box with a thickness of 5 cm, which
effectively reduced the intensity of gamma rays. This arrangement enabled the
observation of the effect of fast neutrons on the plastic scintillator.
However, due to the limited size of the scintillator, it was insufficient to
fully detect the high-energy part of neutrons, resulting in only a few events
being visible in the high-energy region of the spectra. The implementation of
the pulse shape and zero crossing discrimination methods [30] proved
ineffective in distinguishing between neutrons and gamma rays. Despite efforts
to utilize this technique, the response of the scintillator did not exhibit
the required differentiation between the two types of radiation. As a result,
the pulse shape and zero crossing discrimination methods were not effective
for distinguishing neutrons from gamma rays in this particular scintillator
setup. Further investigations and alternative methods may be necessary to
achieve the desired discrimination capability in future studies.

##  4  Conclusion

The gamma ray, beta particle and neutron detection performance of the newly
fabricated plastic scintillator by TENMAK was investigated. The plastic
scintillator exhibits two absorption bands at 225 and 340 nm, with an
absorption edge observed at 410 nm. The wavelength of maximum emission for the
plastic scintillator is $\sim$ 427 nm, with the light yield of 6134 photons
for 1 MeV incident radiation. The obtained results showed that the Compton
edge of detected gamma rays was observed to increase linearly with increasing
the energy of gamma rays. The plastic scintillator demonstrated effective
detection of beta particles emitted by 137Cs and 90Sr sources, with an energy
resolution of 25 $\%$ for monoenergetic (626 keV) electrons. Furthermore, this
type of plastic scintillator exhibited sensitivity to gamma rays and fast
neutrons. The scintillator did not demonstrate the capability to discriminate
between neutrons and gamma rays using the pulse shape discrimination method.
Considering these advantageous characteristics, the combination of this
plastic scintillator with SiPM technology makes it well-suited for detecting
beta particles, gamma rays and neutrons in radiation monitoring and
applications related to public safety. Further investigations are essential to
attain optimal light yields and enhance the properties of the scintillator.

## ACKNOWLEDGMENTS

This project has received funding from the European Union’s Horizon 2021
Research and Innovation Programme under the Marie Sklodowska-Curie grant
agreement 101086178.

## References

  * [1] E. Saito et al., _Light yield, long-term stability, and attenuation length of a new plastic scintillator cured at room temperature_ , _Nucl. Instrum. Meth. A_ 953 (2020) 162885.
  * [2] S. Min et al., _Optimization of Plastic Scintillator for Detection of Gamma-Rays: Simulation and Experimental Study_ , _Chemo sensors_ 9 (2021) 239.
  * [3] F. Ahmadov, _Investigation of Their Physical Properties and Application in Radiation Spectrometer_ , Doctoral thesis in Radiation materials science Baku (2022).
  * [4] Syed Naeem Ahmed, _Physics and Engineering of Radiation Detection_ , Second Edition-London: ELSEVIER, (2015).
  * [5] G.F. Knoll, _Radiation Detection and Measurement_ , fourth edition, Wiley Press (2010).
  * [6] S. Berns et al., _A novel polystyrene-based scintillator production process involving additive manufacturing_ , 2020 _JINST_ 15 P10019.
  * [7] J. Kenneth, E. Faw Richard, _Fundamentals of Nuclear Science and Engineering_ , Second Edition - New York: CRC Press (2007).
  * [8] M. Rizzi , D. Aloia and B. Castagnolo, _Semiconductor Detectors and Principles of Radiation-matter Interaction_ , _Journal of Applied Sciences_ 10 (2010) p. 3141-3155.
  * [9] Ogawa Tatsuhiko, Yamaki Tetsuya, Sato Tatsuhiko, _Analysis of scintillation light intensity by microscopic radiation transport calculation and Förster quenching model_ , _PLoS ONE_ 2018 13(8) :e0202011.
  * [10] Thibault A. Laplace et al., _Modeling ionization quenching in organic scintillators_ , _Mater. Adv._ 3 (2022) 5871-5881.
  * [11] _https://eljentechnology.com/products/plastic-scintillators_
  * [12] E. V. van Loef, U. Shirwadkar, M. McClish and K. S. Shah, _Solution growth and scintillation properties of 9-phenylcarbazole_ , 2017 _J. Cryst. Growth_ 480 145\.
  * [13] A. Sadigov et al., _Improvement of parameters of micro-pixel avalanche photodiodes_ , 2022 _JINST_ 17 P07021.
  * [14] F. Ahmadov et al., _Investigation of parameters of new MAPD-3NM silicon photomultipliers_ , 2022 _JINST_ 17 C01001.
  * [15] S. Nuriyev et al., _Performance of a new generation of micropixel avalanche photodiodes with high pixel density and high photon detection efficiency_ , _Nucl. Instrum. Methods Phys. Res. A_ 912 (2018) 320.
  * [16] S. Nuruyev et al., _Performance of silicon photomultipliers at low temperature_ , 2020 _JINST_ 15 C03003.
  * [17] A. Sadigov et al., _A new detector concept for silicon photomultipliers_ , _Nucl. Instrum. Methods Phys. Res. A_ 824 (2016) 135.
  * [18] A. Sadigov et al., _A new method to improve multiplication factor in micro-pixel avalanche photodiodes with high pixel density_ , _Nucl. Instrum. Methods Phys. Res. A_ 824 (2016) 137.
  * [19] F. Ahmadov et al., _New phoswich detector based on LFS and p-terphenyl scintillators coupled to micropixel avalanche photodiode_ , _Func. materials_ 24 (2017) 341.
  * [20] F. Ahmadov et al., _Alpha particle detector based on micropixel avalanche photodiodes_ , _Phys. Part. Nucl. Lett._ 10 (2013) 778.
  * [21] G. Ahmadov et al., _Gamma-ray spectroscopy with MAPD array in the readout of LaBr3:Ce scintillator_ , 2021 _JINST_ 16 P07020.
  * [22] M. Holik et al., _Miniaturized read-out interface "Spectrig MAPD" dedicated for silicon photomultipliers_ , _Nucl. Instrum. Methods Phys. Res. A_ 978 (2020) 164440.
  * [23] _https://www.perkinelmer.com/_
  * [24] _https://www.agilent.com/_
  * [25] M. Holik et al., _Gamma ray detection performance of newly developed MAPD-3NM-II photosensor with LaBr3(Ce) crystal_ , _Scientific Reports_ 12 (2022) 15855.
  * [26] M. Holik et al., _Investigation of the possibility of a new detector based on SiPM in nuclear forensics_ , 2023 _JINST_ 18 C01015.
  * [27] P. Limkitjaroenporrn, J. Kaewkhao, P. Limsuwan, W. Chewpraditkul, _Nonproportionality of electron response using CCT: Plastic scintillator_ , _Applied Radiation and Isotopes_ 68 (2010) 1780.
  * [28] P. L. Swiderski, R. Marcinkowski, M. Moszynski et al., _Electron response of some low-Z scintillators in wideenergy range_ , _JINST_ 7 (2012) P06011.
  * [29] Graeme Turkington, Kelum A.A. Gamage, James Graham, _Beta detection of strontium-90 and the potential for direct in situ beta detection for nuclear decommissioning applications_ , _Nucl. Instrum. Methods Phys. Res. A_ 911 (2018) 55.
  * [30] R. Akbarov et al., _Fast neutron detectors with silicon photomultiplier readouts_ , _Nucl. Instrum. Methods A_ 936 (2019) 549.


11institutetext:  Dept. of Computer Science and Engineering, PES University,
Bengaluru, India
11email<EMAIL_ADDRESS>22institutetext: Dept. of Mathematics, IIT
Madras, Chennai, India
22email<EMAIL_ADDRESS>33institutetext: Dept. of Computer Science and
Engineering, Manipal University, Jaipur, India
33email<EMAIL_ADDRESS>44institutetext: Dept. of Data
Science and Engineering, Indian Institute of Science Education and Research,
Bhopal, India
44email<EMAIL_ADDRESS>
# HeCiX: Integrating Knowledge Graphs and Large Language Models for Biomedical
Research

Prerana Sanjay Kulkarni1* Muskaan Jain2* Disha Sheshanarayana3* Srinivasan
Parthiban4

###### Abstract

Despite advancements in drug development strategies, 90% of clinical trials
fail. This suggests overlooked aspects in target validation and drug
optimization. In order to address this, we introduce HeCiX-KG, Hetionet-
Clinicaltrials neXus Knowledge Graph, a novel fusion of data from
ClinicalTrials.gov and Hetionet in a single knowledge graph. HeCiX-KG combines
data on previously conducted clinical trials from ClinicalTrials.gov, and
domain expertise on diseases and genes from Hetionet. This offers a thorough
resource for clinical researchers. Further, we introduce HeCiX, a system that
uses LangChain to integrate HeCiX-KG with GPT-4, and increase its usability.
HeCiX shows high performance during evaluation against a range of clinically
relevant issues, proving this model to be promising for enhancing the
effectiveness of clinical research. Thus, this approach provides a more
holistic view of clinical trials and existing biological data.

###### Keywords:

Knowledge Graph Large Language Model LangChain Clinical Trials.

††footnotetext: * These authors contributed equally to this work.

##  1 Introduction

The pharmaceutical industry faces significant challenges in drug discovery,
with alarming clinical trial failure rates of almost 90% [1]. The rise in
attrition rates reflects not only huge financial losses but also delayed
implementation of life-saving treatments for patients.

One of the major reasons underlying this is the fragmented nature of the
available data. Hetionet [2] contains vast domain knowledge about diseases,
genes, and anatomy, yet it lacks sufficient information about previously
conducted clinical trials and experiments. Conversely, ClinicalTrials.gov [3]
houses extensive information about clinical trials and experiments conducted
worldwide, including details about Principal Investigators of studies which
can be useful in identifying Key Opinion Leaders (KOLs). However, it offers
limited insights into the diseases themselves. This disparity between our
understanding of fundamental biology and clinical trial outcomes hinders
effective drug development.

To address this, we propose a novel knowledge graph, HeCiX-KG (Hetionet-
Clinicaltrials neXus Knowledge Graph), which integrates information from
clinicaltrials.gov [3] and Hetionet [2]. HeCiX-KG is a single knowledge graph,
connecting biological knowledge with clinical trial data. This integration can
provide better understanding, revealing linkages and patterns that were
previously missed, but are vital for the effective repurposing and discovery
of new drugs.

Building upon HeCix-KG, we introduce HeCiX, a system that utilizes OpenAI’s
GPT-4 [4] and LangChain [5] to enable seamless interaction with the knowledge
graph. HeCiX translates natural language queries into CQL (Cypher Query
Language) queries, which makes it possible to efficiently retrieve relevant
context from the knowledge graph. Subsequently, the system displays the result
in human-understandable format, rendering the information available to
clinical and biomedical researchers.

We evaluate HeCiX’s performance against a wide array of question-answering
tasks relevant to the domain. The results show notable advancements in the
scope and depth of data obtained, thus providing a helpful tool to enhance the
efficiency of clinical research, and improve the success rates of drug
repurposing and development. HeCiX overcomes significant shortcomings in the
current resources by offering a holistic view of disease biology, clinical
trial history, and expert knowledge at the user’s fingertips.

The structure of the paper is as follows. We describe the background of the
work and the knowledge graph construction in Section 2 and Section 3
respectively. Section 4 talks about the detailed methodology. Experimentation
is described in Section 5 which is followed by results and discussions in
Section 6. Finally, Section 7 discusses the conclusion.

##  2 Background and Related Work

###  2.1 Knowledge Graphs in Biomedical Research

Knowledge graphs have emerged as powerful tools for representing and
integrating complex biomedical information, facilitating efficient data
integration and knowledge discovery. Notable examples include Bio2RDF [6],
CTKG by Chen et al. [7], Hetionet [2], and others. They have been used to
enhance drug discovery, understand disease mechanisms, and identify biomedical
relationships.

###  2.2 Large Language Models in Healthcare

Large Language Models (LLMs) have significantly improved medical literature
analysis, clinical note interpretation, and diagnosis support. Important
models in this field include DeepMind’s MedIC, Microsoft’s BioGPT [8],
TrialGPT [9], BioBart [10], and BioMistral [11], among others. They have
impacted domains such as drug discovery, patient data analysis, and clinical
decision support.

###  2.3 Hetionet: A Comprehensive Biomedical Knowledge Graph

Hetionet is a heterogeneous network of biomedical knowledge, integrating data
such as genes, compounds, diseases, and their interrelationships [2]. It has
been used in predicting drug-target interactions, identifying disease
mechanisms, and supporting drug repurposing efforts, making it ideal for
integration with other data sources to enhance clinical research capabilities.

###  2.4 ClinicalTrials.gov: A Repository of Clinical Trial Data

ClinicalTrials.gov is an extensive data source providing information about
clinical trials, studies conducted on various diseases, principal
investigators, study locations, and trial outcomes [3]. It thereby supports
clinical research and drug development.

##  3 Knowledge Graph Construction

HeCiX-KG is constructed from two primary sources of data, Hetionet [2] and
ClinicalTrials.gov [3]. It combines their data into a single knowledge source
and includes data related to six specific diseases, namely Vitiligo, Atopic
Dermatitis, Alopecia Areata, melanoma, Epilepsy, and Hypothyroidism.

###  3.1 Hetionet

Hetionet is a highly interconnected knowledge base, which combines data from
29 distinct databases. It comprises a total of 47,031 nodes across 11 types:
Disease, Compound, Gene, Symptom, Side Effect, Biological Process, Molecular
Function, Anatomy, Cellular Component, Pathway, and Pharmacologic Class [2].
For the purpose of constructing HeCiX-KG, we have extracted a subgraph of
Hetionet, consisting of data related to the six chosen diseases. This
comprises a total of 1071 nodes and 1125 relations.

###  3.2 ClinicalTrials.gov

ClinicalTrials.gov provides massive amounts of information about clinical
trials and studies on various diseases and conditions [3]. While the total
number of records in ClinicalTrials.gov exceeds 500,000, our research focuses
on a selected subset of 1,200 records, spanning the six selected diseases.
This subset when constructed as a knowledge graph, contains 5,454 nodes and
11,466 edges. The nodes in this subset are classified into 9 types: Disease,
Principal Investigators (PI), Study, Conditions, Phases, Locations,
Interventions, Age, and Sex. There are 10 types of relationships connecting
these nodes.

###  3.3 Schema

By taking inspiration from the schemas of the knowledge graphs constructed in
Hetionet [2] and the work of Devarakonda et al. [12], we have constructed a
comprehensive schema for HeCiX-KG. The ‘Disease’ node serves as the main
connecting point between the two integrated databases. This schema has been
illustrated in Figure 1. By populating the schema with our data, we have
obtained a knowledge graph consisting of 6,509 nodes and 14,377 edges.

Figure 1: HECIX-Schema

##  4 Methodology

Our approach consists of two major stages, construction of HeCiX-KG, and its
integration with GPT-4 using LangChain.

###  4.1 Knowledge Graph Construction

HeCiX-KG is constructed by extracting and integrating relevant data from
Hetionet [2] and ClinicalTrials.gov [3] for six specific diseases. The
resulting knowledge graph has 6,509 nodes and 14,377 edges. The construction
process involves data extraction, schema design, entity-relationship mapping,
and graph population.

###  4.2 LLM Integration using LangChain

To enhance the usability of HeCiX-KG, we developed HeCiX, a system that
integrates our knowledge graph with GPT-4 using LangChain. Specifically, we
utilized the GraphCypherQAChain component from the LangChain ecosystem for
this integration. As indicated in Figure 2, our query processing pipeline is
as follows:

  1. 1.

User Query Input: A user submits a natural language prompt to LangChain.

  2. 2.

Query and Prompt Processing: The user’s question is combined with a set prompt
template, and then sent to GPT-4

  3. 3.

Cypher Query Generation: GPT-4 generates a Cypher query based on the user’s
input and sends it back to LangChain.

  4. 4.

Database Querying: LangChain executes the generated Cypher query on HeCiX-KG.

  5. 5.

Raw Results Retrieval: HeCiX-KG returns the raw query results (the ‘Full
Context’) to LangChain.

  6. 6.

Context Forwarding: LangChain forwards the full context to GPT-4 for
interpretation and conversion into a human-readable format.

  7. 7.

Human-Readable Response Generation: GPT-4 generates a human-readable response
based on the full context sent to it, and sends it to LangChain.

  8. 8.

User Response: Finally, LangChain returns the human-readable response to the
user, thereby providing the user with the answer to their query.

Figure 2: Query processing pipeline

This integration allows users to interact with the complex HeCiX-KG using
natural language queries, significantly enhancing its accessibility and
usability for clinical researchers.

##  5 Experimentation

###  5.1 Experimental Setup

Our experimental setup consisted of the following major steps:

  1. 1.

We set up our AuraDB [13] server to host HeCiX-KG.

  2. 2.

We constructed the schema of our knowledge graph, based on the individual
structures of Hetionet and clinical trials data.

  3. 3.

We populated the schema with data from both Hetionet and ClinicalTrials.gov to
create HeCiX-KG.

  4. 4.

We integrated HeCiX-KG with GPT-4 using LangChain’s GraphQueryQAChain
component.

###  5.2 Evaluation Methodology

To properly assess and analyse the question-answering abilities of our system,
we designed a set of question-answering tasks for HeCiX to answer. These tasks
were carefully crafted to cover various aspects of clinical research,
including drug discovery, identification of KOLs, and analysis of biomedical
data, among others.

##  6 Results and Discussion

###  6.1 Experimentation Results

We utilized the RAGAS [14] framework to evaluate our model’s performance. The
framework calculates several key metrics: faithfulness, answer relevance,
context precision, and context recall. The results of our experimentation have
been displayed in Table 1.

Table 1: Performance metrics of HeCiX using the RAGAS framework Metric | Score
---|---
Faithfulness | 0.8572
Answer Relevance | 0.9340
Context Precision | 0.9202
Context Recall | 0.6654

###  6.2 Comparative Analysis

In addition to evaluating our model on the evaluation metrics seen earlier, we
also carried out a comparative analysis with OpenAI’s GPT-4 and Anthropic’s
Claude 3 Sonnet [15] in a zero-shot scenario. We tested these models on the
same set of questions used to evaluate HeCiX. We discovered that while these
models excel at answering questions based on domain-specific knowledge, they
struggle with queries related to clinical trials.

Our experiments with HeCiX demonstrated promising results in enhancing
clinical research capabilities. Figure 3 shows a sample interaction with
HeCiX, where HeCiX is answering questions which test both domain expertise and
knowledge about clinical trials.

Figure 3: Sample interaction with HeCiX demonstrating QA capabilities

###  6.3 Limitations

While HeCiX has shown promising results, it is important to acknowledge its
limitations.

  * •

Uncertainty of the model’s performance as the knowledge graph expands.

  * •

Additional testing on a wider range of diseases to ensure robustness of the
system.

##  7 Conclusion

This paper introduces HeCiX, an innovative system that connects knowledge
graphs from clinical trials data and Hetionet with large language models to
address major challenges in clinical research. Our experiments showcase that
HeCiX enhances drug discovery processes and uses existing scattered biomedical
data effectively. HeCiX uncovers all possible relationships between diseases,
genes, and treatments, potentially accelerating drug development which can
lead to unexpected discoveries.

HeCiX is a major advancement in using AI for clinical research. We believe
that HeCiX will play a crucial role in shaping the future of biomedical
research, providing a connected and innovative ecosystem in the field.

The scope of future enhancements for HeCiX could include adding SNOMED CT to
the knowledge graph for better clinical term standardization and employing
meta-path matching techniques to identify complex relationships more
effectively.

## References

  * [1] Sun D, Gao W, Hu H, Zhou S. Why 90% of clinical drug development fails and how to improve it? Acta Pharm Sin B. 2022 Jul;12(7):3049-3062. doi: 10.1016/j.apsb.2022.02.002. Epub 2022 Feb 11. PMID: 35865092; PMCID: PMC9293739.
  * [2] Himmelstein, Daniel Scott, et al. "Systematic integration of biomedical knowledge prioritizes drugs for repurposing." Elife 6 (2017): e26726.
  * [3] U.S. National Library of Medicine, “ClinicalTrials.gov,” https://www.clinicaltrials.gov/, accessed July 14, 2024.
  * [4] OpenAI, “GPT-4 Model Documentation,” https://openai.com/research/gpt-4/, accessed July 14, 2024.
  * [5] Harrison Chase, “LangChain,” https://github.com/langchain-ai/langchain, accessed July 14, 2024.
  * [6] Belleau, François, et al. "Bio2RDF: towards a mashup to build bioinformatics knowledge systems." Journal of biomedical informatics 41.5 (2008): 706-716.
  * [7] Chen, Ziqi, et al. "Ctkg: A knowledge graph for clinical trials." medRxiv (2021): 2021-11.
  * [8] Luo, Renqian, et al. "BioGPT: generative pre-trained transformer for biomedical text generation and mining." Briefings in bioinformatics 23.6 (2022): bbac409.
  * [9] Jin, Qiao, et al. "Matching patients to clinical trials with large language models." ArXiv (2023).
  * [10] Yuan, Hongyi, et al. "BioBART: Pretraining and evaluation of a biomedical generative language model." arXiv preprint arXiv:2204.03905 (2022).
  * [11] Labrak, Yanis, et al. "Biomistral: A collection of open-source pretrained large language models for medical domains." arXiv preprint arXiv:2402.10373 (2024).
  * [12] Devarakonda, Murthy V., et al. "Clinical trial recommendations using Semantics-Based inductive inference and knowledge graph embeddings." Journal of biomedical informatics 154 (2024): 104627.
  * [13] Neo Technology, “Neo4j AuraDB Documentation,” https://neo4j.com/docs/aura/current/, accessed July 14, 2024.
  * [14] Es, Shahul, et al. "Ragas: Automated evaluation of retrieval augmented generation." arXiv preprint arXiv:2309.15217 (2023).
  * [15] Anthropic, “Claude 3 Sonnet,” https://www.anthropic.com/news/claude-3-5-sonnet, accessed July 14, 2024.