# Emergence of Structure in Columns of Grains and Elastic Loops

Arman Guerra and Douglas P. Holmes  Mechanical Engineering, Boston University,
Boston, MA, 02215, USA

###### Abstract

It is possible to build free-standing, load-bearing structures using only
rocks and loops of elastic material. We investigate how these structures
emerge, and find that the necessary maximum loop spacing (the critical
spacing) is a function of the frictional properties of the grains and the
elasticity of the confining material. We derive a model to understand both of
these relationships, which depends on a simplification of the behavior of the
grains at the edge of a structure. We find that higher friction leads to
larger stable grain-grain and grain-loop contact angles resulting in a simple
function for the frictional critical spacing, which depends linearly on
friction to first order. On the other hand, a higher bending rigidity enables
the loops to better contain the hydrostatic pressure of the grains, which we
understand using a hydroelastic scale. These findings will illuminate the
stabilization of dirt by plant roots, and potentially enable the construction
of simple adhesion-less structures using only granular material and fiber.

Ensembles of dry, adhesion-less grains can flow like a fluid Lun et al.
(1984); Silbert et al. (2001); Kamrin and Koval (2012), or jam into a solid-
like state Majmudar et al. (2007); Song et al. (2008). There are two notable
mechanisms by which grains jam. The first is the densification induced by an
external force Liu and Nagel (1998); van Hecke (2009), e.g. coffee beans will
jam when vacuum packed in an air-tight bag. Jamming occurs when the number of
inter-grain contacts reaches a critical value – the aggregate becomes
isostatic, since each grain is held in place by its neighbors or the container
van Hecke (2009). As anyone who has run their hand through the sand on a beach
knows, however, these states are fragile, and a small change in the direction
of the external forces can often induce flow. The second mechanism is the
entanglement or interlocking of particles. This occurs either when the
particles have a high aspect ratio and are flexible, such as when birds build
a nest made from sticks Hansell (2000); Weiner et al. (2020), fibers are spun
into felt or cotton balls Kabla and Mahadevan (2007); Picu (2011); Gravish et
al. (2012), or when particles geometrically interlock, such as when ants
assemble themselves into bridges and rafts Anderson et al. (2002), or staples
are mixed together Franklin (2012); Murphy et al. (2016). These jammed states
can be very stable in the absence of any additional, external confinement.

Both of these mechanisms appear when plant roots penetrate into soil, i.e. an
elastogranular interaction Schunter Jr et al. (2018) between a slender elastic
structure and a granular material. The elastica will encircle and confine
grains, and at the same time entangle with themselves, securing themselves and
the grains that they contact in place Reubens et al. (2007); Ghestem et al.
(2014). Entangled plant roots aid in preventing landslides and stabilizing the
banks of rivers Gray (1981); Nilaweera and Nutalaya (1999). Recently this
method of stabilizing granular matter has been used to build structures –
researchers have combined different kinds of fibers and grains and shaped them
into load bearing walls and columns Aejmelaeus-Lindström et al. (2016, 2017);
Cohen et al. (2020).

Figure 1:  (a) The emergence of structure of a column of grains with exterior
loops of string. Each frame is a column prepared with a different initial
spacing $S_{i}$. The frame indicated with the blue curly bracket is prepared
with $S_{i}$ just below the critical spacing $S_{c}$. (b) The column
underlined in blue in (a) loaded to over $10^{5}$ Pa

In this Letter, we consider a simplified version of these structures, columns
made from grains and loops of elastic rods that provide external confinement
(Figure 1). These loops act as containers for the grains, but in contrast to
many other structures made from confined grains Wu and Hong (2008); Ramli et
al. (2013), the spacings between the loops can be larger than the grain
diameter. Therefore, we ask the question: what are the minimum criteria to
form a stable, elastogranular column?

##  I Methods

We will parameterize the positions of the loops as the spacing between them,
$S$. In the limit of very large initial loop spacing
($S_{i}\rightarrow\infty$), a column will collapse into a pile, the shape of
which is determined by the properties of the grains (Figure 1a, left) Al-
Hashemi and Al-Amoudi (2018). As $S_{i}\rightarrow 0$, the column will retain
most of its initial shape. We define the critical spacing $S_{c}$ as the
maximum loop spacing that will allow the columns to stand up with minimal
reduction in their height (schematic in Figure 2a) i.e. the final height
$H_{f}$ is within $95\%$ of $H_{i}$ (to account for grain settling, more
detail in Appendix A). When $S_{i}$ crosses below $S_{c}$, the column may
still lose rattler particles, i.e. particles which do not contribute to the
stability of the jammed state Baule et al. (2018), but retain its initial
imposed shape. Decreasing $S_{i}$ further leads to no qualitative change in
the final shape (Figure 1a, right). We note that a column prepared at
$S_{i}\approx S_{c}$ can handle a uniaxial compressive stress of over $10^{5}$
Pa (Figure 1b).

###  I.1 Experiments

To determine $S_{c}$, we built elastogranular columns which varied in $S_{i}$,
and considered a wide variety of granular matter (glass beads, plastic sous-
vide balls, peanut M&M’s, and ceramic rocks – diameters ranging between 1.0
and 1.9 cm) which we selected based on their varying frictional properties. We
note that all of the granular material we considered is approximately
spherical making it otherwise challenging to stabilize, i.e we do not consider
any grains with large aspect ratios Weiner et al. (2020), interlocking parts
Murphy et al. (2016) etc., which may otherwise become kinematically trapped
into a structural form. First, we poured the granular particles into an 8 cm
diameter hollow, hard, cylindrical slip-cast mold (Clear Cast Acrylic Tube,
3-1/2” OD x 3-1/8” ID, McMaster) until they reached a prescribed initial
height $H_{i}$ of 16cm. At regular intervals in this pouring process we
leveled the grains and placed an 8cm diameter loop around the exterior of the
grains. For our first experiments we used string (type 18 Twisted Mason’s
Twine, McMaster). We then removed the slip cast mold vertically in a quasi-
static manner (20mm/s) using a linear actuator (Zaber Technologies T-LSR300B),
and measured the final height $H_{f}$ of the elastogranular column when it
came to rest. We use the angle of repose $\alpha$, commonly defined as the
angle that a quasi-statically heaped pile of grains makes with the ground Al-
Hashemi and Al-Amoudi (2018), to account for the friction between grains, as
well as any slightly non-spherical geometrical features. In the case that the
grains are made from an ideal, cohesion-less Coulomb material, this angle is
related to the coefficient of friction by $\alpha=\phi=\arctan{\mu}$, where
$\phi$ is the angle of internal friction and $\mu$ is the coefficient of
static sliding friction Nedderman (2005). This assumption is not perfect – in
reality $\alpha$ is a complicated function of the rolling and sliding
friction, gradation, shape, etc. of the particles Zhou et al. (2002); Robinson
and Friedman (2002); Al-Hashemi and Al-Amoudi (2018). There are in fact many
definitions and ways to measure the shear strength and frictional properties
of granular materials, which may be applicable to different particle sizes,
shapes, gradations, and loading scenarios Al-Hashemi and Al-Amoudi (2018);
Chakraborty and Salgado (2010). However, since our grains have a comparably
high sphericity and are of uniform size and shape, and further since they are
under a self-load on a flat surface, we believe that the angle of repose is a
sufficient and appropriate metric to estimate the role of friction and
geometry in this work.

Intuition would suggest that if the grains are approximately spherical,
$S_{c}$ would be on the order of the grain diameter $d$, and that $S_{c}$ will
increase with the grain–grain friction. However, the strings are flexible, in
addition to being able to translate and rotate, and as such we find that if
the friction between the grains is low, the grains tend to push the string out
of the way and escape even when $S_{i}<d$ (Figure 2b, (i)). In the limit of
large $\alpha$, we find that grains will tend to collocate into stable
arrangements between strings (Figure 2b (iii)), allowing $S_{c}$ to reach as
high as $1.5d$. We plot $S_{c}/d$ vs $\alpha$ for our experiments with string
in Figure 2c (circles).

