# A New Order Parameter for the Higgs Transition in
$SU(2)$-Higgs Theory

Mridupawan Deka<EMAIL_ADDRESS>Bogoliubov Laboratory of Theoretical
Physics, JINR, 141980 Dubna, Russia  Institute of Mathematical Sciences,
Chennai 600113, India  Sanatan Digal<EMAIL_ADDRESS>Institute of
Mathematical Sciences, Chennai 600113, India

###### Abstract

We investigate the Higgs transition within the four dimensional $SU(2)-$
gauge-Higgs model in search for an order parameter as a function of the Higgs
field hopping parameter, $\kappa$, using Lattice technique. We measure the
Higgs condensate after applying Landau Gauge Fixing and study the
corresponding susceptibility, magnetization and fourth order Binder cumulant
using four different spatial volumes with $N_{\tau}\\!=\\!2$. The computation
is carried out with gauge coupling, $\beta_{g}\\!=\\!8$, for a range of scalar
self-coupling, $\lambda\\!=\\!\\{0.00010,0.00350\\}$, with emphasis near the
critical end-point. Finite size scaling analysis of the gauge fixed condensate
and its cumulants agree with the standard $3$d Ising values $\nu=0.62997$,
$\beta/\nu=0.518$, $\gamma/\nu=1.964$ at $\lambda=0.00150$. These results are
in agreement with previous studies suggesting $3$d Ising universality class.
The numerical results also indicate that, at the transition point, the gauge
fixed condensate vanishes in the infinite volume limit.

###### pacs:

11.10.Wx,11.15.Ha,11.15.-q

##  I Introduction

The understanding of the phase structure of $SU(2)-$Higgs theory has crucial
cosmological and experimental consequences Kuzmin:1985mm ; Matveev:1988pj ;
Rubakov:1996vz ; Fodor:1998if . Over the years, several non-perturbative
lattice studies have been devoted in this regard Damgaard:1986qe ;
Evertz:1986af ; Bunk:1992xt ; Farakos:1994xh ; Kajantie:1995kf ;
Rummukainen:1998as ; Bunk:1992kf ; Bunk:1994xs ; Kajantie:1996mn ;
Gurtler:1997hr ; Gurtler:1997ki ; Buchmuller:1994qy ; Karsch:1996aw ;
Karsch:1996yh ; Fodor:1994dm ; Fodor:1994sj ; Csikor:1995jj ; Aoki:1996cu ;
Aoki:1999fi ; Bonati:2009pf ; DOnofrio:2014rug ; DOnofrio:2015gop ;
DOnofrio:2012yxq ; Laine:2012jy ; Gould:2022ran  (For reviews, please see
Kajantie:1994mt ; Jansen:1995yg ; Rummukainen:1996sx ). These studies find
that the nature of Higgs transition varies with Higgs quartic coupling
($\lambda$). With increase in $\lambda$, the strength of this transition
weakens and eventually becomes a crossover for large $\lambda$. The $1$st
order transition and crossover regions are separated by a $2$nd order critical
end-point at $\lambda_{c}$ Aoki:1999fi . There are similar examples in QCD,
spin-models, condensed matter systems etc. where lines/surfaces of second
order transition that separate $1$st order and crossover regions in the phase
diagram Pisarski:1983ms ; Karsch:2000xv ; Alonso:1993tv . Many attempts have
been made to study the universality class of such critical end-points in
different theories Pisarski:1983ms ; Alonso:1993tv ; Rummukainen:1998as ;
Karsch:2000xv . These end-points are expected to belong to the universality
class of the Ising model Alonso:1993tv ; Rummukainen:1998as ; Janke:1996qb ;
Karsch:2000xv ; Wilding:1994zkq ; Wilding:1994ud ; Wilding:1996 ; Rehr:1973zz
. In the electroweak theory, the critical end-point in dimensionally reduced
$SU(2)-$Higgs theory is studied in three dimensions Rummukainen:1998as . The
analysis of the critical behavior is done by identifying energy and
magnetization like observables Alonso:1993tv . The distribution of these
observables, higher cumulants and the extracted critical exponents, clearly
show that the end-point is in the $3$d Ising universality class.

In ref. Rummukainen:1998as , a linear combinations of various gauge invariant
observables plays the role of magnetization. Another orthogonal linear
combination plays the role of energy like observable. These linear
combinations change with the bare parameters of the theory. Though gauge
invariant observables can be used to study nature of transitions, it is
desirable to have an order parameter for the Higgs transition. In the absence
of gauge symmetry, the Higgs condensate (volume average of the Higgs field)
plays the role of an order parameter. On other hand, gauge symmetry renders
the condensate unphysical. However, it can be made well behaved by a suitable
choice of a gauge. It is expected that the choice of a gauge will not affect
physical quantities such as critical exponents. Note that the necessity of a
suitable order parameter exists in other gauge theories, such as in high
density QCD to describe the normal to color superconducting transition.

In the present work, we propose that the Higgs condensate in the Landau gauge,
denoted by $\Phi^{g}$, as an order parameter for the Higgs transition in
$SU(2)-$Higgs theory. In $3+1$ dimensions, $\Phi^{g}$ and its cumulants are
computed using Monte Carlo simulation of the partition function by varying the
Higgs hoping parameter ($\kappa$) and the quartic coupling. It is observed
that the partition function average of $\Phi^{g}$ behaves similar to that of
magnetization in spin models. For small $\lambda$, the Higgs transition is of
$1$st order, and $\Phi^{g}$ is found to vary discontinuously across the
transition point. With the increase in $\lambda$, the Higgs transition weakens
which is also seen in $\Phi^{g}$ behavior. The jump in $\Phi^{g}$ across the
transition point decreases with $\lambda$. At a critical point $\lambda_{c}$,
$\Phi^{g}$ varies continuously while its various cumulants show singular
behavior. The observed scaling behavior is consistent with $3$d Ising
universality class and deviates from the universality class of the $3$d $O(4)$
spin models Kanaya:1994qe ; Toussaint:1996qr . The Finite Size Scaling (FSS)
of $\Phi^{g}$ suggests that it vanishes at the critical point in the infinite
volume limit. It also vanishes in the high temperature or Higgs symmetric
phase at least for $\lambda<\lambda_{c}$. We would like to mention here that
$\Phi^{g}$, as an order parameter, reasonably reproduces the previous results.
A linear combination along the lines of previous studies may further fine tune
the results.

The draft is organized as follows. In Sec. II, we provide the lattice
definition of the $SU(2)$-Higgs theory. In Sec. III, we discuss the numerical
procedures for our work. This includes the Landau Gauge Fixing in Lattice. In
Sec. V, we analyze our data using FSS, while we present the magnetization
study at various $\lambda$’s in Sec. VI. Finally, we present the discussion
and conclusion of our study in Sec. VII.

##  II Lattice Formalism

The continuum action of $SU(2)$ gauge-Higgs theory is discretized in Lattice
as,

$\displaystyle\displaystyle S$ $\displaystyle=$
$\displaystyle\sum_{x}\left[\sum_{\mu>\nu}\frac{\beta}{2}\,\mbox{Tr}\,U_{x,\mu\nu}+\sum_{\mu}\kappa\,\mbox{Tr}\,\left(\phi^{\dagger}_{x}\,U_{x,\mu}\,\phi_{x+{\hat{\mu}}}\right)-\frac{1}{2}\,\mbox{Tr}\left(\phi^{\dagger}_{x}\phi_{x}\right)-\lambda\left(\frac{1}{2}\,\mbox{Tr}\left(\phi^{\dagger}_{x}\phi_{x}\right)-1\right)^{2}\right]\,,$

where $U_{x,\mu}\in SU(2)$ is a link variable, and $U_{x,\mu\nu}$ is a product
of four link variables which form a plaquette, and $\phi_{x}$ is $2\otimes 2$
matrix in isospin space describing the Higgs scalar field. The bare parameters
$\beta_{g}\equiv 4/g^{2}$ is the gauge coupling, $\lambda$ is the scalar
quartic coupling, and $\kappa$ is the scalar hopping parameter which is
related to the bare mass square by the relation,
$\mu^{2}_{0}=(1-2\lambda)\,\kappa^{-1}-8$.

In order to search for an order parameter, we apply Landau gauge fixing,
$\partial_{\mu}\,A_{\mu}\\!=\\!0$, and construct the condensate, $\Phi^{g}$.
We then study the behavior of the magnetization, susceptibility and fourth
order Binder cumulant of $\Phi^{g}$ (defined below) at various values of
$\kappa$ and $\lambda$ while $\beta_{g}\\!=\\!8$ is kept fixed.

Magnetization $\displaystyle:$ $\displaystyle\hskip
14.22636pt\Phi^{g}\,\equiv\,\langle\tilde{\Phi}^{g}\rangle\,=\,\left\langle\frac{1}{N_{s}^{3}N_{t}}\sqrt{\sum_{i=1}^{4}\left(\displaystyle\sum_{x}\phi^{g}_{i,x}\right)^{2}}\right\rangle,$
(2) $\displaystyle\mbox{Susceptibility,}\,\,\chi^{c}\hskip 14.22636pt$
$\displaystyle:$ $\displaystyle\hskip
14.22636ptV\,\left(\langle{\tilde{\Phi}^{g}}{}^{2}\rangle-\langle\tilde{\Phi}^{g}\rangle^{2}\right),$
(3) $\displaystyle\mbox{Total Susceptibility,}\,\,\chi\hskip 14.22636pt$
$\displaystyle:$ $\displaystyle\hskip
14.22636ptV\,\langle{\tilde{\Phi}^{g}}{}^{2}\rangle,$ (4)
$\displaystyle\mbox{Binder Cumulant,}\,\,B_{4}\hskip 14.22636pt$
$\displaystyle:$ $\displaystyle\hskip
14.22636pt1-\frac{\langle{\tilde{\Phi}^{g}}{}^{4}\rangle}{3\,\langle{\tilde{\Phi}^{g}}{}^{2}\rangle^{2}}.$
(5)

where $\langle\cdots\rangle$ stands for gauge average, and $V=N_{s}^{3}$ is
the spatial volume of the lattice. We use the superscript, $c$, in Eq. (3) to
specify it as the Connected Susceptibility.

We devote a significant part of our study near the end-point for $1$st order
transition, namely at $\lambda=0.00116,0.00132$ and $0.00150$. Note that the
first value of $\lambda$ is found to be the critical end-point in Ref.
Aoki:1999fi . The two higher $\lambda$ values lie around $\sim+2\sigma$ of the
first one. The reason to include these higher values of $\lambda$ is discussed
in the Sec. V.1. We list the run parameters in Table 1.

##  III Numerical Procedure and Parameters

For Monte Carlo simulations, we use the pure gauge part of the publicly
available MILC code milc  and modify it to accommodate the Higgs fields. To
update the gauge fields, we first use the standard heat bath algorithm
Creutz:1980zw ; Cabibbo:1982zn , and then update Higgs fields using pseudo
heat bath algorithm Bunk:1994xs . We then again update the gauge fields using
$4$ overrelaxation steps Whitmer:1984he  after which Higgs fields are updated
again using pseudo heat bath algorithm. To reduce auto-correlation between
successive configurations, we carry out $50$ cycles of this updating procedure
between two subsequent measurements. We test our codes to reproduce histogram
profiles similar to Ref. Aoki:1996cu . For our simulations, we use four
spatial volumes, $N_{s}\\!=\\!20,24,32,40$ while the temporal lattice extent
is kept at $N_{\tau}=2$ Aoki:1999fi .

To implement Landau gauge fixing, we maximize the following quantity