Figure 2:  (a) A schematic of our experiment. We build columns inside of an
acrylic slip-cast mold to a height $H_{i}$, then quasi-statically remove the
mold and measure $H_{f}$. The critical spacing $S_{c}$ is defined as the
minimum spacing such that the column does not reduce in height (rigorously
defined in Appendix A). (b) Edges of columns made from (i) Marbles
($\alpha=0$) and string, (ii) Ceramic ($\alpha\sim 30$) and acrylic rings, and
(iii) Ceramic and string. White lines are added to draw the eye to the loops
(c) $S_{c}$ vs $\alpha$ for experiments (strings– circles, acrylic – squares)
and simulations (diamonds). Data is colored by the bending rigidity of the
loops ($\sim Eh^{3}$) divided by the density of the grains $B/\rho$ (Pa
$\cdot$ m6/kg). The yellow data point has $B/\rho\sim 30\times 10^{-4}$ but we
have cut off the colorbar for visualization.

For the same range of granular material properties, we also investigated the
opposite limit – where the strings are replaced with rigid acrylic rings of
the same diameter and thickness (Figure 2c, squares). We found that $S_{c}$
for these experiments followed a similar trend in $\alpha$, and was always
higher than experiments using strings. We attribute the difference between the
behavior of columns made with string and with acrylic to the negligible
bending rigidity of the string, which is noted in Figure 2b (ii) – the grains
do not bend acrylic rings much, and as such the grains can be stable even when
they only overlap with the rings slightly. On the other hand when the loops
are flexible, the grains tend to bend them out of the way, thereby providing a
means for escape. One might expect that if the grains are perfectly
frictionless they may slip between the flexible strings no matter how small
the spacing is, whereas no rigid grain can move between acrylic loops which
are spaced less than the grain diameter.

###  I.2 Numerical Simulations

We next investigated the influence of bending rigidity of the confining loops
on the stability of the columns. Experimentally it is difficult to vary the
bending rigidity of the confining loops without changing their material and
geometric properties. To complement our experimental data, we repeated the
experiments in the Large-scale Atomic/Molecular Massively Parallel Simulator
(LAMMPS) Plimpton (1995) using hard (Young’s Modulus
$E_{g}=10^{9}~{}\mathrm{Pa}$, Poissons ratio $\nu=0.4$, coefficient of
restitution $(e)=0.35$) spheres (diameter $d$=1.15 cm). We use tangential and
rolling stiffnesses $K_{t}=4E_{g}/2(2-\nu)(1+\nu)$, $K_{r}=0.1K_{t}$ and
rolling damping coefficient $\gamma_{r}=1$ Jiang et al. (2005); Luding (2008);
Horabik and Molenda (2016). We varied the sliding friction coefficient
$\mu_{s}$ from 0.1-0.55 and we took the rolling friction coefficient
$\mu_{r}=\mu_{s}$ Jiang et al. (2005); Luding (2008); Horabik and Molenda
(2016). To make confining loops we simulate many small spheres (Figure 4e)
with nearest-neighbor potentials

$U=Y(r-r_{0})^{2}+B(1+\cos(\theta))$ (1)

Which includes a stretching term, a harmonic function of the distance between
adjacent loop particles $r$ and the equilibrium distance $r_{0}$, and a
bending term, which is a function of the angle between groups of three
respective loop particles $\theta$. If we take the stretching modulus
$Y=E_{l}\pi h/8$ and the bending modulus $B=E_{l}\pi h^{3}/64$ where $E_{l}$
is the Young’s modulus of the loops and $h$ is the thickness of the confining
loop, we will recover the formula for the energy of a cylindrical elastic loop
(Appendix B). We fix the stretching rigidity and vary the bending rigidity of
the loops in simulations (Figure 2c, diamonds). We will note that, in addition
to enabling us to systematically vary the bending rigidity of the loops, the
simulations were complementary to our experiments in that they allowed us to
study an experimentally inaccessible range of $\alpha$, specifically
$\alpha\in(1,10)\cup(16,26)$, for which were unable to find grains in our
target diameter range, and enabled a systematic study of the dependance of the
bending of the loops on each experimental parameter, thereby helping to
justify a scaling law for the system (see Figure 4).

##  II Analysis

These experiments and simulations so far indicate that the critical spacing at
which columns stand up is a function of the angle of repose $\alpha$ of the
grains and the bending rigidity $B$ of the confining loops. To understand how
these factors play a role, we will consider a reduced-order model of the grain
arrangement at the edge of a column (which takes into account the possible
contributions of grain-grain friction, grain-loop friction, and loop bending)
and find the conditions for local stability. This simplification will allow us
to establish some guidance on what combination of material and geometric
parameters will enable the emergence of a column that can bear its own weight.

Consider an arrangement of three grains near the edge of a column, one
interior grain contacting two exterior grains both of which contact the
confining loops (Figure 3a). We will assume that the grains are spherical and
monodisperse, and we note that since we have illustrated a 2 dimensional
lateral cross section of a 3 dimensional scenario, the exterior grains may
vertically overlap with one another. In the absence of the exterior grains,
the interior grain would escape the column radially, and as such there must be
an outward radial force acting on the interior grain which is balanced by the
exterior grains. Consequently, there must be a force acting inwards on the
exterior grains from the confining loops – to understand the role of grain-
grain friction we will first assume that this force results from a hoop
stress, that is, it acts radially inwards. We illustrate the outward radial
force as $F_{r}$ and the inward hoop force as $F_{h}$.

Figure 3:  (a) A schematic of our mathematical model. An interior grain is
subject to an outward radial force but held in place by two exterior grains
which contact the loops. According to the Coulomb constraint, for the grain-
grain contact point to be stable, the angle between the centers of the
interior and exterior grains can be at most $\alpha$, so the maximum stable
distance between loops $S_{c}=2d\sin{\alpha}+2y_{c}$. (b) We denote
$y_{c}^{f}$ as the maximum stable height of the loops with respect to the
grains due to frictional stability, which can be found with a similar
constraint as above. (c) (i) An elastic rod loaded with a force $P$ will bend
a distance $\widetilde{y}$ proportional to the ratio of the applied force and
the bending modulus. (ii) The radial force from the interior grains threatens
to bend the loops out of the way. (iii) We denote
$y_{c}^{B}=d/2-\widetilde{y}$ as the maximum stable height of the loops with
respect to the grains due to the bending rigidity of the loops.

From here, we can determine whether the contact point between the interior and
exterior grains will be stable or will slip based on their relative positions.
Using the familiar Coulomb constraint, we can say that the contact will be
stable if $\mu F_{h,n}\geq F_{h,t}$ where $F_{h,n}$ and $F_{h,t}$ are the
components of the hoop force that are normal and tangential, respectively, to
the tangent line of contact between the spheres, illustrated in Figure 3a.
This provides an upper bound for the angle $\theta$ between the horizontal and
the vector which points from the centers of the interior and exterior grains –
we can say
$\alpha=\arctan{\mu}\geq\arctan{F_{h,t}/F_{h,n}}=\arctan{F_{h}\sin{\theta}/F_{h}\cos{\theta}}=\theta$.
Therefore the maximum distance between the centers of the exterior grains is
$2d\sin{\alpha}$.

The only further unknown is the maximum vertical distance between the loops
and the centers of the exterior grains. To find this we will separate out the
influence of the friction from the influence of the bending rigidity – one
could imagine that if the loops had negligible bending rigidity, but there was
a high degree of friction between the grains and the loops, the loops may
constrain the grains because the contact point between the loops and the
grains does not slip. Alternatively, in the case of infinitely slippery grains
and loops, if the loops were stiff, the grain-loop contact point may slip but
the loops may not bend out of the way enough for the grains to escape. We will
call $y_{c}^{f}$ the maximum critical height of the loops due to frictional
stability, and $y_{c}^{B}$ the maximum critical height of the loops due to
bending rigidity. We will then find two versions of the largest stable
spacing, $S_{c}^{f}=2d\sin{\alpha}+2y_{c}^{f}$ and
$S_{c}^{B}=2d\sin{\alpha}+2y_{c}^{B}$. The actual critical spacing will be the
overall largest stable spacing, and therefore is the larger of the two –
$S_{c}/d=\mathrm{max}(S_{c}^{f}/d,S_{c}^{B}/d)$