$\displaystyle H$ $\displaystyle=$
$\displaystyle\displaystyle\sum_{x}\sum^{4}_{\mu=1}\mbox{Re}\,\left[\mbox{Tr}\,U_{x,\mu}\right]\,.$
(6)

In order to maximize $H$, we separately apply two different methods as
described in Aoki:2003ip . One of the methods is the standard $SU(2)$ subgroup
method, and the other one is Overrelaxed Steepest Descent method. This was
done to check whether both the methods arrive at the same global maxima. In
addition, we use single, double and quadruple precisions in our simulations so
that we are able to set our convergence conditions accordingly. We find that
both the methods for maximization of $H$ as well as the different precisions
produce similar results. Eventually, as described in Aoki:2003ip , we combine
both $SU(2)$ subgroup and Overrelaxed Steepest Descent methods to speed up
gauge fixing. First, we apply $SU(2)$ subgroup method to bring the
configuration close to the maximum, and then apply Overrelaxed Steepest
Descent method for the final convergence. We set our overrelaxation parameter,
$\omega=1.98$, and we use double precision for all our simulations. For
convergence of $SU(2)$ subgroup method, we measure

$\displaystyle h$ $\displaystyle\sim$
$\displaystyle\frac{H}{N^{3}_{s}\,N_{\tau}},$ (7)

We set the convergence condition at the $i$th iteration as,

$\displaystyle\left|h_{i}-h_{i-1}\right|$ $\displaystyle<$ $\displaystyle
10^{-12}\,.$ (8)

For the Steepest Descent method, we measure

$\displaystyle\Delta$ $\displaystyle=$
$\displaystyle\displaystyle\frac{1}{N^{3}_{s}\,N_{\tau}}\sum_{x}\frac{1}{2}\sum^{4}_{\mu=1}\frac{1}{4}\,\mbox{Tr}\,\left[\Delta_{x,\mu}\Delta^{\dagger}_{x,\mu}\right],$
(9)

where

$\displaystyle\Delta_{x,\mu}$ $\displaystyle=$ $\displaystyle
U_{x-\hat{\mu},\mu}-U_{x,\mu}-\mbox{h.c.}-\mbox{trace}\,.$ (10)

We set the convergence condition for $\Delta_{j}$ at the $j$th iteration as,

$\displaystyle\left|\Delta_{j}\right|$ $\displaystyle<$ $\displaystyle
10^{-15}.$ (11)

Simulations are performed for a wide range of $\kappa$ and $\lambda$ values.
We focus our studies at $\lambda\\!=0.00116,0.00132$ and $0.00150$. For each
of them, we generate $2.5$ million trajectories out of which gauge fixing is
performed on every $50$-th trajectory at each $\kappa$. This gives us a total
$50\,000$ gauge-fixed configurations for each $\kappa$. The parameters are
listed in Table 1. For each of other values of $\lambda$, we generate
$500\,000$ trajectories instead out of which every $50$-th trajectory is
selected for gauge fixing giving $10\,000$ gauge-fixed configurations for each
$\kappa$. The reduced statistics is due to the limitations in computational
resources. The parameters are listed in Table 2. All the error analyses in
this study are carried out using Jackknife method with bin size determined
such a way that the total number of bins is $100$.

Table 1: Simulation parameters near $\lambda_{c}$. $\beta_{g}$ | $\lambda$ | $N_{s}$ | $N_{\tau}$ | $\kappa$ | Confs. | Trajectories
---|---|---|---|---|---|---
between
consecutive
confs.
8.0 | 0.00116 | 20, 24, 32, 40 | 2 | {0.129390, 0.129550} | 50 000 | 50
0.00132 | 20, 24, 32, 40 | 2 | {0.129420, 0.129820} | 50 000 | 50
0.00150 | 20, 24, 32, 40 | 2 | {0.129500, 0.130000} | 50 000 | 50

##  IV $M$-like and $E$-like observables

The $SU(2)-$Higgs theory is similar to the case of the liquid-gas system, i.e
none of the gauge invariant observables can be directly identified with the
$M$-like and $E$-like observables. In the previous studies of the critical
endpoint, observables such as
$S_{K}=\displaystyle\sum_{x,\mu}\mbox{Tr}\left(\phi^{\dagger}_{x}\,U_{x,\mu}\,\phi_{x+{\hat{\mu}}}\right)$,
$S_{\phi}=\displaystyle\sum_{x}\left(\frac{1}{2}\,\mbox{Tr}\
\left(\phi^{\dagger}_{x}\phi_{x}\right)-1\right)^{2}$ etc. have been used to
define the $M$-like and $E$-like observables. It has been pointed out that
$S_{K}$, $S_{\phi}$ etc. are not necessarily orthogonal in the vicinity of the
critical point, i.e. $\langle\Delta S_{K}\Delta S_{\phi}\rangle\neq 0$ where
$\Delta S_{K}=S_{K}-\langle S_{K}\rangle$ and $\Delta
S_{\phi}=S_{\phi}-\langle S_{\phi}\rangle$ Rummukainen:1998as ; Alonso:1993tv
; Karsch:2000xv ; Wilding:1994zkq ; Wilding:1994ud ; Wilding:1996 . However,
linear and orthogonal combinations of them have been found to behave as
$M$-like and $E$-like observables. In the present work, we include the Landau
gauge fixed Higgs condensate, $\Phi^{g}$, along with other gauge invariant
observables to define $M$-like and $E$-like observables. The studies of
correlations between $\Phi^{g}$ and other gauge invariant observables seem to
suggest that $\langle\Delta\Phi^{g}\Delta S_{K}\rangle\approx 0$ and
$\langle\Delta\Phi^{g}\Delta S_{\phi}\rangle\approx 0$, are satisfied within
the uncertainties.

We compute the correlation relations, $\langle\Delta\Phi^{g}\Delta
S_{K}\rangle$ and $\langle\Delta\Phi^{g}\Delta S_{\phi}\rangle$, over the
whole range of respective $\kappa$ values at $\lambda=0.00116$ and $0.00150$.
They are plotted in Figs. 1 and 2, respectively. We clearly see that there is
no correlation between the $\Phi^{g}$ and other observables. We also observe
that the sizes of uncertainties decrease as the volume increases. This
suggests that $\Phi^{g}$ is a $M$-like observable.

We also study the behaviour of the $\Phi^{g}$ distributions vs. $S_{K}$ as
well as $S_{\phi}$ distributions. In Fig. 3, we plot $\Phi^{g}$ vs. $S_{K}$
and $S_{\phi}$ for $\lambda=0.00116$ and $\kappa=0.129465$. In Fig. 4, the
same quantities are plotted for $\lambda=0.001500$ and $\kappa=0.129765$. From
these figures, we see that these distributions behave similarly to those of
energy vs. magnetization distributions of the $3$d Ising Model. Note that
$\Phi^{g}$ has four components. The individual components have no physical
significance, as there are no Goldstone modes in this case. The magnitude of
$\Phi^{g}$ is physical and as a consequence, we can only reproduce the part of
the Ising model figure for which $M>0$. In Figs. 3 and 4, we observe that the
distributions are not linear, unlike the distribution of $S_{k}$ vs.
$S_{\phi}$. In this situation, no unique linear combination is possible.

Figure 1: Correlations between (a) $\Phi^{g}$ and $S_{K}$ (b) $\Phi^{g}$ and
$S_{\phi}$ at $\lambda=0.00116$.

Figure 2: Correlations between (a) $\Phi^{g}$ and $S_{K}$ (b) $\Phi^{g}$ and
$S_{\phi}$ at $\lambda=0.00150$.

Figure 3:  $\Phi^{g}$ vs. (a) $S_{K}$ and (b) $S_{\phi}$ at $\lambda=0.001160$
and $\kappa=0.129465$.

Figure 4:  $\Phi^{g}$ vs. (a) $S_{K}$ and (b) $S_{\phi}$ at $\lambda=0.001500$
and $\kappa=0.129765$.

##  V Finite Size Scaling Analysis

Figure 5:  (a) Magnetization (b) Connected susceptibility (c) Total
susceptibility, and (d) Binder cumulant as functions of $\kappa$ at
$\lambda=0.00116$. The dotted lines are not fitted and are for eye guides
only.

Figure 6: The values of ${\kappa_{\chi}}_{\mbox{\scriptsize max}}$ as a
function of $N_{s}^{-1/\nu}$ at $\lambda=0.00116$ for $N_{s}=20,24,32$ and
$40$. We set $\nu=0.629971$. The intersection of the fitted straight line
provides the value of $\kappa_{c}$.

Figure 7: Scaled (a) Magnetization (b) Connected susceptibility (c) Total
susceptibility, and (d) Binder cumulant as functions of $\kappa$ at
$\lambda=0.00116$.

###  V.1 $\lambda=0.00116$

In Fig. 5, we plot magnetization, susceptibility, total susceptibility and
Binder cumulant for $\lambda=0.00116$. We see that $\kappa$ acts like an
“inverse temperature” Wurtz:2009gf  for all the four quantities, and there are
indications that these four quantities behave similar to that of Ising Model.
However, it is necessary to study further by scaling all the four quantities
with the standard Ising exponents and see whether the scaled quantities follow
the FSS behavior.

In order to carry out the scaling procedure, the corresponding critical value
of $\kappa$ needs to be evaluated. To find $\kappa_{c}$, we first estimate the
value of $\chi^{c}_{\mbox{\scriptsize max}}$ for each volume and obtain the
corresponding ${\kappa_{\chi}}_{\mbox{\scriptsize max}}$. The standard
procedure to find $\chi^{c}_{\mbox{\scriptsize max}}$ is the Rewieghting
method. It can be seen from Fig. 5 that we have a reasonable amount of data
near the peak point for every volume. Thus, we instead use the Cubic Spline
Interpolation method to generate a few hundred points close to
${\kappa_{\chi}}_{\mbox{\scriptsize max}}$ for every Jackknife sample. From
these sets of interpolated points, we find the values of
$\chi^{c}_{\mbox{\scriptsize max}}$ and the corresponding
${\kappa_{\chi}}_{\mbox{\scriptsize max}}$ for each volume. The value of
$\kappa_{c}$ is then obtained by using the following FSS relation given by,

$\displaystyle{\kappa_{\chi}}_{\mbox{\scriptsize max}}^{-1}$ $\displaystyle=$
$\displaystyle\kappa_{c}^{-1}+a\,N_{s}^{-1/\nu}.$ (12)

By using the standard value of $\nu=0.629971$ for $\\!3$d Ising Model, we get
the value of $\kappa_{c}$ from the linear fit of
${\kappa_{\chi}}_{\mbox{\scriptsize max}}$ as a function of $N_{s}^{-1/\nu}$
as (See Fig. 6),

$\displaystyle\kappa_{c}$ $\displaystyle=$ $\displaystyle 0.1294655(5).$ (13)

Using the value of $\kappa_{c}$ from Eq. (13), and the standard $\\!3$d Ising
values of $\gamma=1.237075$ and $\beta=0.326419$, the quantities are scaled.
They are plotted in Fig. 7.

It can be clearly seen from Fig. 7 that quantities do not properly obey the
scaling behavior. The differences are more visible in the cases of
susceptibility and Binder cumulant. Further, no improvement in the scaling is
seen by changing the universality class. We thus extend our study to other
$\lambda$ values for which scaling can be seen. We considered two other values
of $\lambda$, namely $\lambda=0.00132$ and $0.00150$, which are around $\sim
1\sigma$ and $\sim 2\sigma$ from $\lambda=0.00116$, respectively.