Figure 4:  (a) We fix $\alpha\sim 17$deg and measure the change in
$\widetilde{y}$ as we individually vary the height of the columns, (b) radius
of the columns, (c) diameter of the particles and (d) thickness of the loops
in simulations. We normalize each input variable and $\widetilde{y}$ by the
values in the original simulations
($H_{0},R_{0},d_{0},h_{0},\widetilde{y}_{0}$). The black guide lines have the
slopes from our mathematical model: 1, 4, 0, and -1 respectively, and the blue
lines indicate the range of the parameters in experiment. (e) Two images from
our simulations, Stukowski (2009) and a 3 dimensional graph of a ring selected
from each. The z-axis is scaled up to emphasize the curvature of the ring. The
bending length $L$ of the rings is a function of both $R$ and $d$. We ignored
the $d$ contribution, which is why we overestimate the dependance of
$\widetilde{y}$ at high $R$ and underestimate at low $d$. (f) Comparison
between
$S_{c}^{theory}/d=\mathrm{max}(3\sin{\alpha},\frac{2\sin{\alpha}+1}{1+2\beta_{he}})$
and our experimental values of $S_{c}$. We color by the inverse of the
hydroelastic scale, $1/\beta_{he}$ and fill based on whether experiments are
in a friction controlled or bending controlled region of $(\alpha,\beta_{he})$
space. The yellow data point has $1/\beta_{he}\sim 150$ but we have cut off
the colorbar for visualization. Error bars are taken from Figure 2 and do not
include the error in $x$. In the inset we show the data from Figure 2c with
our analytical results plotted. The thicker line represents $S_{c}^{f}$ and
the thinner lines represent $S_{c}^{B}$ for different $\beta_{he}$. The dotted
sections represent when $S_{c}^{B}$ is below $S_{c}^{f}$.

We start with $y_{c}^{f}$. If we still consider the hoop force acting between
the exterior grains and the loops, and we assume that the grain-grain
coefficient of friction is the same as the grain-loop coefficient of friction,
we find ourselves with an almost identical constraint as before – the maximum
value of the angle between the horizontal and the vector pointing between the
centers of the grains and the loops is $\alpha$ and we find that
$y_{c}^{f}=\frac{d}{2}\sin{\alpha}$

Now, to find $y_{c}^{B}$ we must consider the force which acts to bend the
loop out of the way. We will set aside the grain-loop friction and relax the
hoop-force assumption, such that now there is a vertical force between the
grains and the loops. When a force $P$ acts on an Euler–Bernoulli beam with a
cylindrical cross-section $h$, length $L$, bending modulus $B$, and clamped
edges, it will bend with a maximum amplitude

$\widetilde{y}=\frac{PL^{3}}{192Bh}$ (2)

Illustrated in Figure 3c (i). If the loop is to constrain a grain, then it
must not bend out of the way so much as to clear the top of the grain (Figure
3c (ii)). We can therefore take the height of the top of the grain $d/2$ and
subtract the bending of the loop $\widetilde{y}$ to find the maximum height of
the loop such that it will not bend out of the way of the particle,
$y_{c}^{B}=d/2-\widetilde{y}$ illustrated in Figure 3c (iii). To put
$y_{c}^{B}$ in terms of the variables of our problem, we will assume that the
force on each loop is due to the hydrostatic pressure from the grains, which
will be $\phi_{rlp}\rho gD$ at a depth $D$, density $\rho$, and random loose
packing fraction $\phi_{rlp}\sim 0.55$. The area associated to each loop is
$2\pi RS_{c}$ where R is the radius of a column. Therefore the average force
on each loop is $P\sim\phi_{rlp}\rho g(H/2)2\pi RS_{c}$ where $H$ is the
height of a column. Depending on how the force is distributed on each loop,
the bending length L will either scale with the total length of the loop $\sim
R$ or on the fluctuations in the force applied to the loop $\sim d$. We find
qualitatively from our simulation data that, in the range of $d/R$ that we
study, most of the force on a given loop can be attributed to a small number
of grains, which implies that $L$ will be limited by $R$. This gives the
result

$y_{c}^{B}=\frac{d}{2}-\frac{\phi_{rlp}\pi}{192}\left(\frac{\rho
gHRS_{c}}{Bh/R^{3}}\right)$ (3)

Where we have expressed the second term as a ratio between the hydrostatic
force on the loops and their bending rigidity. We will separate $S_{c}$ from
this second term and wrap the rest of the term into a hydroelastic scale,
$\beta_{he}=\frac{\phi_{rlp}\pi}{192}\left(\frac{\rho gHR}{Bh/R^{3}}\right)$.

We note that in our experiments and simulations so far the only variables
which have changed significantly are $S_{c}$, $\rho$, and $B$. To further test
our scalings of $\widetilde{y}$ in the variables that we had not yet varied,
we performed some additional simulations, the results of which are shown in
Figure 4. We note that when $d/R$ is small, $\widetilde{y}$ starts to scale
with $d$ and the $R$ dependence reduces. This occurs because as the number of
contacts with the loop increases, the force fluctuations start to limit the
bending of the loops, and the bending wavelength shortens (Figure 4e). But,
within our experimental range (marked in blue), our scalings seem appropriate.

We can now plug in our values of $y_{c}^{f}$ and $y_{c}^{B}$ to find the
critical spacing. Rearranging, we find that

$\begin{split}&S_{c}^{f}/d=3\sin{\alpha}\\\
&S_{c}^{B}/d=\frac{2\sin{\alpha}+1}{1+2\beta_{he}}\end{split}$ (4)

The expected stability is therefore set by the maximum of these two values,
one which is solely dependent on the friction ($\alpha$) and a second that is
dependent on both friction and the bending rigidity of the loops ($\alpha$ and
$\beta_{he}$ where $\beta_{he}\propto B^{-1}$, see Equation 3). We plot
$S_{c}/d$ against our theoretical finding in Figure 4f and find good agreement
between our experiments and the behavior derived from our simplified model,
indicating that this model provides good intuition for the interactions at
play. This reduced order model describes the minimum conditions necessary for
a column of loops encased by a finite number of elastic loops to retain its
shape when placed in a gravitational field. However, the simplifications in
the model make it difficult to comment on the stability of the resulting
equilibrium shapes. Indeed, these columns may be fragile to eccentric loading,
shear, or bending, and may be sensitive to defects and imperfections which are
effectively averaged over by our inherent assumptions. These are important
considerations left to future work, however we reiterate that, once these
structures are jammed in response to uniaxial loading, they are capable of
bearing significant loads (Figure 1).

##  III Adaptable Forms

Figure 5:  (a) (i) Snapshots of an experiment with ceramic grains where
columns are built with an initial height of 33cm and $S_{i}\sim 3.6S_{c}$.
(ii) The final state of the columns for different runs of the experiment shown
in (i). As grains escape between the strings, the columns shrink until they
end up at a $H_{f}$ which depends on the number of spacings $N_{S}$ (which
includes the spacing between the bottom (top) of the column and the bottom
(top) string). (b) Left: Height of columns made from ceramic grains as the
slip cast-mold is removed. We normalize the height by the product
$S_{c}N_{S}$. If each spacing was exactly $S_{c}$ this value would be 1.
Right: We find $S_{f}$ at multiple points around the circumference of the
cylinder and plot the average (blue circles) and median (red diamonds) value,
normalized by $S_{c}$, for different $\alpha$. (c) We can use these ideas to
make structures that are stable in compression and can be dismantled by a
simple user input.

Now finally we ask the question: what happens when we prepare columns with
$S_{i}>S_{c}$? When constructed on flat ground, the columns become encompassed
by a pile of the grains which escape, as shown in Figure 1a. This pile
obscures the final state of the loops. If instead, we set the columns on a
pedestal that is the same diameter as the diameter of the columns, the
escaping grains fall away and we can observe the dynamics of the loops.

The evolution of the average spacing of strings $S$ for various columns with
$S_{i}>S_{c}$ is shown in Figure 5b. Remarkably, as the columns settle,
$S\rightarrow S_{c}$. In other words, a string will fall until the spacing
between it and the string below it (or the ground) is, on average, _equal to_
$S_{c}$, irrespective of the initial spacing. The result is a shorter column
with a final spacing $S_{f}=S_{c}$. A striking example of this phenomenon is
shown in Figure 5a. This behavior emphasizes the robustness of a $S_{c}$ which
is impervious to dynamical effects, and indicates the possibility that this
technology could be used to construct collapse-mitigating structures.

Another interesting feature of these aggregate structures is that there is no
adhesion between respective elements. This allows them to respond in
drastically different ways to slightly different inputs, and quickly and
easily change shape. One could imagine building stable elastogranular
structures that could handle large compressive loads, but then be shortened or
demolished in seconds with the appropriate input. One example of such a
structure is shown in Figure 5c. We 3-d printed a plastic (ABS P430) helix
which we filled with rocks to form a column. This column was stable in
compression, however when the helix was uncoiled, which could easily be done
by hand, it collapsed layer-by-layer, allowing for some control over its final
height.

##  IV Conclusions

In this work, we set out to understand how elastic rods constrain and jam
granular matter. We have found that the elastogranular interaction in this
regime is predictable and robust, allowing us to quickly and reliably form
structure which can bear significant load (although the limits of this load,
and the intricacies of the mechanical response of the columns to compression
is beyond our current scope). The stability of these structures depends on the
size and frictional properties of the grains, as well as flexibility of the
rods. We have so far only considered external loops of elastic, although it
has been shown that internal elastic rods are also sufficient to form
structure from granular material Aejmelaeus-Lindström et al. (2017). We expect
that this will be related to the mechanics which govern the knotting of ropes
Bayman (1977); Maddocks and Keller (1987) and, similar to the jamming of
chains of beads Brown et al. (2012); Dumont et al. (2018) could be understood
through a similar mathematical framework as governs the entanglement of
polymers Doi (1996).

## Conflicts of interest

There are no conflicts to declare.

## Appendix A: critical spacing

To find the critical spacing for both experiments and simulation, we build
columns with a fixed initial height $H_{i}$ and vary the number of strings
with equal spacings $S_{i}$ between them. We linearly fit the data for which
$H_{f}<H_{i}$ and extrapolate to find the point where $H_{f}=\gamma H_{i}$. We
use $\gamma=0.95$ for experiments and $\gamma=0.9$ for simulations.

The reason that we use a tolerance factor $\gamma$ which is less than 1 is
because we found that, in practice, when we removed the slip-cast mold, the
grains in the columns had a tendency to settle, and even if $S_{i}<<S_{c}$
some of the grains from the top of the columns would dislodge, shortening the
column slightly. We chose $\gamma$ empirically as the fraction of $H_{i}$ that
most columns tended to surpass when $S_{i}<<S_{c}$.

Figure 6:  To find the critical spacing we plot the final height of the
columns vs the number of strings, which is linear when $S_{i}>S_{c}$ (filled
circles) based on the data presented in Figure 5. We fit this data and find
$S_{c}$ as the spacing at which $H_{f}>\gamma H_{i}$ where $\gamma$ is an
empirical tolerance factor close to 1. The error on the linear fit gives us
the error in $S_{c}$.

## Appendix B: simulation of loops

We would like to replicate the physics of a thin cylindrical loop with Young’s
modulus $E$, circumference $l$, and cross sectional diameter $h$ (Figure 7a).
We will connect particles of diameter $h$ in a loop (Figure 7b) with the
following potentials.

To add a stretching rigidity to the loop we will apply a harmonic “Bond”
between adjacent particles in the beam using the LAMMPS `bond_style harmonic`
pair potential. We define an equilibrium distance between adjacent particles,
or equilibrium Bond length $r_{0}$ and calculate the energy of a Bond between
particle $i$ and particle $j$ which has a Bond length $r_{ij}$ as

$U_{s}=Y(r_{ij}-r_{0})^{2}=Y(\Delta r_{ij})^{2}$ (5)

Where we define $\Delta r_{ij}$ to be the change from equilibrium of the Bond
length. The force on the particles are as such $F=dU_{s}/dq=2Y\Delta r_{ij}$.
In all of the simulations above, $r_{0}$ is equal to the diameters of the
particles $h$, that is, the potential acts to keep the particles of the beam
just in contact. This is not necessary, one could imagine a beam made up of
more spherical particles with some neighbor overlap, or fewer spherical
particles with a larger distance between each pair.

We take the definition of the Young’s modulus $E=\sigma/\epsilon=(F/A)/(\Delta
l/l_{0})$ where A is the cross sectional area of the loops $\pi(h/2)^{2}$. If
we statically compress or stretch a beam, the force on the ends will be equal
to the force between any adjacent particle, which we found above to be
$F=2Y\Delta r$ (we have dropped the subscript of $r$ because the Bond lengths
will all be the same). The strain of the whole loop $\Delta l/l_{0}$ will be
equal to the strain of each Bond $\Delta r/r_{0}$ so we have

$E=\frac{F/A}{\Delta l/l_{0}}=\frac{2Y\Delta r}{\pi(h/2)^{2}}\frac{h}{\Delta
r}=Y\frac{8}{\pi h}\implies Y=\frac{\pi}{8}Eh$ (6)

Where we have used the fact that $r_{0}=h$.

To derive the pair potential needed to induce the correct bending rigidity in
our LAMMPS beams we will consider first a continuous beam (Figure 7c inset,
left). We will bend the beam and find the resultant bending energy, and we
will use that as the target energy for a discrete beam bent in the same
orientation.

For a continuous elastic beam with a moment $I$, the energy due to bending is

$U_{b}^{c}=\frac{EI}{2}\int\kappa^{2}ds$ (7)

Where $\kappa$ is the curvature of the beam, and we have used the superscript
$c$ to indicate that the beam is continuous. To achieve an energy of bending
in LAMMPS we will use what is called an “Angle.” A LAMMPS Angle is like a Bond
(which we used in the previous section on the stretching potential) except
instead of considering the interaction between two particles, each Angle
applies a potential based on the relative positions of three particles. The
specific Angle potential that we will use is `angle_style cosine` which
applies a potential

$U_{b}^{d}=B(1+\cos\theta_{ijk})$ (8)

Where $\theta_{ijk}$ is the angle between the three particles (Figure 7c
inset, right) and we have used the superscript $d$ to indicate that the beam
is discrete. We will bend both the continuous and the discrete beams to a
radius of curvature $R$ such that for both beams $\kappa=1/R$ everywhere
(Figure 7c, inset). The energy in a chunk of the continuous beam of width $h$
will be

$U_{b}^{c}=\frac{EI}{2R^{2}}\int ds=\frac{EI}{2R^{2}}h=\frac{E\pi
h^{5}}{128R^{2}}$ (9)

Figure 7:  (a) Theoretical predictions for the shape of a loop of
circumference $l$ loaded with a force $F$ at two points opposite one another.
(b) Results of a LAMMPS simulation of the setup in (a). (c) Non-
dimensionalized Force $Fl^{2}/EI$ vs displacement $\Delta/l$ for a theoretical
Euler-Bernoulli loop and a loop simulated in LAMMPS. In the inset we
illustrate two beams bent with a radius of curvature R – (left) a continuous
cylindrical beam of diameter h, and (right) a discrete cylindrical beam made
up of particles of diameter h. The energy of bending of the discrete beam will
depend on the angle $\theta$ between each triplet of particles.

Where we have taken the moment of the cylindrical beam $I=\pi h^{4}/64$. In
the discrete case, in a chunk of width $h$ in the bulk of the beam there is
one particle (since the particles have width $h$). The energy assigned to that
particle because of the bending of the beam can be found by summing up the
contributions of each of the Angles that it is a part of. There are three
particles in each Angle, so to each constituent particle we will assign one
third of the energy in that Angle. Furthermore each particle is a part of
three angles, so the energy in each particle in the column due to bending, and
therefore the energy in a chunk of the beam of width $h$ is
$3*(1/3)B(1+\cos\theta)=B(1+\cos\theta)$ where we have used the fact that the
entire beam is bent to the same curvature so each angle $\theta_{ijk}$ between
all sets of particles is the same.

The only thing that is left is to connect $(1+\cos\theta)$ to the radius of
curvature $R$. We will do this geometrically referring to the angles and
lengths defined in Figure 7c inset, right. By the law of cosines, we have that

$h^{2}=2R^{2}-2R^{2}\cos\phi=2R^{2}(1-\cos\phi)$ (10)

Since the interior angles of a triangle must sum to 180 degrees, we have that
$\phi=180+\theta/2+\theta/2\implies\theta=180-\phi$ so
$1+\cos\theta=\frac{h^{2}}{2R^{2}}$. Setting the energies equal we have

$U_{b}=\frac{Bh^{2}}{2R^{2}}=\frac{E\pi h^{5}}{128R^{2}}\implies B=\frac{E\pi
h^{3}}{64}$ (11)

Thus the total energy of the loops is

$U=\frac{\pi}{8}Eh\sum(r_{ij}-r_{0})^{2}+\frac{E\pi
h^{3}}{64}\sum(1+\cos(\theta_{ijk}))$ (12)

To validate this result we simulate a loop compressed by two point forces and
plot the distance between the points $\Delta$ as the force on the loop
increases (Figure 7). We find good agreement between the simulation and the
theoretical prediction.

## Acknowledgements