Figure 8:  (a) Magnetization (b) Connected susceptibility (c) Total
susceptibility, and (d) Binder cumulant as functions of $\kappa$ at
$\lambda=0.00132$. The dotted lines are for eye guides only.

Figure 9: Scaled (a) Magnetization (b) Connected susceptibility (c) Total
susceptibility, and (d) Binder cumulant as functions of $\kappa$ at
$\lambda=0.00132$.

Figure 10:  (a) Magnetization (b) Connected susceptibility (c) Total
susceptibility, and (d) Binder cumulant as functions of $\kappa$ at
$\lambda=0.00150$. The dotted lines are for eye guides only.

Figure 11: Scaled (a) Magnetization (b) Connected susceptibility (c) Total
susceptibility, and (d) Binder cumulant as functions of $\kappa$ at
$\lambda=0.00150$.

###  V.2 $\lambda=0.00132$ and $0.00150$

We compute magnetization, susceptibility, total susceptibility and Binder
cumulant for $\lambda=0.00132$ and $0.00150$. They are plotted in Figs. 8 and
10, respectively. We again see Ising like behavior of these quantities. Note
that the value of $\kappa_{c}$ changes with the value of $\lambda$. We follow
the same procedure as described in Sec. V.1 to determine $\kappa_{c}$ for
$\lambda=0.00132$ and $0.00150$. Thus, we get,

$\displaystyle\kappa_{c}$ $\displaystyle=$ $\displaystyle
0.1296071(8)\,\,\,\,\,\mbox{for $\lambda=0.00132$},$ (14)
$\displaystyle\kappa_{c}$ $\displaystyle=$ $\displaystyle
0.1297642(7)\,\,\,\,\,\mbox{for $\lambda=0.00150$}.$ (15)

Using these values of $k_{c}$ and the standard values of the exponents for
$\\!3$d Ising Model, we scale the quantities as in Sec. V.1. They are plotted
in Figs. 9 and 11, respectively. By comparing Figs. 7, 9 and 11, we see
scaling for $\lambda=0.00150$. This is an indication that $\Phi^{g}$ is a good
order parameter. Also this suggests that the critical end-point is close to
$\lambda=0.00150$.

##  VI Magnetization at various other $\lambda$’s

Table 2: Run parameters for other $\lambda$’s. $\beta_{g}$ | $\lambda$ | $N_{s}$ | $N_{\tau}$ | $\kappa$ | Confs. | Trajectories
---|---|---|---|---|---|---
between
consecutive
confs.
8.0 | 0.00010 | 20, 24, 32, 40 | 2 | {0.127800, 0.128890} | 10 000 | 50
0.00050 | 20, 24, 32, 40 | 2 | {0.127800, 0.129200} | 10 000 | 50
0.00080 | 20, 24, 32, 40 | 2 | {0.129000, 0.129400} | 10 000 | 50
0.00095 | 20, 24, 32, 40 | 2 | {0.129000, 0.129400} | 10 000 | 50
0.00100 | 20, 24, 32, 40 | 2 | {0.129120, 0.129500} | 10 000 | 50
0.00170 | 20, 24, 32, 40 | 2 | {0.129690, 0.130300} | 10 000 | 50
0.00220 | 20, 24, 32, 40 | 2 | {0.129450, 0.130600} | 10 000 | 50
0.00270 | 20, 24, 32, 40 | 2 | {0.130000, 0.131340} | 10 000 | 50
0.00310 | 20, 24, 32, 40 | 2 | {0.129500, 0.135000} | 10 000 | 50
0.00350 | 20, 24, 32, 40 | 2 | {0.127500, 0.137000} | 10 000 | 50

Figure 12:  Magnetization at various values of $\lambda$. (a)
$\lambda=0.00010$, (b) $\lambda=0.00050$, (c) $\lambda=0.00080$, (d)
$\lambda=0.00095$, (e) $\lambda=0.00100$, (f) $\lambda=0.00170$, (g)
$\lambda=0.00310$, and (h) $\lambda=0.00350$.

We further investigate whether $\Phi^{g}$ displays the behavior of a
conventional order parameter as a function of $\lambda$. We study the
magnetization for various $\lambda$ values between $\\{0.0001,0.0035\\}$. The
simulation parameters for this study are listed in Table 2. As mentioned in
Sec. III, the number of configurations is restricted to $10\,000$ only for
each combination of $\lambda$ and $\kappa$.

From Figs. 12 and 12 at $\lambda=0.00010$ and $0.00050$, respectively, we
clearly see that $\Phi^{g}$ displays a $1$st order transition as a function of
$\kappa$. Finally, we see from Figs. 12 and 12 that the volume dependence
nearly disappears at the values of $\lambda=0.00310$ and $0.00350$,
respectively, indicating a cross-over.

We also evaluate ${\kappa_{\chi}}_{\mbox{\scriptsize max}}$’s from the maximum
value of susceptibility for $N_{s}=40$ for $\lambda=\\{0.0008,0.0035\\}$. They
are plotted in Fig. 13.

Figure 13: The behavior of ${\kappa_{\chi}}_{\mbox{\scriptsize max}}$ for
$N_{s}=40$ as a function of $\lambda$.

##  VII Discussion and Conclusion

In this work, we study the behavior of the volume average of the Landau gauge-
fixed $\Phi^{g}$ and related quantities in $SU(2)-$Higgs theory. The aim of
our study is to see if $\Phi^{g}$ behaves like an order parameter for the
Higgs transition. We consider $\beta_{g}=8$ as in Aoki:1999fi , and study
$\Phi^{g}$ at various values of $\lambda$, from $\lambda=0.00010$ to
$\lambda=0.0035$. Our results show that, $\Phi^{g}$ for different $\lambda$
accurately identifies the transitions points ($\kappa_{c}$’s). Across first
order transition, $\Phi^{g}$ varies discontinuously, as expected from an order
parameter. This is clearly seen at $\lambda=0.00010$ and$\lambda=0.00050$ in
agreement with previous study Aoki:1999fi . From $\lambda=0.0008$ onward, we
see weakening of transition, with the jump in $\Phi^{g}$ decreasing with
$\lambda$. At much higher values of $\lambda$, such as $\lambda=0.00310$ and
$\lambda=0.00350$, we see a gradual disappearance of volume dependence
signifying a possible crossover. It is interesting to note that a different
lattice study Wurtz:2013ova  at zero temperature performed at
$\beta_{g}=8,\lambda=0.0033,\kappa=0.131$ with a lattice size of $24^{3}\times
48$ has found the physical Higgs mass to be $123(1)$ GeV which is closer to
the experimental value.

A significant part our numerical simulations are devoted to study behavior of
$\Phi^{g}$ near the end-point. The magnetization, susceptibility, total
susceptibility and binder cumulant etc. of $\Phi^{g}$ as functions of $\kappa$
give an impression of an order parameter like behavior at $\lambda_{c}$. The
scaling of these quantities for $\nu=0.62997$, $\beta/\nu=0.518$ and
$\gamma/\nu=1.964$, shows that the end-point belongs to the $3$d Ising
universality class. From the scaling of the Binder cumulant from Fig. 11, we
find its value at $(\lambda_{c},\kappa_{c})$ to be consistent with the
standard $3$d Ising value of $\sim 0.47$ Ferrenberg:1991 ; Hasenbusch:1998gh .

We mention here that the FSS of $\Phi^{g}$ and its various cumulants suggests
that the $\\!3$d Ising like behavior seen close to $\lambda=0.00150$ which is
within $2\sigma$ of the end-point found in Aoki:1999fi . The FSS of $\Phi^{g}$
seen in Fig. 11 suggests that exactly at the critical point
$(\lambda_{c},\kappa_{c})$ $\Phi^{g}$ vanishes. Also the numerical results
suggest that for $\lambda<0$ in the Higgs symmetric phase and at $\kappa_{c}$,
$\Phi^{g}$ vanishes in the infinite volume limit. This is in contrast to
behavior of magnetization in $3$-state Potts model, where the magnetization
vanishes in the symmetric phase only in the absence of external field. In the
presence of external field, in the symmetric phase at transition point,
magnetization increases with external field. This behavior of magnetization
$vs$ external field can be be obtained by adding a simple explicit breaking
term to Ginzburg-Landau type of mean field free energy. Such an effective
description of $\Phi^{g}$ require non-trivial term(s) such that $\Phi^{g}$ is
zero in the Higgs symmetric phase at least for $\lambda<\lambda_{c}$.

###### Acknowledgements.

We thank Saumen Datta, Mikko Laine and Kiyoshi Sasaki for helpful discussions.
M. Deka thanks for the hospitality of the Institute of Mathematical Sciences,
India where a part of this work has been completed. All our numerical
computations have been performed at Govorun super-cluster at Joint Institute
for Nuclear Research, Dubna and Annapurna super-cluster based at the Institute
of Mathematical Sciences, India. We have used the MILC collaboration’s public
lattice gauge theory code (version 6) milc  as our base code.

REFERENCES

## References

  * (1) V. A. Kuzmin, V. A. Rubakov and M. E. Shaposhnikov, Phys. Lett. B 155, 36 (1985) doi:10.1016/0370-2693(85)91028-7
  * (2) V. A. Matveev, V. A. Rubakov, A. N. Tavkhelidze and M. E. Shaposhnikov, Sov. Phys. Usp. 31, 916-939 (1988) doi:10.1070/PU1988v031n10ABEH005633
  * (3) V. A. Rubakov and M. E. Shaposhnikov, Usp. Fiz. Nauk 166, 493-537 (1996) doi:10.1070/PU1996v039n05ABEH000145 [arXiv:hep-ph/9603208 [hep-ph]].
  * (4) Z. Fodor, KEK-TH-558.
  * (5) P. H. Damgaard and U. M. Heller, Nucl. Phys. B 294, 253 (1987) doi:10.1016/0550-3213(87)90582-7
  * (6) H. G. Evertz, J. Jersak and K. Kanaya, Nucl. Phys. B 285, 229-252 (1987) doi:10.1016/0550-3213(87)90336-1.
  * (7) K. Farakos, K. Kajantie, K. Rummukainen and M. E. Shaposhnikov, Nucl. Phys. B 442, 317-363 (1995) doi:10.1016/0550-3213(95)80129-4 [arXiv:hep-lat/9412091 [hep-lat]].
  * (8) K. Kajantie, M. Laine, K. Rummukainen and M. E. Shaposhnikov, Nucl. Phys. B 466, 189-258 (1996) doi:10.1016/0550-3213(96)00052-1 [arXiv:hep-lat/9510020 [hep-lat]].
  * (9) K. Rummukainen, M. Tsypin, K. Kajantie, M. Laine and M. E. Shaposhnikov, Nucl. Phys. B 532, 283-314 (1998) doi:10.1016/S0550-3213(98)00494-5 [arXiv:hep-lat/9805013 [hep-lat]].
  * (10) M. Gurtler, E. M. Ilgenfritz and A. Schiller, Phys. Rev. D 56, 3888-3895 (1997) doi:10.1103/PhysRevD.56.3888 [arXiv:hep-lat/9704013 [hep-lat]].
  * (11) B. Bunk, E. M. Ilgenfritz, J. Kripfganz and A. Schiller, Phys. Lett. B 284, 371-376 (1992) doi:10.1016/0370-2693(92)90447-C
  * (12) B. Bunk, E. M. Ilgenfritz, J. Kripfganz and A. Schiller, Nucl. Phys. B 403, 453-474 (1993) doi:10.1016/0550-3213(93)90043-O
  * (13) B. Bunk, Nucl. Phys. Proc. Suppl. 42, 566 (1995). doi:10.1016/0920-5632(95)00313-X.
  * (14) M. Gurtler, E. M. Ilgenfritz, A. Schiller and C. Strecha, Nucl. Phys. B Proc. Suppl. 63, 563-565 (1998) doi:10.1016/S0920-5632(97)00834-7 [arXiv:hep-lat/9709020 [hep-lat]].
  * (15) W. Buchmuller and O. Philipsen, Nucl. Phys. B 443, 47-69 (1995) doi:10.1016/0550-3213(95)00124-B [arXiv:hep-ph/9411334 [hep-ph]].
  * (16) K. Kajantie, M. Laine, K. Rummukainen and M. E. Shaposhnikov, Phys. Rev. Lett. 77, 2887-2890 (1996) doi:10.1103/PhysRevLett.77.2887 [arXiv:hep-ph/9605288 [hep-ph]].
  * (17) Z. Fodor, J. Hein, K. Jansen, A. Jaster, I. Montvay and F. Csikor, Phys. Lett. B 334, 405-411 (1994) doi:10.1016/0370-2693(94)90706-4 [arXiv:hep-lat/9405021 [hep-lat]].
  * (18) Z. Fodor, J. Hein, K. Jansen, A. Jaster and I. Montvay, Nucl. Phys. B 439, 147-186 (1995) doi:10.1016/0550-3213(95)00038-T [arXiv:hep-lat/9409017 [hep-lat]].
  * (19) F. Csikor, Z. Fodor, J. Hein and J. Heitger, Phys. Lett. B 357, 156-162 (1995) doi:10.1016/0370-2693(95)00886-P [arXiv:hep-lat/9506029 [hep-lat]].
  * (20) F. Karsch, T. Neuhaus, A. Patkos and J. Rank, Nucl. Phys. B 474, 217-234 (1996) doi:10.1016/0550-3213(96)00224-6 [arXiv:hep-lat/9603004 [hep-lat]].
  * (21) F. Karsch, T. Neuhaus, A. Patkos and J. Rank, Nucl. Phys. B Proc. Suppl. 53, 623-625 (1997) doi:10.1016/S0920-5632(96)00736-0 [arXiv:hep-lat/9608087 [hep-lat]].
  * (22) Y. Aoki, Phys. Rev. D 56, 3860 (1997). doi:10.1103/PhysRevD.56.3860 [hep-lat/9612023].
  * (23) Y. Aoki, F. Csikor, Z. Fodor and A. Ukawa, Phys. Rev. D 60, 013001 (1999). doi:10.1103/PhysRevD.60.013001 [hep-lat/9901021].
  * (24) C. Bonati, G. Cossu, M. D’Elia and A. Di Giacomo, Nucl. Phys. B 828, 390-403 (2010) doi:10.1016/j.nuclphysb.2009.12.003 [arXiv:0911.1721 [hep-lat]].
  * (25) M. D’Onofrio, K. Rummukainen and A. Tranberg, Phys. Rev. Lett. 113, no.14, 141602 (2014) doi:10.1103/PhysRevLett.113.141602 [arXiv:1404.3565 [hep-ph]].
  * (26) M. D’Onofrio and K. Rummukainen, Phys. Rev. D 93, no.2, 025003 (2016) doi:10.1103/PhysRevD.93.025003 [arXiv:1508.07161 [hep-ph]].
  * (27) M. D’Onofrio, K. Rummukainen and A. Tranberg, PoS LATTICE2012, 055 (2012) doi:10.22323/1.164.0055 [arXiv:1212.3206 [hep-ph]].
  * (28) M. Laine, G. Nardini and K. Rummukainen, JCAP 01, 011 (2013) doi:10.1088/1475-7516/2013/01/011 [arXiv:1211.7344 [hep-ph]].
  * (29) O. Gould, S. Güyer and K. Rummukainen, [arXiv:2205.07238 [hep-lat]].
  * (30) K. Kajantie, Nucl. Phys. B Proc. Suppl. 42, 103-112 (1995) doi:10.1016/0920-5632(95)00192-C [arXiv:hep-lat/9412072 [hep-lat]].
  * (31) K. Jansen, Nucl. Phys. B Proc. Suppl. 47, 196-211 (1996) doi:10.1016/0920-5632(96)00045-X [arXiv:hep-lat/9509018 [hep-lat]].
  * (32) K. Rummukainen, Nucl. Phys. B Proc. Suppl. 53, 30-42 (1997) doi:10.1016/S0920-5632(96)00597-X [arXiv:hep-lat/9608079 [hep-lat]].
  * (33) R. D. Pisarski and F. Wilczek, Phys. Rev. D 29, 338-341 (1984) doi:10.1103/PhysRevD.29.338.
  * (34) F. Karsch and S. Stickan, Phys. Lett. B 488, 319-325 (2000) doi:10.1016/S0370-2693(00)00902-3 [arXiv:hep-lat/0007019 [hep-lat]].
  * (35) J. L. Alonso, V. Azcoiti, I. Campos, J. C. Ciria, A. Cruz, D. Iniguez, F. Lesmes, C. Piedrafita, A. Rivero and A. Tarancon, et al. Nucl. Phys. B 405, 574-592 (1993) doi:10.1016/0550-3213(93)90560-C [arXiv:hep-lat/9210014 [hep-lat]].
  * (36) W. Janke and R. Villanova, Nucl. Phys. B 489, 679-696 (1997) doi:10.1016/S0550-3213(96)00710-9 [arXiv:hep-lat/9612008 [hep-lat]].
  * (37) N. B. Wilding, Phys. Rev. E 52, 602-611 (1995) doi:10.1103/PhysRevE.52.602 [arXiv:cond-mat/9503145 [cond-mat]].
  * (38) N. B. Wilding and M. Mueller, J. Chem. Phys. 102, 2562 (1995) doi:10.1063/1.468686 [arXiv:cond-mat/9410077 [cond-mat]].
  * (39) Nigel B Wilding 1997 J. Phys.: Condens. Matter 9 585.
  * (40) J. J. Rehr and N. D. Mermin, Phys. Rev. A 8, 472-480 (1973) doi:10.1103/PhysRevA.8.472.
  * (41) K. Kanaya and S. Kaya, Phys. Rev. D 51, 2404-2410 (1995) doi:10.1103/PhysRevD.51.2404 [arXiv:hep-lat/9409001 [hep-lat]].
  * (42) D. Toussaint, Phys. Rev. D 55, 362-366 (1997) doi:10.1103/PhysRevD.55.362 [arXiv:hep-lat/9607084 [hep-lat]].
  * (43) http://www.physics.utah.edu/~detar/milc/.
  * (44) M. Creutz, Phys. Rev. D 21, 2308 (1980). doi:10.1103/PhysRevD.21.2308.
  * (45) N. Cabibbo and E. Marinari, Phys. Lett. B 119, 387 (1982). doi:10.1016/0370-2693(82)90696-7.
  * (46) C. Whitmer, Phys. Rev. D 29, 306 (1984). doi:10.1103/PhysRevD.29.306.
  * (47) S. Aoki et al. [CP-PACS Collaboration], Phys. Rev. D 70, 034503 (2004). doi:10.1103/PhysRevD.70.034503 [hep-lat/0312011].
  * (48) M. Wurtz, R. Lewis and T. G. Steele, Phys. Rev. D 79, 074501 (2009) doi:10.1103/PhysRevD.79.074501 [arXiv:0902.1167 [hep-lat]].
  * (49) M. Wurtz and R. Lewis, Phys. Rev. D 88, 054510 (2013) doi:10.1103/PhysRevD.88.054510 [arXiv:1307.1492 [hep-lat]].
  * (50) Alan M. Ferrenberg and D. P. Landau, Phys. Rev. B 44, 5081, 1991.
  * (51) M. Hasenbusch, K. Pinn and S. Vinti, Phys. Rev. B 59, 11471-11483 (1999) doi:10.1103/PhysRevB.59.11471 [arXiv:hep-lat/9806012 [hep-lat]].


# Cost Management on Commercial Cloud Platforms

G. Bruce Berriman1, William O’Mullane2, Arik Mitschang3, and Ivelina Momcheva
4

###### Abstract

Commercial cloud platforms are a powerful technology for astronomical
research. Despite the benefits of cloud computing—such as on-demand
scalability and reduction of systems management overhead—confusion over how to
manage costs remains, for many, one of the biggest barriers to entry. This
confusion is exacerbated by the rapid growth in services offered by commercial
providers, by the growth in the number of these providers, and by storage,
compute, and I/O metered at separate rates—all of which can change without
notice. As a rule, processing is very cheap, storage is more expensive, and
downloading is very expensive. Thus, an application that produces large image
data sets for download will be far more expensive than an application that
performs extensive processing on a small data set. This Birds of a Feather
(BoF) session aimed to quantify the above statement by presenting case studies
of the costing of astronomy applications on commercial clouds that covered a
range of processing scenarios; these presentations were the basis for
discussion by the attendees.

1Caltech/IPAC-NExScI, Pasadena, CA 91125, USA<EMAIL_ADDRESS>
2Rubin Observatory Project Management Office, Tucson, AZ 85721, USA

3Johns Hopkins University, Baltimore, MD 21218, USA

4 Space Telescope Science Institute, Baltimore, MD 21218, USA

##  1 Introduction: Why Is Cost Management on the Cloud Difficult?

Cost management of cloud services can be complex and difficult because pricing
is variable, with separate rates for compute, storage, data movement and use
of APIs; because pricing changes without notice; because an abundance of
providers now offer an exceptionally diverse range of services; and because
billing practices can appear opaque. The authors of this paper have, however,
acquired considerable collective expertise in managing costs. This paper aims
to pass on this expertise for the benefit of the broader community, and
represents the findings of a Birds of a Feather (BoF) session held as part of
the ADASS XXX virtual conference.

The BoF presented five case studies of the costing of astronomy projects on
commercial clouds:

  * •

W. O’Mullane "Hosting the Rubin Observatory Interim Data Facility on a cloud
platform" (hereafter, WOM)

  * •

A. Mitschang "SciServer on a cloud platform" (hereafter, AM)

  * •

I. Momcheva "Managing cloud services at STScI" (hereafter, IMO)

  * •

B. Berriman "Creating TESS Mosaics on AWS" (hereafter, BBT)

  * •

B. Berriman "Cost Management Workshop at IPAC" (hereafter, BBC)

These presentations were followed by a discussion of cost management among the
100 persons attending. What follows summarizes the conclusions of the
presentations and of the discussion. All prices quoted below refer to prices
in effect at the dates of provisioning. There is no guarantee that such prices
remain in effect; nor does any pricing quoted constitute an endorsement of a
provider.

##  2 How Much Will It All Cost: Getting Started

So far, at least in the US, funding agencies for astronomy and space science
have not completed negotiations with cloud providers to offer special pricing
for their grantees. At present, individual institutions or projects must
negotiate themselves (WOM, IMO). They may also be able to take advantage of
third-party "holding organizations" such as Cloud Bank. This pilot project
manages cloud accounts for grantees for some NSF awards to avoid the high
overheads imposed by many US universities. Such overheads are generally not
imposed on purchases of on-premises hardware and can make total cloud costs
prohibitively high.

It is important to understand what compute, storage, backup, database, and API
or toolkit resources an application requires during the period it will run on
the cloud. While many resources such as computing are inexpensive, others
resources such as Amazon Elastic File Storage and egress of very large volumes
out of the cloud can prove expensive (BBT, BBC). So can moving data across
regional cloud data centers if held in different accounts, but this can be
mitigated by using a third-party to manage them. Services such as Amazon
Virtual Private Cloud (VPC) endpoints connect to the cloud without an internet
gateway and can considerably reduce data movement costs.