The authors gratefully acknowledge the financial support from DARPA
(#HR00111810004) and from NSF CMMI–CAREER through Mechanics of Materials and
Structures (#1454153), and the computing resources of the Boston University
Shared Computing Cluster. We also thank Kate Flanagan and Xin Jiang for
initial experimental design and Skylar Tibbits for initial discussions.

## References

  * Lun et al. (1984) C. Lun, S. B. Savage, D. Jeffrey, and N. Chepurniy, Journal of fluid mechanics 140, 223 (1984).
  * Silbert et al. (2001) L. E. Silbert, D. Ertaş, G. S. Grest, T. C. Halsey, D. Levine, and S. J. Plimpton, Physical Review E 64, 051302 (2001).
  * Kamrin and Koval (2012) K. Kamrin and G. Koval, Physical Review Letters 108, 178301 (2012).
  * Majmudar et al. (2007) T. Majmudar, M. Sperl, S. Luding, and R. P. Behringer, Physical review letters 98, 058001 (2007).
  * Song et al. (2008) C. Song, P. Wang, and H. A. Makse, Nature 453, 629 (2008).
  * Liu and Nagel (1998) A. J. Liu and S. R. Nagel, Nature 396, 21 (1998).
  * van Hecke (2009) M. van Hecke, Journal of Physics: Condensed Matter 22, 033101 (2009).
  * Hansell (2000) M. Hansell, _Bird nests and construction behaviour_ (Cambridge University Press, 2000).
  * Weiner et al. (2020) N. Weiner, Y. Bhosale, M. Gazzola, and H. King, Journal of Applied Physics 127, 050902 (2020).
  * Kabla and Mahadevan (2007) A. Kabla and L. Mahadevan, Journal of the Royal Society Interface 4, 99 (2007).
  * Picu (2011) R. Picu, Soft Matter 7, 6768 (2011).
  * Gravish et al. (2012) N. Gravish, S. V. Franklin, D. L. Hu, and D. I. Goldman, Physical review letters 108, 208001 (2012).
  * Anderson et al. (2002) C. Anderson, G. Theraulaz, and J.-L. Deneubourg, Insectes sociaux 49, 99 (2002).
  * Franklin (2012) S. V. Franklin, Physics Today 65, 70 (2012).
  * Murphy et al. (2016) K. A. Murphy, N. Reiser, D. Choksy, C. E. Singer, and H. M. Jaeger, Granular Matter 18, 26 (2016).
  * Schunter Jr et al. (2018) D. J. Schunter Jr, M. Brandenbourger, S. Perriseau, and D. P. Holmes, Physical review letters 120, 078002 (2018).
  * Reubens et al. (2007) B. Reubens, J. Poesen, F. Danjon, G. Geudens, and B. Muys, Trees 21, 385 (2007).
  * Ghestem et al. (2014) M. Ghestem, G. Veylon, A. Bernard, Q. Vanel, and A. Stokes, Plant and Soil 377, 43 (2014).
  * Gray (1981) D. H. Gray, _Forest vegetation removal and slope stability in the Idaho Batholith_ , vol. 271 (US Department of Agriculture, Forest Service, Intermountain Forest and Range …, 1981).
  * Nilaweera and Nutalaya (1999) N. Nilaweera and P. Nutalaya, Bulletin of engineering geology and the environment 57, 337 (1999).
  * Aejmelaeus-Lindström et al. (2016) P. Aejmelaeus-Lindström, J. Willmann, S. Tibbits, F. Gramazio, and M. Kohler, Granular Matter 18, 28 (2016).
  * Aejmelaeus-Lindström et al. (2017) P. Aejmelaeus-Lindström, A. Mirjan, F. Gramazio, M. Kohler, S. Kernizan, B. Sparrman, J. Laucks, and S. Tibbits, Architectural Design 87, 82 (2017).
  * Cohen et al. (2020) Z. Cohen, N. Elberfeld, A. Moorman, J. Laucks, S. Kernizan, H. D. P., and S. Tibbits, Technology — Architecture + Design 4:2, 211 (2020).
  * Wu and Hong (2008) C.-S. Wu and Y.-S. Hong, Geotextiles and Geomembranes 26, 302 (2008).
  * Ramli et al. (2013) M. Ramli, T. Karasu, and E. T. Dawood, Alexandria Engineering Journal 52, 705 (2013).
  * Al-Hashemi and Al-Amoudi (2018) H. M. B. Al-Hashemi and O. S. B. Al-Amoudi, Powder Technology 330, 397 (2018).
  * Baule et al. (2018) A. Baule, F. Morone, H. J. Herrmann, and H. A. Makse, Reviews of Modern Physics 90, 015006 (2018).
  * Nedderman (2005) R. Nedderman, _Statics and Kinematics of Granular Materials_ (Cambridge University Press, 2005), ISBN 9780521019071, URL https://books.google.com/books?id=52XuyGlewh8C.
  * Zhou et al. (2002) Y. Zhou, B. H. Xu, A.-B. Yu, and P. Zulli, Powder technology 125, 45 (2002).
  * Robinson and Friedman (2002) D. Robinson and S. Friedman, Physica A: Statistical Mechanics and its Applications 311, 97 (2002).
  * Chakraborty and Salgado (2010) T. Chakraborty and R. Salgado, Journal of geotechnical and geoenvironmental engineering 136, 527 (2010).
  * Plimpton (1995) S. Plimpton, Journal of computational physics 117, 1 (1995).
  * Jiang et al. (2005) M. Jiang, H.-S. Yu, and D. Harris, Computers and Geotechnics 32, 340 (2005).
  * Luding (2008) S. Luding, Granular matter 10, 235 (2008).
  * Horabik and Molenda (2016) J. Horabik and M. Molenda, Biosystems Engineering 147, 206 (2016).
  * Stukowski (2009) A. Stukowski, Modelling and Simulation in Materials Science and Engineering 18, 015012 (2009).
  * Bayman (1977) B. F. Bayman, American Journal of Physics 45, 185 (1977).
  * Maddocks and Keller (1987) J. H. Maddocks and J. B. Keller, SIAM Journal on Applied Mathematics 47, 1185 (1987).
  * Brown et al. (2012) E. Brown, A. Nasto, A. G. Athanassiadis, and H. M. Jaeger, Physical review letters 108, 108302 (2012).
  * Dumont et al. (2018) D. Dumont, M. Houze, P. Rambach, T. Salez, S. Patinet, and P. Damman, Physical review letters 120, 088001 (2018).
  * Doi (1996) M. Doi, _Introduction to polymer physics_ (Oxford university press, 1996).


# Could Geoneutrinos Interact With the Geomagnetic Field?

C. A. B. Quintero1 and J.A. Helayël Neto2

1 _Departamento de Geofísica, Observatório Nacional (ON),_
_Bairro Imperial de São Cristovão, Rio de Janeiro, RJ, 20921-400, Brazil_
2 _Centro Brasileiro de Pesquisas Físicas (CBPF),_
_Urca, Rio de Janeiro, RJ, 22290-180, Brazil_
∗E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In the present paper, we consider the possibility of interaction between
geoneutrinos and the geomagnetic field, by adopting an approach based on the
Dirac’s equation with a non-minimal coupling that accounts for the magnetic
interaction of the massive neutrinos. In our approach, we see that the
magnetic interaction is controlled by a dimensionless parameter, $f\simeq
10^{-1}$, and we estimate the mean value of this interaction to be of the
order of $10^{-14}\ MeV^{2}$.

##  1 Introduction

Heat is the engine of the dynamics of the planet and its magnetic field, so
that it becomes impossible to split these two quantities [1] [2]. At present,
Earth is releasing heat from its surface at a rate of about $47\ TW$ [3] [4].
The principal contribution to the Earth loss are the secular cooling of the
Earth and the decay of long-lived radioactive isotopes of uranium, thorium and
potassium.
The Bulk Silicate Earth (BSE) model establishes that the Earth’s chemical
composition must be the same as the one of meteorite chemical composition from
primordial cloud. The chondritic ratio, $m_{Th}/m_{U}$, is one of the main
characteristic of primordial cloud chemical composition and varies from $2.6$
to $4.2$ [5] [6]. BSE model gives values of mass for U, Th and K such that
$m_{U}=0.81\times 10^{17}\ kg$, $m_{Th}=3.16\times 10^{17}\ kg$ and
$m_{K}=0.49\times 10^{21}\ kg$. These amounts are distributed only in the
Crust and Upper Mantel, this is, the BSE establishes that U and Th are not
present in the core, although alternative hypothesis have been studied [7]
[8]. The main relevance of the chondritic ratio is that it is directly related
to the flux of electronic antineutrinos, or geoneutrinos, from beta-decay of
long-lived radioactive isotopes.
Geoneutrinos are detected via inverse beta-decay,
$\bar{\nu_{e}}+p\longrightarrow e^{+}+n$, with an energy threshold of $1.806\
MeV$. The prompt scintillation light from $e^{+}$ gives a measure of the
$\bar{\nu_{e}}$-energy, $E_{\bar{\nu_{e}}}\simeq E_{p}+\bar{E_{n}}+0.8MeV$,
where $E_{p}$ is the prompt event energy including the positron kinetic and
annihilation energies, and $\bar{E}$ is the average neutron recoil energy,
$\sigma\left(10KeV\right)$ [9]. Actually, neutrinos are produced only in
electron capture of ${}^{40}K$. In contrast to the Sun, Earth shines
essentially in antineutrinos, ie, in geoneutrinos.
(Anti)neutrinos in the Standard Model (SM) do not have mass. As it is well-
known, a massless, chiral neutrino cannot have a nonzero magnetic (and
electric) dipole moment. However, beyond the SM, a massive Dirac or Majorana
neutrino will, in general, exhibit magnetic moment [10]. Therefore, the
geoneutrino feels magnetic fields. Low [11] [12] has studied the interaction
between particles through non-minimal couplings, opening up the possibility of
understanding another kind of interaction, in this case, the spin of the
neutrino with magnetic fields.
The purpose of this work is to study the possibility of an interaction of the
massive geoneutrinos with the Geomagnetic field. For this purpose, we adopt a
model from theoretical physics that is a relativistic generalization of the
conventional spin-current, using non-minimal couplings. However, here, in our
approach, we shall focus the discussion on the non-relativistic limit and then
we will adapt our theoretical description for the case of geoneutrinos.
This work is organized as follows. In Section 2, we describe the formalism and
calculate the Lorentz force for geoneutrinos. In Section 3, we analyze the
propagation of geoneutrinos through the matter and we derive our approximateve
expressions for the interaction between geoneutrinos and Geomagnetic field.
Finally,our Final Considerations and Future Prospects are cast in Section 3.

##  2 Presentation of the Model

The interaction between the geoneutrinos and matter is described by the weak
nuclear force. geoneutrinos are electronic antineutrinos, hence, they are
fermions and their energy values are known. We are going to consider them as
relativistic particles that obey the Dirac equation for spin-$\frac{1}{2}$
relativistic particles. In this approach, we are taking neutrinos as massive
particles. Here, we present, for the sake of completeness, the Dirac’s
equation for a Dirac or Majorana fermion with non-minimal couplings in
addition. However, for our practical purposes in this paper, since we are
concerned with neutrinos, we are going to finally set $e=0$, and focus on the
parameters that govern the non-minimal couplings. Our investigation is not
committed with the neutrino being a Majorana or a Dirac particle. In view of
that, the covariant Dirac equation with non-minimal couplings is given by,

$i\gamma^{\mu}\partial_{\mu}\psi-m\psi-eA_{\mu}\gamma^{\mu}\psi-
fZ_{\mu}\gamma^{\mu}\gamma_{5}\psi+igZ_{\mu\nu}\Sigma^{\mu\nu}\gamma_{5}\psi=0.$
(1)

Through this paper, we shall to work in natural unitis ($\hbar=c=1$). Wherever
necessary to recover the correct dimensions, we re-insert $\hbar$ and $c$ in
the right places. The non-minimal couplings $f$ and $g$ are both nontrivial
for either type of neutrino, once the neutrino field is described by Grassman
value spinors. From the equation above, we get a Hamiltonian in which there
appear new terms that correspond to a non-minimal interactions,

$\displaystyle H\psi=$
$\displaystyle\alpha^{i}\left[-i\partial_{i}-eA^{i}-fZ^{i}\gamma_{5}\right]\psi+m\beta\psi+e\phi\psi$
(2)
$\displaystyle+fZ_{0}\gamma_{5}\psi+ga^{i}\gamma^{i}\gamma_{5}\psi+2cgb^{i}\beta\gamma_{5}\Sigma^{i}\psi.$

Before going on, it is necessary to set up some definitions for a better
understanding of the parameters of the system. The field
$\mathbf{Z_{\mu\nu}}$, present in the Hamiltonian, is a background field that
can be split into $SO\left(3\right)$ irreducible representations. We define
them as $\vec{a}$ such that $a_{i}$, and $\vec{b}$ such that $b_{k}$, thus:

$\displaystyle a_{i}$ $\displaystyle=$
$\displaystyle-\partial_{t}Z_{i}-\partial_{i}Z_{0},$ $\displaystyle b_{k}$
$\displaystyle=$ $\displaystyle\epsilon_{ijk}\partial_{i}Z_{j}.$ (3)

It is important to point out that $\vec{b}$ plays the rôle of a kind of
background magnetic field that we decide to interpret as the geomagnetic
field.

###  2.1 The Lorentz-like Force

We now consider the situation of geoneutrinos in presence of the background
fields, $\vec{a}$ and $\vec{b}$, that we shall relate to the wave geomagnetic
field. After some algebraic manipulation, we arrive at the Lorentz force in
terms of the external fields for the Hamiltonian of Eq(1),

$\displaystyle\mathfrak{F}^{i}$
$\displaystyle=f\gamma_{5}\left(\partial_{i}Z_{0}-\partial_{t}Z^{i}\right)-f\Sigma^{j}\left(\partial_{j}Z^{i}-\partial_{j}Z^{i}\right)$
(4)
$\displaystyle+g\partial_{i}a^{j}\gamma^{j}\gamma_{5}\psi-2g\partial_{i}b^{j}\gamma^{j}-2ifmZ^{i}\gamma_{5}\beta$
$\displaystyle+4icgfZ^{i}b^{j}\beta\Sigma^{j}+2igfZ^{i}a^{j}\gamma^{j},$

where we do not see the conventional Lorentz force for the neutrinos
chargeless (we have set $e=0$ in our Dirac’s equation). On the other hand, we
can see that there appears in the Lorentz force a new kind of interaction
associated to the non-minimal couplings with the background fields $\vec{a}$
and $\vec{b}$. Here, the $\Sigma^{i}$’s are the well-known spin matrices which
appear in a covariant Dirac’s equation as the spatial components of
$\Sigma^{\mu\nu}=\frac{i}{4}\left[\gamma^{\mu},\gamma^{\nu}\right]$.
Since we are concerned with antineutrinos from beta-decay inside the Earth, we
have to neglect all the terms in the Lorentz-like force related to the charge,
and further on we shall focus on the terms that describe the interaction
between the background magnetic field, that we are interpreting as a
geomagnetic field, and the geoneutrinos, i.e., all the other terms are
considered zero,

$\mathfrak{F}^{i}=cf\gamma_{5}\left(\partial_{i}Z_{0}-\frac{1}{c}\partial_{t}Z^{i}\right)-cf\Sigma^{j}\left(\partial_{j}Z^{i}-\partial_{j}Z^{i}\right).$

The second term
$cf\Sigma^{i}\left(\partial_{i}Z^{j}-\partial_{j}Z^{i}\right)$, describes the
interaction between the spin of the geoneutrino and the background magnetic
field, $\vec{b}$.
Our next step is to work out its value, for this purpose, we shall calculate
the mean value of the term

$cf\left\langle\Sigma^{i}b^{i}\right\rangle.$ (5)

There are many problems related to calculating the mean value of the Eq.(5).
One of them is to know the value of the coupling constant $f$. Another
important problem is to find out the most convenient wave function to describe
geoneutrinos in this approach.

###  2.2 Experimental data and estimation of parameters

The wave function for a relativistic neutrino of mass $m$ and momentum
$\vec{p}$ through vacuum can be written as [13]

$\nu\left(t\right)_{vacuum}=\psi_{0}e^{-it\left(\frac{m^{2}}{2p}\right)};$ (6)

however, when the geoneutrios propagate through matter, there is a change in
the phase that affects even its oscillation. In the approximation $\Delta
m_{31}^{2}\sim\Delta m_{32}^{2}>>\Delta m_{21}^{2}$, the survival probability
$P_{ee}$ for a electronic antineutrino in the vacuum is given by [13]:

$P_{ee}=\sin^{4}\theta_{13}+\cos^{4}\theta_{13}\left[1-\sin^{2}2\theta_{12}\sin^{2}\left(\frac{1.267\Delta
m_{21}^{2}L}{4E}\right)\right],$ (7)

where $\theta$ is the angle between the mass eigenstates and the weak
eigenstates and $\Delta m^{2}$ is the difference of the squared mass
eigenvalues. This equation is valid in the 2-flavor approximation. Taking that
the diameter of the Earth is about $13000\ km$, for antineutrinos with energy
$\sim 3\ MeV$, the oscillation length,

$L_{0}\sim\pi c\hbar\frac{4E}{\Delta m^{2}_{21}},$ (8)

is of the order of $\sim 100\ km$, which is very small compared to Earth’s
diameter, and the effect of the neutrino oscillation to the total neutrino
flux is well-averaged, given an overall survival probability of $P_{ee}\sim
0.54$ [14] [15].
Matter is composed by quarks and electrons. The electron neutrino has a
especial behavior; it interacts with the electron via the exchange of the
charged boson $W^{+}$. Neutrino-neutrino interactions constitute neutral
currents and only take place with $Z^{0}$-exchange. The fact that electronic
neutrino interacts in a different form has important consequences; one of them
is that it feels a potential due to electrons and nucleons. Thus, the time-
dependent Hamiltonian for a neutrino propagating through matter gets an extra
term that modifies its phase and its related wave function as below:

$\nu\left(t\right)_{matter}=\psi_{0}e^{-it\left(\frac{m^{2}}{2p}+\sqrt{2}G_{F}N_{e}\right)},$
(9)

where $G_{F}$ is the Fermi constant, $N_{e}=\frac{\rho N_{0}Z}{A}$,
$\frac{Z}{A}$ is the average charge to mass ratio of the electrically neutral
matter and $N_{0}$ is the Avogadro number [13]. Since the electron density in
the Earth is not constant and it shows moreover drastic changes in
correspondence with boundaries of different Earth’s layers, the behavior of
the survival probability is not trivial and the equations have to be solved by
numerical tracing. It has been set in [16] [17] that this so-called matter
effect contribution to the average survival probability corresponds to an
increase of about $2\%$ and the spectral distorsion is below $1\%$. Thus, the
effect of the flavor oscillation on the total flux of geoneutrinos during
propagation is $\sim 0.55$ with a very small spectral distortion, completely
negligible for the precision of the current geoneutrino experiments.
Notwithstanding the description above, we are interested in setting up a
scenario that includes the geomagnetic field. We are going to assume that the
magnetic field is much smaller than $m^{2}_{p}/c$, where $m_{p}$ is the proton
mass, with typical values $\simeq 10^{-13}\ MeV^{2}$, that is, the mean value
of the geomagnetic field at Rio de Janeiro, Brazil, for example. As it has
been mentioned above, we wish to point out the rôle of the geomagnetic field
and its interaction through non-minimal coupling with antineutrinos from the
inverse beta-decay. For this purpose, we are going to use the formalism
developed in [18].
There is an important reason to adopt this approach: when we analyze the
energy spectrum, we can immediately see that the typical values of the
geomagnetic field are very tiny, so that it is reasonable to consider that the
effect of the geomagnetic field on geoneutrinos is perturbative,

$E^{2}=m^{2}+p_{z}^{2}+2ne\mathfrak{B}^{z},$ (10)

with $m$ the electron mass, and $p_{z}$ the electron momentum in the
$z-$direction and $\mathfrak{B}^{z}\sim 10^{-13}MeV^{2}$ stands for magnetic
field, as presented in [18]. In this approach, the magnetic field is
considered along the $z$\- direction. The maximum value of the parameter $n$
is

$n_{max}=int\left\\{\frac{1}{2e\mathfrak{B}}\left(\left[\left(m_{n}-m_{p}\right)+E_{\bar{\nu}}\right]^{2}-m^{2}\right)\right\\}\simeq
3\times 10^{16},$ (11)

where $m_{p}$ is the proton mass, $m_{n}$ is the neutron mass, $E_{\bar{\nu}}$
is the energy of the antineutrino. Since $n$ can take many values, we are
going to consider, in this work, only the maximum value, since in some cases
the minimum value is zero. Thus, the electron moment in the $z$-direction is
[18],

$p_{z}=\sqrt{\left[\left(M_{n}-M_{p}\right)+E_{\bar{\nu}}\right]^{2}-m^{2}-2ne\mathfrak{B}}.$
(12)

Replacing the values in the equation Eq.(12) above, we obtain $p_{z}\simeq
11.277\ MeV$. Now, we can calculate, from Eq.(10), the energy value
$E_{max}\simeq 15.9637MeV$. We calculate the currents $J_{i}$ and we get the
analytic form given below:

$\displaystyle J_{x}$ $\displaystyle=$ $\displaystyle 0$ $\displaystyle J_{y}$
$\displaystyle=$ $\displaystyle 0$ $\displaystyle J_{z}$ $\displaystyle=$
$\displaystyle\frac{2\left(e\mathfrak{B}^{z}\right)^{\frac{3}{2}}}{E_{n}+m}p_{z},$

with $J_{z\ max}\simeq 2\times 10^{-22}MeV$. Finally, we have calculated a set
of parameters that will help us to calculate the mean value for geomagnetic
field. In this formalism, the analytic expression for the mean value of the
geomagnetic field projection is,

$\displaystyle f<\Sigma^{x}\mathfrak{B}^{x}>_{x}$ $\displaystyle=$
$\displaystyle 0$ $\displaystyle f<\Sigma^{y}\mathfrak{B}^{y}>_{y}$
$\displaystyle=$ $\displaystyle 0$ $\displaystyle
f<\Sigma^{z}\mathfrak{B}^{z}>_{z}$ $\displaystyle=$ $\displaystyle
f\mathfrak{B}^{z}\left[1+\frac{1}{\left(E_{n}+m\right)^{2}}p_{z}^{2}-\frac{2ne\mathfrak{B}^{z}}{\left(E_{n}+m\right)^{2}}\right].$
(13)

Replacing the values in the expression above, we obtain for the mean value the
maximum and minimum values, namely:

$\displaystyle f<\Sigma^{z}\mathfrak{B}^{z}>_{z\ max}$ $\displaystyle=$
$\displaystyle 4.55921\times 10^{-14}\ MeV^{2},$ $\displaystyle
f<\Sigma^{z}\mathfrak{B}^{z}>_{z\ min}$ $\displaystyle=$ $\displaystyle
1.26148\times 10^{-14}\ MeV^{2};$ (14)

hence, in the spirit of the perturbative approach, we estimate the value of
the parameter to be $f\simeq 10^{-1}$.
There are three classes of bulk silicate earth (BSE) models: the
cosmochemical, geochemical and geodynamical [19]. These models provide
estimations for the global geoneutrino signal from the crust and mantle for
${}^{238}U$ and ${}^{232}Th$. The maximal value of the signal in the above
description is $35\ TNU$. We follow Fiorentini [20] and estimate the signal in
this approach arrive to $S_{T}\simeq 0.66\times 10^{48}\ TNU$ that is a unreal
value for the signal.

##  3 Final Considerations and Future Prospects

A formulation based on non-minimal couplings allows us to explore the
possibility of an interaction between geoneutrinos and the geomagnetic field.
Attempts to understand how neutrinos interacts in Nature are pursued in [21]
[22] [23] [24]. We investigate one face of the problem: their possible
interaction with the geomagnetic field. In this context, we take an approach
that highlights this aspect, and we place the interaction with the matter in a
second plane, once it is related to oscillation length. But, immediately,
there arise a question: is it possible that the combined action of the matter
and the geomagnetic field has a more important rôle in the description of this
problem? It is known now a days, that only oscillations between two families
is considered in geoneutrinos; this is the same scenario when we consider the
effects of the matter and the geomagnetic field together, even with a low
value of the geomagnetic field. This is, the consideration about the two
families still kept allow us to go deeper in this issue?
In the scenario that we have proposed, we can see that the interaction is
supported by a coupling parameter estimated to be of order $f\simeq 10^{-1}$,
and we conclude that the mean value of the magnetic interaction is $\sim
10^{-22}\ MeV^{2}$, a very small value. It is clear that, for to study this
kind of interaction, in this approach, we have to consider the possibility of
obtaining an enhanced flux of geoneutrinos in zones where the magnetic field
is more intense, ie, the geomagnetic poles. Possibly, a phenomenon like this
would have been detected already by the neutrinos detectors; but this is not
happening.
Measurement of geoneutrino obtained from $1353$ days at Laboratory Nazionali
del Gran Sasso (LNGS), in Italia, reports a signal of $38.8\pm 12.0\ TNU$ with
just a $6\times 10^{-6}$ probability for a null geoneutrino measurement [25].
Different analyses using cosmochemical, geochemical and geodynamical approach
seem to agree with the above value in the context of the BSE model. This model
is based on the supposition of the a chondritic Earth, we could ask us is this
measurement support only the consideration of a chondritic Earth with a ratio
$m_{Th}/m_{U}\sim 4$ [26] [27] [28]. We obtain a value for the signal of the
geoneutrinos very high, an unrealistic value. We interpret that the
overestimation of the geoneutrino signal value is due to the fact that the
model only uses the geomagnetic field value for the calculus of the current
and the flux of the geoneutrinos and not use any other geophysical parameter.
It is an issue to be pursued is to consider the specific case of massive
Majorana neutrinos, the family oscillations and the magnetism (actually, the
magnetic dipole moment) of massive neutrinos. It would be interesting, as a
follow-up of this paper, to discuss those matters in a geomagnetism scenario.
These are open problems that might be the object of future investigations.

## References

  * [1] N. Tolich, 2012. Geo-neutrino review. Nucl. Phys. Proc. Suppl. 229-232, 407.
  * [2] Buffett, Bruce A., 2002. Estimates of Heat Flow in the Deep Mantle Based on the Power Requirements for the Geodynamo. Geophysical Research Letters 29, no. 12: 7–1. doi:10.1029/2001GL014649.
  * [3] Davies, J. H., and D. R. Davies., 2010. Earth’s Surface Heat Flux. Solid Earth 1, no. 1: 5–24. doi:10.5194/se-1-5-2010.
  * [4] C. Jaupart, S. Labrosse, and J. C. Mareschal, 2007. “Temperatures, Heat and Energy in the Mantle of the Earth,” in Treatise of Geophysics, D. J. Stevenson, Ed., pp. 1-53, Elsevier, Amsterdam, The Netherlands, 2007.
  * [5] Lovering, J. F., and J. W. Morgan. 1964. Uranium and Thorium Abundances in Stony Meteorites: 1. The Chondritic Meteorites. Journal of Geophysical Research 69, no. 10: 1979–88. doi:10.1029/JZ069i010p01979.
  * [6] Shinotsuka, Kazunori, Hiroshi Hidaka, and Mitsuru Ebihara. 1995. Detailed Abundances of Rare Earth Elements, Thorium and Uranium in Chondritic Meteorites: An ICP-MS Study. Meteoritics 30, no. 6: 694–99. doi:10.1111/j.1945-5100.1995.tb01166.x.
  * [7] Herndon, J M. 1996. Substructure of the Inner Core of the Earth. Proceedings of the National Academy of Sciences of the United States of America 93, no. 2: 646–48.
  * [8] Murthy, V. Rama, Wim van Westrenen, and Yingwei Fei. 2003. Experimental Evidence That Potassium Is a Substantial Radioactive Heat Source in Planetary Cores. Nature 423, no. 6936: 163–65. doi:10.1038/nature01560.
  * [9] Abe, S., T. Ebihara, S. Enomoto, K. Furuno, Y. Gando, K. Ichimura, H. Ikeda, et al. 2008. Precision Measurement of Neutrino Oscillation Parameters with KamLAND. Physical Review Letters 100, no. 22: 221803. doi:10.1103/PhysRevLett.100.221803.
  * [10] Fujikawa, Kazuo, and Robert E. Shrock. 1980. Magnetic Moment of a Massive Neutrino and Neutrino-Spin Rotation. Physical Review Letters 45, no. 12: 963–66. doi:10.1103/PhysRevLett.45.963.
  * [11] Low, F. E. 1965. Heavy Electrons and Muons. Physical Review Letters 14, no. 7: 238–39. doi:10.1103/PhysRevLett.14.238.
  * [12] Barshay, S. 1965. CP and Non-Minimal Electromagnetic Interactions. Physics Letters 17, no. 1: 78–80. doi:10.1016/0031-9163(65)90656-6.
  * [13] Felix Bohem and Petr Vogel. 1987. Physics of Massive Neutrinos., Cambridge University Press., Great Britain, Cambridge.
  * [14] Ludhova, L., and S. Zavatarelli. 2013. Studying the Earth with Geoneutrinos. Advances in High Energy Physics: e425693. doi:10.1155/2013/425693.
  * [15] Fogli, G. L., E. Lisi, A. Marrone, D. Montanino, A. Palazzo, and A. M. Rotunno. 2012. Global Analysis of Neutrino Masses, Mixings, and Phases: Entering the Era of Leptonic $CP$ Violation Searches. Physical Review D 86, no. 1: 013012. doi:10.1103/PhysRevD.86.013012.
  * [16] S. Enomoto. 2005. “Neutrino geophysics and observation of Geoneutrinos at KamLAND” [PhD. thesis], Tohoku University, Honshu, Japan.
  * [17] Honda, Minako, Yee Kao, Naotoshi Okamura, and Tatsu Takeuchi. 2006. A Simple Parameterization of Matter Effects on Neutrino Oscillations. arXiv:hep-ph/0602115. http://arxiv.org/abs/hep-ph/0602115.
  * [18] Bhattacharya, Kaushik, and Palash B. Pal. 1999. Inverse Beta-Decay in Magnetic Fields. arXiv:hep-ph/9911498. http://arxiv.org/abs/hep-ph/9911498.
  * [19] Šrámek, Ondřej, William F. McDonough, Edwin S. Kite, Vedran Lekić, Stephen T. Dye, and Shijie Zhong. 2013.Geophysical and Geochemical Constraints on Geoneutrino Fluxes from Earth’s Mantle. Earth and Planetary Science Letters 361: 356–66. doi:10.1016/j.epsl.2012.11.001.
  * [20] Fiorentini, Gianni, Marcello Lissia, and Fabio Mantovani. 2007. Geo-Neutrinos and Earth’s Interior. Physics Reports 453, no. 5–6 : 117–72. doi:10.1016/j.physrep.2007.09.001.
  * [21] Mukhopadhyay, Banibrata. Neutrino Oscillations under Gravity: Mass Independent Oscillation. 2003. arXiv:hep-ph/0307167, July 12. http://arxiv.org/abs/hep-ph/0307167.
  * [22] Miller, Jonathan, and Roman Pasechnik. 2013. Quantum Gravity Effect on Neutrino Oscillations in a Strong Gravitational Field. arXiv:1305.4430 [gr-Qc, Physics:hep-Ph, Physics:hep-Th], May 19. http://arxiv.org/abs/1305.4430.
  * [23] Plyatsko, Roman, and Mykola Fenyk. Highly Relativistic Spin-Gravity Coupling for Fermions. 2015. arXiv:1503.08415 [astro-Ph, Physics:gr-Qc, Physics:hep-Ex, Physics:hep-Th], March 29. http://arxiv.org/abs/1503.08415.
  * [24] Dvornikov, Maxim. 2010. Field Theory Description of Neutrino Oscillations. arXiv:1011.4300 [hep-Ph, Physics:hep-Th], November 18. http://arxiv.org/abs/1011.4300.
  * [25] Belini et al. 2013. Measurement of Geo-Neutrinos from 1353 Days of Borexino. Physics Letters B 722, no. 4–5: 295–300. doi:10.1016/j.physletb.2013.04.030.
  * [26] Bezrukov, L., and V. Sinev. 2014. Geoneutrinos and Hydridic Earth (or Primordially Hydrogen-Rich Planet). 2014. arXiv:1405.3161 [astro-Ph], May 13. http://arxiv.org/abs/1405.3161.
  * [27] Rusov, V. D., V. N. Pavlovich, V. N. Vaschenko, V. A. Tarasov, T. N. Zelentsova, V. N. Bolshakov, D. A. Litvinov, S. I. Kosenko, and O. A. Byegunova. 2007. Geoantineutrino Spectrum and Slow Nuclear Burning on the Boundary of the Liquid and Solid Phases of the Earth’s Core.” Journal of Geophysical Research: Solid Earth 112, no. B9: B09203. doi:10.1029/2005JB004212.
  * [28] Meijer, R. J. de, and W. van Westrenen. 2008. The Feasibility and Implications of Nuclear Georeactors in Earth’s Core-Mantle Boundary Region.” South African Journal of Science 104, no. 3–4: 111–18.