The type of application needs to be considered as well. Cloud platforms are
built with commodity hardware and are most suitable for embarrassingly
parallel applications, rather than HPC or HTC applications, which require
high-speed network connections. Moreover, be aware that experimenting with
data types not yet well supported can lead to considerable expense; the
example of memory mapping n-dimensional files is a case in point discussed
during the BoF.

##  3 Enterprise Price Negotiations

It is crucial to perform due diligence and estimate the costs for all aspects
of a data center. BBC described the need to consider overhead costs,
availability marginal cost, sizing and flexibility, data transfer (ingress and
egress), networking and firewalls, monitoring, and systems administration.
There is no one-size-fits-all solution. For the case of the IPAC data center,
with over 1 PB, Amazon Web Services (AWS) proved superior for availability,
for sizing and flexibility, and for networks and firewalls, but IPAC was
cheaper for overheads, marginal cost, and data egress; monitoring and systems
administration were a tie.

In spring 2020, Wikipedia reported there were 285 cloud platforms offering a
wealth of services, each with their own pricing. Some, such as the "big
three"—AWS, Microsoft Azure and Google Cloud—can be considered "full service
providers," while others may specialize in niche services. If, however, all
that is needed is off-premises storage, a dedicated cloud storage provider may
be the best choice; Wikipedia listed 143 of these in spring 2020. BBC pointed
out that, for example, Wasabi offers "S3"-type object storage at a rate much
cheaper than AWS, without egress charges, but it does not offer compute
services. BBC also pointed out the cost of back-up storage for IPAC’s data on
AWS Glacier was comparable to on-premises backup, but may prove cheaper in
total because on-premises maintenance costs are eliminated.

The cost estimates, once complete, provide a basis for negotiations with
providers. This is best done by a dedicated cost management team familiar with
the technical needs as well as business and financial management, and by
engaging directly with providers.

WOM and IMO presented two examples of the results of negotiations conducted by
the Rubin Observatory and the Space Telescope Science Institute (STScI). WOM
described how, after thorough prototyping experiments with Google Cloud and
AWS, the Rubin Observatory set up an interim data facility with Google Cloud
as the provider for three years.

IMO described how STScI used AWS through a third-party provider until 2018,
and then engaged in major negotiations to establish an enterprise agreement
with AWS. The best practices that emerged in negotiating discounts included:
an upfront payment discount (“savings plan”) of 50-60% (best for predictable
workloads and a plan for storage costs upfront with commitment for a 10-20%
discount. The option to add credits can be added to discounts and options to
waive some egress costs—typically 15% of annual spending (this is now provided
as standard by Google to Internet2111https://internet2.edu/ members) can be
negotiated.

##  4 What of Lonely Astronomers?

Astronomers without access to negotiated prices have opportunities to avoid
paying retail rates. IMO and BBT described these practices in detail:

  * •

Do an on-premises pilot project if possible, to understand the data and
resources needed (BBT).

  * •

Understand workflow needs and how they translate to cloud services. How long
will you need cloud services? Do you need to reprocess periodically? The TESS
example demonstrated the need for large storage volume of 150 TB over the life
of the project.

  * •

Use free or low-cost tiers to do the set-up.

  * •

Use the AWS Command Line Interface (CLI) to investigate compute resource
management before, for example, building compute clusters.

  * •

Consider requesting research credits to defray part of the costs.

  * •

Consider using public data sets programs, which provide public serving of data
of wide interest (though not for private data).

  * •

Use "spot" pricing where possible, which offers very cheap rates for unused
capacity. Jobs may stop if demand rises; jobs are saved and can be restarted.

  * •

If there is adequate demand, use reserved server instances. AM showed that for
the SciServer science platform, this can reduce annual charges by over 45%
compared with on-demand instances. item Consider infrastructure-as-code tools
(such as Terraform) to simplify creation and destruction of cloud resources.

  * •

Use serverless mode, whereby users submit their code and the provider runs and
optimizes the server. Costs only accrue when jobs are running.

  * •

Use S3 object storage whenever possible. Allow enough time to overcome the
latency in S3 storage. Explore different S3 tiers for data that are not
accessed frequently.

  * •

Learn how to use cost calculators offered by providers. With some practice,
they yield estimates that are good to 80%.

  * •

Ask support for help provisioning the resources you need.

A number of attendees proposed holding cost management workshops at national
astronomy conferences. A repository of useful cost management tools and
documentation was also thought valuable; an organization such as Cloud Bank
may ultimately provide this function.

A further option is to use national or private clouds, though they may not be
as powerful or full featured as commercial clouds. Examples are European Open
Science Cloud (Europe), Nectar (Australia), and TACC Jetstream (US). They are
free to qualified users, at the expense of writing a proposal.

##  5 Monitoring Costs and Resource usage: Avoiding Fiscal Black Holes

Diligent attention to several procedures can avoid unnecessary costs and
prevent financial headaches. IMO described them:

  * •

Do monitor costs daily.

  * •

Do clean up jobs and do not forget about them.

  * •

Do run a test/prototype first and examine the changes in costs before scaling
up.

  * •

Do understand the long tail of costs; avoid having jobs run for a long time.

  * •

Do use standard instances before customizing; customize or scale up only when
necessary.

  * •

Do think carefully about the storage costs.

  * •

Do download only the data you need. If you are downloading everything, re-
think your cloud approach. Downloading large volumes of data can be very
expensive (BBT, BBC).

  * •

Do set up all the billing alarms.

  * •

Do automate shut down of a workflow once a spending cap is reached.


# Concurrent consideration of cortical and cancellous bone within continuum
bone remodelling

Ina Schmidta, Areti Papastavroua, and Paul Steinmannb CONTACT Ina Schmidt
Author. Email<EMAIL_ADDRESS>aFaculty of Mechanical Engineering,
Nuremberg Tech, Keßlerplatz 12, 90489 Nuremberg, Germany; bChair of Applied
Mechanics, University of Erlangen-Nuremberg, Egerlandstraße 5, 91058 Erlangen,
Germany

###### Abstract

Continuum bone remodelling is an important tool for predicting the effects of
mechanical stimuli on bone density evolution. While the modelling of only
cancellous bone is considered in many studies based on continuum bone
remodelling, this work presents an approach of modelling also cortical bone
and the interaction of both bone types. The distinction between bone types is
made by introducing an initial volume fraction. A simple point-wise example is
used to study the behaviour of novel model options, as well as a proximal
femur example, where the interaction of both bone types is demonstrated using
initial density distributions. The results of the proposed model options
indicate that the consideration of cortical bone remarkably changes the
density evolution of cancellous bone, and should therefore not be neglected.

###### keywords:

bone remodelling; cortical bone; cancellous bone; density change; finite
element method

##  1 Introduction

Tubular bone is characterised by cortical bone surrounding the bone marrow in
the shaft and the cancellous bone at the proximal and distal ends. Due to its
special microstructure composed of beams and plates (trabeculae), cancellous
bone is also denoted as spongy bone. Depending on environmental changes, such
as mechanical loading, bone can adapt its internal microstructure by altering
the density and trabecular orientation. Remodelling of cancellous bone was
first formulated by (Wolff 1892) and is referred to as Wolff’s law of bone
remodelling. Based on continuum mechanics and the finite element method,
simulation approaches to bone remodelling were developed, see (Cowin and
Hegedus 1976), (Huiskes et al. 1987), (Weinans et al. 1992), (Harrigan and
Hamilton 1993), (Jacobs et al. 1995) and many more. Most of these studies
address exclusively the adaptation of cancellous bone. Even if the focus is
also on the adaptation of cortical bone, it is typically considered separately
from spongy bone, see (Carter and Beauprè 2010). One of the few approaches
involving the interaction of both bone types was investigated by (Hambli 2014)
and (Barkaoui et al. 2019) based on the cellular activities of bone. Another
proposal was developed by (Nutu 2018) based on the model of (Huiskes et al.
1987) for dental implants.

This work aims to simulate bone density changes of cancellous and cortical
bone as well as their interplay based on a continuum bone remodelling
approach. Unlike spongy bone, the bone fibrils of cortical bone are densely
packed, which is why this lamellar structure is also called compact bone.
Remodelling of cortical bone can occur along the endosteal and periosteal
surface as well as within the channels. The surface area per unit volume is
considerably reduced compared to cancellous bone and therefore remodelling
requires more time with less change in density. In order to adapt to the load
case, the main mechanism considered here is the adjustment of the
microstructural density. Whereas a change in bone geometry will not be
considered, the simulation of bone adjustments in the endosteal area such as
trabecularisation and thinning of the inner compact bone in old age is made
possible.

Based on previous studies, see (Kuhl and Steinmann 2003), (Liedtke et al.
2017) (Papastavrou et al. 2020a) and (Papastavrou et al. 2020b), an initial
volume fraction is introduced in this work, which is defined as the ratio
between initial nominal density and true density and therefore allows a
distinction between the two bone types. The influence of this new parameter
and the associated modelling of spongy and compact bone will be studied in
this context.

In the following, the quantities $\\{\circ\\}$, which take into account the
effective mixture of solid bone material and the (empty) pore space are
referred to as nominal quantities, whereas true quantities based on the
properties of the solid (bone) material are denoted as $\\{\circ\\}^{s}$.
Furthermore we distinguish between the initial value $\\{\circ\\}^{\star}$ of
a quantity (at time zero), and its current (time-evolving) value
$\\{\circ\\}(t)$. The latter should not be confused with the actual
configuration, which refers to the geometric description of the deformed
continuum body.

This contribution is organised as follows. The governing and constitutive
equations including the proposed modification are detailed in Section 2. The
performance is illustrated by numerical problems in Section 3. Section 4
provides a concluding discussion.

##  2 Continuum framework of bone remodelling

With a view on its continuum modelling, bone is here considered as a porous
material consisting of a solid skeleton and (empty) pore space, see Fig. 1. In
the following, the underlying equations of continuum bone remodelling are
briefly re-iterated, see also (Kuhl and Steinmann 2003) and (Holzapfel 2001)
for additional details.

###  2.1 Kinematics

The kinematic description is characterised by the deformation map
$\boldsymbol{\varphi}$ that maps the referential placement $\mathbf{X}$ of a
physical particle in the material (reference) configuration $\mathcal{B}_{0}$
at time $t_{0}$ to its actual placement $\mathbf{x}$ in the spatial (actual)
configuration $\mathcal{B}_{t}\subset\mathbb{E}^{3}$ at time
$t\in\mathbb{R}_{+}$

$\mathbf{x}=\boldsymbol{\varphi}(\mathbf{X},t):\mathcal{B}_{0}\times\mathbb{R}_{+}\to\mathcal{B}_{t}\>.$
(1)

The referential gradient $\nabla_{0}$ of the nonlinear deformation map denotes
the deformation gradient $\mathbf{F}$, it maps from the material tangent space
$T\mathcal{B}_{0}$ to the spatial tangent space $T\mathcal{B}_{t}$

$\mathbf{F}=\nabla_{0}\boldsymbol{\varphi}(\mathbf{X},t):T\mathcal{B}_{0}\to
T\mathcal{B}_{t}\>.$ (2)

The determinant of the deformation gradient is the Jacobian
$J=\det\mathbf{F}>0$ and relates the volume elements between the two
configurations. The right and left Cauchy-Green tensor
$\mathbf{C}=\mathbf{F}^{T}\cdot\mathbf{F}$ and
$\mathbf{b}=\mathbf{F}\cdot\mathbf{F}^{T}$ are introduced as characteristic
strain measures.

In the following, the expression
$\dot{\overline{\\{\circ\\}}}=\partial_{t}\\{\circ\\}|_{X}$ is used for the
material time derivative of a quantity $\\{\circ\\}$ at fixed material
placement $\mathbf{X}$. The presented formulation refers completely to the
material configuration. Therefore, $\nabla_{0}\\{\circ\\}$ and
$\mbox{Div}\\{\circ\\}$ denote the (referential) gradient and divergence of
any field $\\{\circ\\}$ with respect to the material placement $\mathbf{X}$. A
density quantity $\\{\circ\\}$ may be expressed per unit volume
$\\{\circ\\}_{0}$ or per unit mass $\\{\circ\\}$, both related via the current
referential nominal density $\rho_{0}$ as
$\\{\circ\\}_{0}=\rho_{0}\\{\circ\\}$.

###  2.2 Balance equations

The balance of linear momentum is given by

$\mbox{Div}\mathbf{P}=\mathbf{0}$ (3)

with $\mathbf{P}=\mathbf{P}(\nabla_{0}\boldsymbol{\varphi},\rho_{0},...)$ the
Piola stress tensor. Since the time scale of inertial effects is significantly
smaller than that of the remodelling process, the inertial term is neglected
herein (Frost 1987). The loads applied during locomotion typically correspond
to a multiple of gravity. Therefore, body forces are usually not considered in
bone remodelling as they are not the dominant stimulus (Jacobs et al. 1995).

The balance of mass is given by

$\dot{\rho_{0}}=R_{0}$ (4)

where $R_{0}$ is the scalar-value mass source term. As the actual nominal
density $\rho_{t}={\rho_{0}}/{J}$ is time-dependent also due to the coupling
with the Jacobian, it is useful to capture the evolution of the referential
nominal density $\rho_{0}$. The relative change in (referential) nominal
density $\tilde{\rho_{0}}$ is defined in terms of the initial nominal density
$\rho_{0}^{\star}$ as

$\tilde{\rho_{0}}=\frac{\rho_{0}-\rho_{0}^{\star}}{\rho_{0}^{\star}}.$ (5)

Therefore, density gain results for values $\tilde{\rho_{0}}>0$ whereas
density loss is obtained for $-1<\tilde{\rho_{0}}<0$.

Figure 1: Homogenised macro scale with nominal density $\rho_{0}(t)$ versus
heterogeneous micro scale with true density $\rho_{0}^{s}$

###  2.3 Nominal versus true mass densities

We distinguish between the homogenised macro scale with nominal (time-
evolving, referential) density $\rho_{0}=\rho_{0}(t)\in(0,\rho_{0}^{s}]$
occupying the (time-independent, referential) volume element $\,\mbox{d}V$ and
the heterogeneous micro scale with true (time-independent, referential)
density $\rho_{0}^{s}$ of the solid skeleton material occupying the (time-
evolving, referential) partial volume element
$\,\mbox{d}V^{s}=\,\mbox{d}V^{s}(t)\in(0,\,\mbox{d}V]$, see Fig. 1. The (time-
evolving, referential) volume fraction of solid skeleton material thus follows
as

$\nu_{0}=\nu_{0}(t):=\frac{\,\mbox{d}V^{s}(t)}{\,\mbox{d}V}=\frac{\rho_{0}(t)}{\rho_{0}^{s}}\in(0,1].$
(6)

The (referential) volume fraction $\nu_{0}(t)$ indicates how much space within
a (referential) volume element $\,\mbox{d}V$ at the macro scale is occupied by
solid skeleton material, thus $1-\nu_{0}(t)$ denotes the (referential) volume
fraction of the (empty) pore space.

At time $t=0$ the corresponding initial values for the (referential) volume
fraction $\nu_{0}(0)$, the (referential) partial volume element
$\,\mbox{d}V^{s}(0)$ and the nominal (referential) density $\rho_{0}(0)$ are
related via

$\nu_{0}^{\star}:=\frac{\,\mbox{d}V^{s\star}}{\,\mbox{d}V}=\frac{\rho_{0}^{\star}}{\rho_{0}^{s}}\in(0,1].$
(7)

and are denoted the initial (referential) volume fraction
$\nu_{0}^{\star}:=\nu_{0}(0)$, the initial (referential) partial volume
$\,\mbox{d}V^{s\star}:=\,\mbox{d}V^{s}(0)$ and the initial nominal
(referential) density $\rho_{0}^{\star}:=\rho_{0}(0)$, respectively.
Typically, $\nu_{0}^{\star}$ is significantly higher for cortical bone as
compared to cancellous bone.

Furthermore, we may distinguish between the nominal relative (referential)
density and its true counterpart

$\frac{\rho_{0}(t)}{\rho_{0}^{\star}}\equiv\frac{\nu_{0}(t)}{\nu_{0}^{\star}}\in(0,\infty)\quad\mbox{versus}\quad\frac{\rho_{0}(t)}{\rho_{0}^{s}}\equiv\nu_{0}(t)\in(0,1]$
(8)

when phenomenologically modelling the homogenised mechanical properties of
bone based on its skeleton structure at the micro scale, i.e. its
microstructure. In the following, the present contribution aims to rationalise
various modelling options stemming from this distinction.

###  2.4 Nominal versus true free energy densities

Bone possesses a microstructure characterised by the spatially heterogeneous
distribution of solid skeleton and pore space. The resulting structural bone
properties of the heterogeneous solid skeleton at the micro scale may then be
homogenised into nominal bone properties of the effective material at the
macro scale.

A simplistic, however efficient and effective homogenisation rule – typically
employed for cellular materials – relates the true material parameters of the
solid skeleton at the micro scale to the corresponding effective, nominal
parameters of the homogenised material at the macro scale via a power $N$ of
the initial (referential) volume fraction $\nu_{0}^{\star}$. Let for example
$\lambda^{s}$ and $\mu^{s}$ denote the Lamé parameters of the solid skeleton
material. Then, the homogenised, effective bone material is characterised by
the nominal material parameters

$\lambda^{\star}=[\nu_{0}^{\star}]^{N}\lambda^{s}\quad\mbox{and}\quad\mu^{\star}=[\nu_{0}^{\star}]^{N}\mu^{s}\quad\mbox{with}\quad
N=\left\\{1,n\right\\}\>.$ (9)

Typically, the concrete values assumed for the homogenisation exponent $N$
capture different microstructures and homogenisation approaches. In this
contribution, either the most elementary volumetric mixture rule with $N=1$ or
a relation with $N=n$ shall be considered. The latter is particularly
motivated by the proposal of (Gibson 2005) postulating a dependence of the
nominal Young’s modulus of a homogenised (linear elastic) porous material on
the Young’s modulus of the solid skeleton material via a power of the true
relative density.

In what follows we assume the concept of strain equivalence as also employed
in damage mechanics. To elucidate this aspect more concretely, we employ,
without loss of generality, a (referential) volume-specific Neo-Hookean true
free energy density for the solid skeleton material

$\Psi_{0}^{s}=\frac{1}{2}\lambda^{s}\ln^{2}J+\frac{1}{2}\mu^{s}\left[\mathbf{F}:\mathbf{F}-3-2\ln
J\right].$ (10)

Note that strain equivalence postulates identical ’strain’ at the macro and
the micro scale. Thus, based on the nominal material parameters, the
corresponding (referential) volume- specific initial nominal free energy
follows from its true counterpart as

$\Psi_{0}^{\star}=[\nu_{0}^{\star}]^{N}\Psi_{0}^{s}\>.$ (11)

Concretely, the (referential) volume-specific Neo-Hookean initial nominal free
energy reads in terms of the nominal Lamé parameters as

$\Psi_{0}^{\star}=\frac{1}{2}\lambda^{\star}\ln^{2}J+\frac{1}{2}\mu^{\star}\left[\mathbf{F}:\mathbf{F}-3-2\ln
J\right]\>.$ (12)

The initial microstructure evolves with time due to remodelling into the
current microstructure, thus necessitating to express the (referential)
volume-specific current nominal free energy density $\Psi_{0}$ in terms of its
initial counterpart.

To this end, again motivated by (Gibson 2005), when taking the microstructure
of cellular materials into account, the (referential) volume-specific current
nominal free energy density $\Psi_{0}$ is determined by weighting its initial
nominal counterpart $\Psi_{0}^{\star}$ with a power $n$ of the nominal
relative density

$\Psi_{0}=\left[\frac{\rho_{0}}{\rho_{0}^{\star}}\right]^{n}\Psi_{0}^{\star}\>.$
(13)

Then the Piola stress entering the equilibrium equation derives from the
(referential) volume-specific current nominal free energy density $\Psi_{0}$
as

$\mathbf{P}:=\frac{\partial\Psi_{0}}{\partial\mathbf{F}}=\left[\frac{\rho_{0}}{\rho_{0}^{\star}}\right]^{n}\frac{\partial\Psi_{0}^{\star}}{\partial\mathbf{F}}=:\left[\frac{\rho_{0}}{\rho_{0}^{\star}}\right]^{n}\mathbf{P}^{\star}$
(14)

and expands concretely for the Neo-Hookean case with the initial Piola stress

$\mathbf{P}^{\star}=\lambda^{\star}\ln
J\,\mathbf{F}^{-t}+\mu^{\star}[\mathbf{F}-\mathbf{F}^{-t}]\>.$ (15)

Incorporating the (referential) volume-specific Neo-Hookean true free energy
density $\Psi_{0}^{s}$ eventually renders $\Psi_{0}$ in terms of the power $n$
of the true relative density and the power $N-n$ of the initial (referential)
volume fraction

$\Psi_{0}=\left[\frac{\rho_{0}}{\rho_{0}^{s}}\right]^{n}[\nu_{0}^{\star}]^{N-n}\,\Psi_{0}^{s}=[\nu_{0}]^{n}[\nu_{0}^{\star}]^{N-n}\,\Psi_{0}^{s}\>.$
(16)

The exponent $n$ is typically set to $n=2$ for trabecular bone with low
nominal density and to $n>2$ for denser cancellous bone, compare to (Rho et
al. 1993) and (Gibson 2005).

Setting finally the homogenisation exponent $N$ in particular to $N=n$ results
in an expression for the (referential) volume-specific current nominal free
energy density $\Psi_{0}$ in terms of its true counterpart and the true
relative density

$\Psi_{0}\stackrel{{\scriptstyle
N=n}}{{=}}\left[\frac{\rho_{0}}{\rho_{0}^{s}}\right]^{n}\Psi_{0}^{s}=[\nu_{0}]^{n}\Psi_{0}^{s}\>.$
(17)

Correspondingly, the Piola stress expands in this case as

$\mathbf{P}\stackrel{{\scriptstyle
N=n}}{{=}}\left[\frac{\rho_{0}}{\rho_{0}^{s}}\right]^{n}\frac{\partial\Psi_{0}^{s}}{\partial\mathbf{F}}=:\left[\frac{\rho_{0}}{\rho_{0}^{s}}\right]^{n}\big{[}\lambda^{s}\ln
J\,\mathbf{F}^{-t}+\mu^{s}[\mathbf{F}-\mathbf{F}^{-t}]\big{]}\>.$ (18)

Setting $N=n$ tacitly assumes that the geometry of the microstructure is not
distorted by the deformation (for example a circular pore remains a circular
pore). In contrast, setting $N\not=n$ conceptually allows to distinguish
between the initial and the (deformation-dependent) current microstructure,
thus in general $n=n(\mathbf{F})$ with $N=n(\mathbf{I})$.

This contribution shall explore the consequences of two limiting cases by
setting $N$ to either $N=1$ (volumetric mixture rule) or to $N=n$ (cellular
material like homogenisation rule). More sophisticated scenarios with
$n=n(\mathbf{F})$ where $N=n(\mathbf{I})$ shall not be analysed here.

###  2.5 Mass density flux

To model changes in the nominal (referential) density $\rho_{0}$ incurred by
mechanical stimuli, the volume-localised mass balance $\dot{\rho}_{0}=R_{0}$
is equipped with a (referential) source term $R_{0}$. Following (Harrigan and
Hamilton 1993), the (referential) mass source $R_{0}$ is expressed by the
(referential) volume-specific current nominal free energy density $\Psi_{0}$
and the nominal relative density as

$R_{0}=c\left[\left[\frac{\rho_{0}}{\rho_{0}^{\star}}\right]^{-m}\,\Psi_{0}-\Psi_{0}^{a}\right]\>.$
(19)

Here, the parameter $c$ controls the speed of density evolution and has
dimension of time divided by length squared. The attractor stimulus
$\Psi_{0}^{a}$ defines the biological equilibrium (homeostasis) to which the
system is driving at constant (mechanical) stimulus. The exponent $m$ is a
material parameter particularising the mechanical stimulus (Harrigan and
Hamilton 1993) in terms of the nominal relative density.

The (referential) mass source may alternatively be expressed by the
(referential) volume-specific initial nominal free energy density
$\Psi_{0}^{\star}$ or, likewise, by its true counterpart $\Psi_{0}^{s}$ as

$R_{0}=c\left[\left[\frac{\rho_{0}}{\rho_{0}^{\star}}\right]^{n-m}\,\Psi_{0}^{\star}-\Psi_{0}^{a}\right]=c\left[\left[\frac{\rho_{0}}{\rho_{0}^{s}}\right]^{n-m}[\nu_{0}^{\star}]^{N-n+m}\,\Psi_{0}^{s}-\Psi_{0}^{a}\right]\>.$
(20)

We may deduce from this representation that the exponent $m$ is also crucial
for the stability of the density evolution: at homeostasis ($R_{0}=0$)
$\rho_{0}>\rho_{0}^{\star}$ for $\Psi_{0}^{\star}>\Psi_{0}^{a}$ only if
$m-n>0$. In conclusion, the exponent $m$ should be chosen as $m>n$ (Harrigan
and Hamilton 1994), (Kuhl and Steinmann 2003).

Setting finally the homogenisation exponent $N$ in particular to $N=n$ results
in alternative expressions for the (referential) mass source

$R_{0}=c\left[\left[\frac{\rho_{0}}{\rho_{0}^{s}}\right]^{n-m}[\nu_{0}^{\star}]^{m}\,\Psi_{0}^{s}-\Psi_{0}^{a}\right]=c\left[\left[\frac{\rho_{0}}{\rho_{0}^{\star}}\right]^{n-m}[\nu_{0}^{\star}]^{n}\,\Psi_{0}^{s}-\Psi_{0}^{a}\right]\>.$
(21)

##  3 Simulation and results

The performance of the proposed model options ($N=n$ versus $N=1$) is
demonstrated using classical benchmark problems. The first example considers a
unit square specimen as defined by (Kuhl and Steinmann 2003) to study temporal
bone density evolution at a single point. Second, we analyse the consequence
of the spatial distribution of bone microstructure using a two-dimensional
proximal femur head as given by (Carter and Beauprè 2010).

###  3.1 Investigating the influence of the initial volume fraction

To differentiate between cancellous and cortical bone tissue, the effect of
the initial volume fraction shall be examined. For this purpose, a simple unit
square sample is studied, which is loaded in tension by a force function that
increases stepwise as illustrated in Figure 2.

Figure 2: Unit square example: specimen and load function

The reference parameters are chosen similar to (Kuhl and Steinmann 2003): true
elastic modulus $E^{s}=1.0$, Poisson’s ratio $\nu^{s}=0.0$, attractor stimulus
$\Psi_{0}^{a}=1.0$, parameter $m=3$ and parameter $n=2$. The homogenisation
exponent is chosen as $N=n$ as motivated by the proposal of (Gibson 2005).

A true density of the solid material $\rho_{0}^{s}=2$, similar to (Gibson
2005), is assumed for both bone microstructure types in the following. It is
known from measurements that the true density of trabeculae is slightly
smaller and the one of cortical bone fibrils slightly larger, see (An and
Draughn 2000). However, as the objective of this analysis is to point out the
general trend in the evolution of the different tissue microstructures, the
assumption of equal true density is sufficient in this context. An initial
nominal density of $\rho_{0}^{\star}=0.5$ leading to an initial volume
fraction of $\nu_{0}^{\star}={\rho_{0}^{\star}}/{\rho_{0}^{s}}={0.5}/{2}=0.25$
is chosen as representative of cancellous bone. A comparatively high value of
$\nu_{0}^{\star}={\rho_{0}^{\star}}/{\rho_{0}^{s}}={1.8}/{2}=0.9$ is set for
the volume fraction representing compact bone.

Figure 3: Unit square example: evolution of relative change in nominal density
and absolute nominal density with a variation in the initial volume fraction
representing cortical and spongy bone

As shown in Figure 3, a small initial volume fraction (gray line) leads to
rapid bone growth while a high value (black line) causes no density change in
the first load step. Only after increasing the load in the second step, and
thus increasing the mechanical stimulus, an adaptation of bone density can be
observed.

Cortical bone in the model therefore requires a significantly higher stimulus
and grows less in relative terms. But despite its rapid growth, the absolute
density of cancellous bone still remains at a lower level.

###  3.2 Investigating the effect of exponent $m$ and the attractor stimulus

Due to its small surface area in relation to the volume, cortical bone adapts
to the load conditions much slower and less than spongy bone. Considering the
mass source term in Equation 19, it appears that there are three parameters
which can influence the behaviour of bone adaptation. By reducing the
parameter $c$, the speed of density evolution can be slowed down, but after a
sufficiently long time the result always converges towards the same saturation
level. Therefore, a parameter adjustment of the exponent $m$ and the attractor
stimulus $\Psi_{0}^{a}$ shall be considered for cortical bone with
$\nu_{0}^{\star}=0.9$ and $N=n=2$.

In order to satisfy the stability criterion, the exponent $m$ must be selected
greater than $n$, see (Harrigan and Hamilton 1993). Thus, $m=3$, $m=2.5$ and
$m=10$ are investigated as depicted in Figure 4.

Figure 4: Unit square example: evolution of relative change in nominal density
and absolute nominal density with variations of $m$ for cortical bone

As can be seen, reducing the exponent leads to an increase in density growth
(dashed line). In contrast, choosing a large value for the exponent $m$ tends
to slow down the density evolution (dotted line).

One possibility to simulate the significantly reduced density evolution of
compact bone compared to spongy bone would therefore be setting a relatively
high exponent $m$.

A variation of the attractor stimulus is analysed in Figure 5.

Figure 5: Unit square example: evolution of relative change in nominal density
and absolute nominal density with variations of the attractor stimulus for
cortical bone

The dotted line in Figure 5 based on a relatively high attractor stimulus
displays low density results. Bone growth can not be achieved until the third
load level. Reducing the attractor stimulus (dashed line) results in an
overall increased gain in nominal density.

Since cortical bone should neither grow nor absorb faster than cancellous
bone, in the following, the same attractor stimulus will be assumed for both
microstructure types.

###  3.3 Investigating the interaction between cancellous and cortical bone

To analyse the local distribution of cortical and cancellous bone and their
interaction, a two-dimensional proximal femur head is examined.

Three typical load cases representing the midstance phase of gait, extreme
abduction and adduction are applied to the model as reported in (Carter and
Beauprè 2010). In the calculations bi-linear element expansions with a 2x2
Gauss quadrature rule are used. The mesh, see Figure 6 (left), consists of
13317 elements and 13602 nodes. The parameters are set as follows: Poisson’s
ratio $\nu^{s}=0.2$, attractor stimulus $\Psi_{0}^{a}=0.01$, exponent $n=2$
and true solid material density $\rho_{0}^{s}=2.0$.

An initial nominal cortical density of $\rho_{0}^{\star}=1.825$, leading to an
initial volume fraction of
$\nu_{0}^{\star}={\rho_{0}^{\star}}/{\rho_{0}^{s}}={1.825}/{2}=0.9125$ is set
for cortical bone (elements at the external boundary), and
$\rho_{0}^{\star}=0.55$ with
$\nu_{0}^{\star}={\rho_{0}^{\star}}/{\rho_{0}^{s}}={0.55}/{2}=0.275$ for
cancellous bone. To avoid an abrupt unrealistic change between elements of
different microstructure types, a small transition zone between both regions
is created using a sigmoid function, see Figure 6 (middle).

Four different models are studied and compared in the following, all using the
homogenisation exponent $N=n=2$. Selected parameters and thus the difference
between the four models are illustrated in Table 3.3.

non-dimensionalised reference parameters of the four proxima femur models with
$N=n=2$ ${\nu_{0}^{\star}}$ of cancellous bone ${\nu_{0}^{\star}}$ of cortical
bone $m$ of cancellous bone $m$ of cortical bone Model 1 0.275 — 3 — Model 2
0.275 0.9125 3 3 Model 3 0.275 0.9125 3 10 Model 4 0.275 0.9125 3 18

The first model represents the full spongy bone as reported in previous
literature with an initial nominal density of $\rho_{0}^{\star}=0.55$. The
second model differentiates between cortical and cancellous bone using the
initial volume fraction or rather the initial nominal density distribution as
depicted in Figure 6 (middle).

Figure 6: Femur example: mesh and the distribution of the initial nominal
density and the nominal Young’s modulus

With a true elastic modulus of $E^{s}=10000$, which is similarly documented
for solid trabecular material by (Rho et al. 1993) and (Yamada et al. 2014),
and based on the quadratic relation $N=n=2$ between stiffness and nominal
density according to (Gibson 2005), this results in a nominal stiffness of
$E^{\star}=756.25$ in the spongy region at an initial nominal density of
$\rho_{0}^{\star}=0.55$, which roughly reflects the magnitude of stiffness as
reported in (Wirtz et al. 2000) and which is also utilised in the full spongy
model. The maximum nominal elastic modulus of cortical bone is approximately
$E^{\star}=8327$ and thus significantly higher, see Figure 6 (right).

Motivated by the results of section 3.2, the third and fourth model include a
variation of the exponent $m$. While the exponent is always set to $m=3$ in
the first and second model, the others include an increase within the cortical
bone area to a maximum of $m=10$ in the third and $m=18$ in the fourth model.
The distribution of the exponent $m$ follows the respective volume fraction
which indicates the type of bone microstructure. A slight transition is
implemented between the two bone types using an hyperbolic tangent function,
see Figure 7.

Figure 7: Femur example: mesh zoom of femur head and distribution of exponent
$m$ with $m=10$ and $m=18$ for cortical bone and $m=3$ for cancellous bone

The results for the evolution of the nominal density are displayed for each
model at different points in time in Figure 8.

Figure 8: Femur example with $N=n=2$: temporal evolution of absolute nominal
density of the full spongy model (top row) and with the differentiation of
cancellous and cortical bone (second row) and a variation in $m$ for cortical
bone: $m=10$ third row, $m=18$ last row

Comparing the first two rows, it becomes clear that the pattern of density
evolution has changed remarkably. The density growth in the femoral head and
neck is significantly less intense when taking cortical bone into
consideration. Density loss increases and slightly shifts in the
intertrochanter region, which leads to an obvious change in density pattern.
In areas with subliminal stress like the medial femur head, cortical bone
decomposes slightly and at the same time ensures a less distinct density
reduction in the adjacent spongy region.

Increasing the exponent $m$ in cortical bone (see row 3,4) suppresses its
resorption and cortical bone remains, although weaker and thinner. Density
loss of cancellous bone in the intertrochanter region is reduced and an
increase in density in the femoral neck, albeit small, can be observed
compared to the second model.

In summary, it appears that cortical bone acts as a kind of stiffer shell.
Therefore, the interior of the femur head bears less load and cancellous bone
does not adapt its density as much as in the first, full spongy model.

###  3.4 Investigating the impact of the homogenisation exponent N

Finally, the alternative volumetric mixture rule with $N=1$ shall be examined.
We first re-analyse the simple square example as described in section 3.1 with
a variation in the initial volume fraction and homogenisation exponent $N=1$,
whereby all the other parameters remain unchanged. The temporal evolution of
the nominal density is shown in Figure 9.

Figure 9: Unit square example: evolution of relative change in nominal density
and absolute nominal density with N=1 and a variation in the initial volume
fraction representing cortical and spongy bone

Comparing the results of $N=1$ with those of section 3.1 in Figure 3, it can
be seen that density growth is reduced for an initial volume fraction of
$\nu_{0}^{\star}=0.25$ (gray lines), whereas a reduction in density growth for
an initial volume fraction of $\nu_{0}^{\star}=0.9$, which even leads to a
small density loss in the first load step, is hardly visible (black lines).

The reason for this difference lies in the modified relationship between
stiffness and nominal density from quadratic to linear. Whereas the initial
nominal stiffness has quadrupled ($1/0.25$) at an initial volume fraction of
$\nu_{0}^{\star}=0.25$, the initial volume fraction of $\nu_{0}^{\star}=0.9$
yields only a factor of $1.111=1/0.9$.

Keeping this in mind, we again consider the proximal femur model, similar to
section 3.3. Setting the true elastic modulus to $E^{s}=10000$ yields an
initial nominal stiffness of $E^{\star}=2750$ for spongy bone with
$\nu_{0}^{\star}=0.275$ via the linear relationship, which is far too high.
Reducing the Young’s modulus of the solid material to $E^{s}=2750$ results in
the same $E^{\star}$ for cancellous bone as in section 3.3, see Figure 10.
However, the stiffness for cortical bone becomes much too small.

Figure 10: Femur example: distribution of nominal Young’s modulus with
exponent $N=1$ and $E^{s}=2750$ (left) versus $N=n=2$ and $E^{s}=10000$
(right)

Nevertheless, Figure 11 shows the results of the four different models
including the full spongy model, the distinction between cancellous and
cortical bone with $m=3$, and two variations of the exponent $m$ for cortical
bone, all calculated with a true elastic modulus of $E^{s}=2750$ and an
homogenisation exponent of $N=1$. Since the true Young’s modulus has been
properly adjusted and thus the nominal stiffness as well as all the other
parameters in the full spongy model correspond to those of section 3.3, the
result is the same as can be seen in the first rows of Figure 11 and 8.

When comparing the different homogenisation exponents for the second model, we
note the significantly higher density growth in the area of the femoral neck,
but also an increased bone loss in the area of the medial femoral head.
Furthermore, the pattern of density evolution more closely resembles that of
the full spongy model, especially in the area of the greater trochanter. This
effect is enhanced when increasing the exponent $m$ for cortical bone, see
Figure 11 row 3 and 4. In this scenario, cortical bone decomposes much more in
subliminally loaded regions with $N=1$, especially with a small exponent $m$.

Figure 11: Femur example with $N=1$: temporal evolution of absolute nominal
density without (top row) and with cortical bone and a variation in $m$ for
cortical bone: $m=3$ second row, $m=10$ third row, $m=18$ last row

Overall, the application of the volumetric mixture rule with a homogenisation
exponent of $N=1$ leads to a nominal density distribution that is more similar
to the full-spongy model. However, ensuring a stiffness-density ratio for
cancellous bone as motivated by literature will not simultaneously achieve a
realistic ratio in the cortical area. A further possible option would
therefore be a separate definition of an individual true Young’s modulus for
cortical bone, but will not be examined in this survey.

##  4 Discussion and conclusion

In this study the classical open-system model for bone remodelling by (Kuhl
and Steinmann 2003) is modified to account for the initial nominal density
distribution especially at the end of long bones, where two microstructure
types, the spongiosa and the compacta can be found. A distinction is made
between the two bone microstructure types by means of the initial volume
fraction under the precondition that the same true bone material is assumed
for both types. The model was applied to a one-dimensional benchmark problem
and the two-dimensional proximal femur head, which is a representative
biomechanical problem. We demonstrated that the additional modelling of
cortical bone influences the density evolution of cancellous bone and should
therefore be considered in future patient-specific simulations.

It should be emphasised that this is a first approach to take the local
initial density distribution, the remodelling of cancellous and cortical bone
and their interaction into account. Still, realistic material parameters have
to be defined via analyses of clinical data, as well as more exact modelling
of the boundary conditions and the influence of muscles tendons and ligaments
in everyday and extreme loading scenarios.

## Disclosure statement

None of the authors has a financial or personal relationship with other people
or organizations that could inappropriately bias his/her work.

## References

  * An and Draughn (2000) An YH, Draughn RA. 2000. Mechanical testing of bone and the bone-implant interface. CRC.
  * Barkaoui et al. (2019) Barkaoui A, Kahla RB, Merzouki T, Hambli R. 2019. Numerical simulation of apparent density evolution of trabecular bone under fatigue loading: effect of bone initial properties. Journal of Mechanics in Medicine and Biology. 19:1950041.
  * Carter and Beauprè (2010) Carter DR, Beauprè GS. 2010. Skeletal function and form. Cambridge University Press.
  * Cowin and Hegedus (1976) Cowin SC, Hegedus DH. 1976. Bone remodeling i: theory of adaptive elasticity. Journal of Elasticity. 6:313–326.
  * Frost (1987) Frost HM. 1987. Bone mass and the mechanostat: a proposal. The Anatomical record. 219:1–9.
  * Gibson (2005) Gibson LJ. 2005. Biomechanics of cellular solids. Journal of Biomechanics. 38:377–399.
  * Hambli (2014) Hambli R. 2014. Connecting mechanics and bone cell activities in the bone remodeling process: An integrated finite element modeling. Frontiers in Bioengineering and Biotechnology. 2.
  * Harrigan and Hamilton (1993) Harrigan TP, Hamilton JJ. 1993. Finite element simulation of adaptive bone remodelling: A stability criterion and a time stepping method. International Journal for Numerical Methods in Engineering. 36:837–854.
  * Harrigan and Hamilton (1994) Harrigan TP, Hamilton JJ. 1994. Necessary and sufficient conditions for global stability and uniqueness in finite element simulations of adaptive bone remodeling. International Journal of Solids and Structures. 31:97–107.
  * Holzapfel (2001) Holzapfel GA. 2001. Nonlinear solid mechanics - a continuum approach for engineering. West Sussex, England: John Wiley and Sons.
  * Huiskes et al. (1987) Huiskes R, Weinans H, Grootenboer HJ, Dalstra M, Fudala B, Slooff TJ. 1987. Adaptive bone-remodeling theory applied to prosthetic-design analysis. Journal of Biomechanics. 20:1135–1150.
  * Jacobs et al. (1995) Jacobs CE, Levenston ME, Beauprè GS, Simo JC, Carter DR. 1995. Numerical instabilities in bone remodeling simulations: The advantages of a node-based finite element approach. J Biomechanics. 28/4:449–459.
  * Kuhl and Steinmann (2003) Kuhl E, Steinmann P. 2003. Theory and numerics of geometrically non-linear open system mechanics. Int J Numer Meth Engng. 58(11):1593–1615.
  * Liedtke et al. (2017) Liedtke H, McBride A, Sivarasu S, Roche S. 2017. Computational simulation of bone remodelling post reverse total shoulder arthroplasty. arXiv preprint arXiv:170508324.
  * Nutu (2018) Nutu E. 2018. Role of initial density distribution in simulations of bone remodeling around dental implants. Acta of bioengineering and biomechanics. 20:23–31.
  * Papastavrou et al. (2020a) Papastavrou A, Schmidt I, Deng K, Steinmann P. 2020a. On age-dependent bone remodeling. Journal of Biomechanics. 103:109701.
  * Papastavrou et al. (2020b) Papastavrou A, Schmidt I, Steinmann P. 2020b. On biological availability dependent bone remodeling. Computer methods in biomechanics and biomedical engineering. 23:432–444.
  * Rho et al. (1993) Rho JY, Ashman RB, Turner CH. 1993. Young’s modulus of trabecular and cortical bone material: Ultrasonic and microtensile measurements. Journal of Biomechanics. 26:111–119.
  * Weinans et al. (1992) Weinans H, Huiskers R, Grootenboer HJ. 1992. The behavior of adaptive bone-remodeling simulation models. J Biomechanics. 25:1425–1441.
  * Wirtz et al. (2000) Wirtz DC, Schiffers N, Pandorf T, Radermacher K, DieternWeichert, Forst R. 2000\. Critical evaluation of known bone material properties to realize anisotropic fe-simulation of the proximal femur. Journal of Biomechanics. 33:1325–1330.
  * Wolff (1892) Wolff J. 1892. Das gesetz der knochentransformation. Berlin: Hirschwald.
  * Yamada et al. (2014) Yamada S, Tadano S, Fukuda S. 2014. Nanostructure and elastic modulus of single trabecula in bovine cancellous bone. Journal of Biomechanics. 47:3482–3487.

###### List of Figures

  1. 1 Homogenised macro scale with nominal density $\rho_{0}(t)$ versus heterogeneous micro scale with true density $\rho_{0}^{s}$
  2. 2 Unit square example: specimen and load function
  3. 3 Unit square example: evolution of relative change in nominal density and absolute nominal density with a variation in the initial volume fraction representing cortical and spongy bone
  4. 4 Unit square example: evolution of relative change in nominal density and absolute nominal density with variations of $m$ for cortical bone
  5. 5 Unit square example: evolution of relative change in nominal density and absolute nominal density with variations of the attractor stimulus for cortical bone
  6. 6 Femur example: mesh and the distribution of the initial nominal density and the nominal Young’s modulus
  7. 7 Femur example: mesh zoom of femur head and distribution of exponent $m$ with $m=10$ and $m=18$ for cortical bone and $m=3$ for cancellous bone
  8. 8 Femur example with $N=n=2$: temporal evolution of absolute nominal density of the full spongy model (top row) and with the differentiation of cancellous and cortical bone (second row) and a variation in $m$ for cortical bone: $m=10$ third row, $m=18$ last row
  9. 9 Unit square example: evolution of relative change in nominal density and absolute nominal density with N=1 and a variation in the initial volume fraction representing cortical and spongy bone
  10. 10 Femur example: distribution of nominal Young’s modulus with exponent $N=1$ and $E^{s}=2750$ (left) versus $N=n=2$ and $E^{s}=10000$ (right)
  11. 11 Femur example with $N=1$: temporal evolution of absolute nominal density without (top row) and with cortical bone and a variation in $m$ for cortical bone: $m=3$ second row, $m=10$ third row, $m=18$ last row