# Radiative first-order phase transitions to next-to-next-to-leading-order

Andreas Ekstedt<EMAIL_ADDRESS>Department of Physics and Astronomy,
Uppsala University, Box 516, SE-751 20 Uppsala, Sweden  II. Institute of
Theoretical Physics, Universität Hamburg, D-22761, Hamburg, Germany  Deutsches
Elektronen-Synchrotron DESY, Notkestr. 85, 22607 Hamburg, Germany  Oliver
Gould<EMAIL_ADDRESS>School of Physics and Astronomy,
University Park, University of Nottingham, Nottingham NG7 2RD, United Kingdom
Johan Löfgren<EMAIL_ADDRESS>Department of Physics and
Astronomy, Uppsala University, Box 516, SE-751 20 Uppsala, Sweden

(May 16, 2022)

###### Abstract

We develop new perturbative tools to accurately study radiatively-induced
first-order phase transitions. Previous perturbative methods have suffered
internal inconsistencies and been unsuccessful in reproducing lattice data,
which is often attributed to infrared divergences of massless modes (the Linde
problem). We employ a consistent power counting scheme to perform
calculations, and compare our results against lattice data. We conclude that
the consistent expansion removes many previous issues, and indicates that the
infamous Linde problem is not as big a factor in these calculations as
previously thought.

##  I Introduction

A first-order phase transition at the electroweak scale would have far-
reaching consequences for the early universe, for example by triggering
baryogenesis Kuzmin  _et al._ (1985), and thereby explain the matter-
antimatter asymmetry. In general such a violent transition would echo across
the universe; leaving a stochastic background of gravitational waves in its
wake. This makes searching for the telltale spectral shape Hindmarsh  _et al._
(2017) a key scientific objective of current and planned gravitational wave
observatories Arzoumanian  _et al._ (2020); Amaro-Seoane  _et al._ (2017);
Kawamura  _et al._ (2011); Ruan  _et al._ (2020); El-Neaj  _et al._ (2020).
Besides giving clear evidence for a first-order transition, a detected signal
would constrain models of dark matter Schwaller (2015); Croon  _et al._
(2018), baryogenesis Vaskonen (2017), inflation Addazi  _et al._ (2019) and
grand unification Croon  _et al._ (2019).

Yet current tools have a hard time making reliable predictions, particularly
so for phase transitions which are _radiatively induced_ ,111We here use the
term “radiative” to denote radiative corrections from energy scales much lower
than the temperature. where perturbative calculations often cannot distinguish
between first-order, second-order, and crossover transitions. Alternatively,
lattice Monte-Carlo simulations can provide quantitatively reliable
predictions, up to (small) statistical uncertainties. But they are slow,
typically requiring thousands of CPU hours to accurately determine the
thermodynamics of a single parameter point.

As a consequence, even with its problems, perturbation theory is the only tool
capable of scanning the high-dimensional parameter spaces of theories beyond
the Standard Model, which makes it imperative to put perturbative calculations
on a solid footing.

In this paper we vie to do just that. In particular, we construct a new
perturbative expansion for models with radiatively-induced phase transitions.
Furthermore, to ensure accuracy we compare the results with lattice data. In
the process we draw attention to prevailing problems with standard
perturbative approaches in this context (for an overview see Ref. Croon  _et
al._ (2021)).

##  II Theoretical challenges

Consistent calculations of physical observables should be renormalization-
scale invariant Gould and Tenkanen (2021), gauge invariant Buchmuller  _et
al._ (1994); Hirvonen  _et al._ (2022); Löfgren  _et al._ (2021), and free
from infrared (IR) divergences Linde (1980); Laine (1994, 1995).

In the present context, most significant renormalization-scale dependence
originates from a hierarchy between the temperature and particle masses. Such
a hierarchy is typically present at phase transitions in weakly coupled
theories because large temperatures are needed for thermal loop-corrections to
overpower the tree-level potential. And so, because different loop orders are
mixed, the expansion is reorganized. For example, two-loop calculations are
required to achieve any parametric cancellation of the scale dependence Gould
and Tenkanen (2021). Additionally, large logarithms arise due to the hierarchy
of scales. Both these issues can be remedied by working with a dimensionally-
reduced effective field theory (EFT) Farakos  _et al._ (1994); Kajantie  _et
al._ (1996a); Braaten and Nieto (1995, 1996). In particular, renormalization-
group equations of the EFT can be used to introduce a second renormalization
scale, which can be chosen independently, thereby eliminating large
logarithms.

This hierarchy of scales also leads to a gauge-dependence problem, which in
the case of tree-level barriers can be resolved by using a dimensionally
reduced EFT Croon  _et al._ (2021); Gould and Tenkanen (2021). Yet this gauge-
dependence problem is more sinister for radiatively induced transitions, where
standard methods yield IR divergences Laine (1994, 1995).

With these considerations in mind, we want to appraise how well perturbation
theory performs. To do so we study first-order phase transitions in the three-
dimensional (3d) $\mathrm{SU}(2)$ Higgs EFT. This model has been extensively
studied using lattice simulations Kajantie  _et al._ (1996b); Gurtler  _et
al._ (1997); Rummukainen  _et al._ (1998); Moore and Rummukainen (2001); Gould
_et al._ (2022), due to its relevance for the electroweak phase transition.
The availability of lattice data allows us to determine the quantitative
reliability of perturbation theory as a function of the expansion parameter.
We also show how the perturbative expansion converges by pushing calculations
to next-to-next-to-leading order (NNLO).

The 3d EFT that we study describes phase transitions of many extensions of the
electroweak sector. This is partly because, in the vicinity of the phase
transition, the degrees of freedom driving it become anomalously light, and
any additional degrees of freedom need only be heavy compared with those of
the EFT. In addition, the $\mathrm{U}(1)$ hypercharge of the electroweak
theory has a relatively small effect on the thermodynamics Kajantie  _et al._
(1997a). In the literature, the 3d $\mathrm{SU}(2)$ Higgs EFT has been used to
describe phase transitions in the Minimally Supersymmetric Standard Model
Laine (1996); Cline and Kainulainen (1996); Farrar and Losada (1997); Cline
and Kainulainen (1998), the two Higgs doublet model Losada (1997); Andersen
(1999); Andersen  _et al._ (2018); Gorda  _et al._ (2019), the singlet scalar
extension of the Standard Model Brauner  _et al._ (2017); Gould  _et al._
(2019), and the real triplet scalar extension of the Standard Model Niemi  _et
al._ (2019). The model also serves as a useful prototype for models with a
radiatively-induced barrier.

##  III A modified expansion

This paper focuses on the 3d EFT, which contains all the phase-dependent
physics of the transition, so the exact form of the four-dimensional (4d)
theory is not of interest. The action of the EFT, barring gauge-fixing terms,
is

$\displaystyle S_{3}=\int
d^{3}x\left[\frac{1}{4}F_{ij}^{a}F_{ij}^{a}+D_{i}\Phi^{\dagger}D_{i}\Phi+V(\Phi)\right],$
(1)

where $i,j$ run over the spatial indices and $a$ runs over the adjoint indices
of $\mathrm{SU}(2)$, $F_{ij}^{a}$ is the gauge field strength tensor with
associated gauge coupling $g_{3}$, $\Phi$ is the Higgs field, in the
fundamental representation of $\mathrm{SU}(2)$, and $D_{i}$ is the gauge-
covariant derivative. Note that the metric for the spatial indices is positive
definite, which implies $D_{i}=D^{i}$. The notation is standard and follows
Ref. Farakos  _et al._ (1994). In addition, we work with a class of $R_{\xi}$
gauges Fujikawa  _et al._ (1972); Fukuda and Kugo (1976), with gauge parameter
$\xi$, in the standard convention of Ref. Peskin and Schroeder (1995).

The tree-level potential is

$\displaystyle
V(\Phi)=m^{2}_{3}\Phi^{\dagger}\Phi+\lambda_{3}(\Phi^{\dagger}\Phi)^{2}.$ (2)

All three parameters, $g_{3}^{2}$, $m^{2}_{3}$ and $\lambda_{3}$ are
effective—temperature dependent—parameters. The couplings are related to their
4d counterparts, and the temperature $T$, as $g_{3}^{2}\approx g^{2}T$ and
$\lambda_{3}\approx\lambda T$. The parameter $m^{2}_{3}$ is positive at high
temperatures and negative at low temperatures—schematically
$m^{2}_{3}\approx-\alpha+\gamma T^{2}$ for positive constants $\alpha$,
$\gamma$.

The phase structure of the model can be studied perturbatively via the
effective potential. As such we make the usual expansion,
$\Phi\to(0,\tfrac{1}{\sqrt{2}}\phi)+\Phi$, where $\phi$ denotes the vacuum
expectation value (vev). In general $\phi$ can either be zero (symmetric-
phase) or non-zero (broken phase). These phases have different free energy
densities; which, when degenerate, identifies a first-order phase transition.

Looking at the tree-level potential in terms of $\phi$, it would appear that
no first-order transition, with a barrier between coexisting phases, is
possible. Yet a barrier can be generated by vector-bosons loops Arnold and
Espinosa (1993). Indeed, incorporating the diagram shown in figure 1(a) gives
an effective potential

$\displaystyle
V_{\text{eff}}(\phi)=\frac{1}{2}m_{3}^{2}\phi^{2}-\frac{g_{3}^{3}}{16\pi}|\phi|^{3}+\frac{\lambda_{3}}{4}\phi^{4}+\ldots$
(3)

Taking a closer look at equation (3), the symmetric minimum ($\phi=0$) has
lowest energy for large $m_{3}^{2}$; while the broken minimum ($\phi\neq 0$)
has lowest energy for small $m_{3}^{2}$. These two minima overlap for some
intermediate $m_{3}^{2}$ value. This occurs when terms in the potential are of
similar size Arnold and Espinosa (1993):

$\displaystyle\phi$ $\displaystyle\sim\frac{g_{3}^{3}}{\lambda_{3}},$
$\displaystyle m_{3}^{2}$ $\displaystyle\sim\frac{g_{3}^{6}}{\lambda_{3}}.$
(4)

Therefore, the vector-boson mass $m_{A}^{2}=\tfrac{1}{4}g_{3}^{2}\phi^{2}$ and
the scalar-boson mass are related by
$m_{3}^{2}/m_{A}^{2}\sim\lambda_{3}/g_{3}^{2}$.

(a) LO

c NNLO

b NLO

Figure 1: Feynman loop diagrams which contribute to the effective potential.
Wavy lines denote vectors, plain lines scalars, and dotted lines denote
ghosts. Parts (a), (b) and (c) show the diagrams contributing to the effective
potential at LO, NLO and NNLO respectively.

Then, if this ratio,

$\displaystyle x\equiv\frac{\lambda_{3}}{g_{3}^{2}},$ (5)

is parametrically small we can formally integrate out the vector bosons, which
is what generates the barrier in equation (3). Because of this, the
perturbative expansion is performed in powers of $x$.

Note that in counting powers of $x$, one can directly identify the
(potentially infinite) classes of Feynman diagrams which contribute at a given
order in $x$, and make appropriate resummations, and subtractions, to avoid
double counting Arnold and Espinosa (1993); Ekstedt and Löfgren (2020).
Alternatively, as we do here, one can use EFT-techniques to systematically
integrate out vector bosons. When expanded strictly in $x$, both approaches
give the same results.

Within this EFT, perturbation theory should work well for $x\ll 1$. However,
as $x$ decreases the vector mass in the broken phase grows as $m_{A}\propto
1/x$, so we cannot let $x$ become too small. If $x$ is as small as $x\sim
g^{2}/(4\pi)^{2}$, the vector mass is of order $m_{A}\sim\pi T$, implying that
the high-temperature expansion is invalid, and hence the original 3d EFT no
longer faithfully describes the infrared physics of the 4d theory.

We therefore assume that $x$ satisfies $g^{2}/(4\pi)^{2}\ll x\ll 1$. A natural
scaling relation satisfying these bounds is the geometric midpoint, $x\sim
g/(4\pi)$, equivalent to $\lambda\sim g^{3}/(4\pi)$ Arnold and Espinosa
(1993); Ekstedt and Löfgren (2020).

In addition to $x$, we use $y\equiv\frac{m_{3}^{2}}{g_{3}^{4}}$ Kajantie  _et
al._ (1996b). With these dimensionless variables, and scaling $\phi\to
g_{3}\phi$, the leading-order potential is

$\displaystyle\frac{V_{\text{LO}}(\phi)}{g_{3}^{6}}=\frac{1}{2}y\phi^{2}-\frac{1}{16\pi}|\phi|^{3}+\frac{x}{4}\phi^{4}.$
(6)

In the broken phase, this is of order $\sim x^{-3}$.

In principle $|y|^{-1/2}$ could also act as an expansion parameter, as it
appears in loop corrections through Feynman integrals. However, close to the
critical temperature $y$ scales as $y\sim x^{-1}$, so it is enough to count
powers of $x$. Hereafter we set $g_{3}^{2}=1$; if necessary, factors of
$g_{3}^{2}$ can be reinstated by dimensional analysis.

In the EFT approach one first integrates out the heavy vector fields, and
highly energetic scalars with $\vec{p}\sim m_{A}$, to construct an EFT for the
light scalar fields Gould and Hirvonen (2021); Löfgren  _et al._ (2021);
Hirvonen  _et al._ (2022); Andreassen  _et al._ (2017); Patel and Radovcic
(2017). At LO this gives equation (6). Subleading corrections to the EFT
action appear at NLO in the $x$ expansion. These corrections come with integer
powers of $x$. After vector bosons are integrated out, the resulting EFT only
contains scalars. These scalars give loop corrections which are suppressed by
powers of $\lambda_{3}/m_{3}\sim x^{3/2}$.

Thus, in the fully coupled gauge-Higgs theory, the perturbative expansion is a
dual expansion in powers of $x$ and $x^{3/2}$ (up to logarithms), or
equivalently an expansion in powers of $x^{1/2}$ starting at order $x$.

That said, the potential to NNLO reads

$\displaystyle
V_{\text{eff}}=V_{\text{LO}}+xV_{\text{NLO}}+x^{3/2}V_{\text{NNLO}}+\ldots$
(7)

Here factors of $x$ only signify the suppression of higher-order terms.
Importantly, the expansion is organized in powers of $x$—not by loops.

This $x$-expansion describes radiatively-induced first-order phase transitions
in an IR-finite and gauge invariant manner. The IR-finiteness of the
$x$-expansion should be contrasted with the IR divergences that start to
appear at two-loops in a strict $\hbar$-expansion Laine (1994, 1995). In
addition, physical quantities would obtain a spurious gauge dependence if
computed by numerically minimizing the effective potential Nielsen (1975);
Fukuda and Kugo (1976) (for an example see Ref. Croon  _et al._ (2021));
whereas the $x$-expansion is gauge invariant at every order.

Achieving order-by-order gauge invariance for observables requires the
$x$-expansion to be performed strictly, so that all quantities are expanded in
powers of $x$, including intermediate quantities such as the Higgs vev. Doing
so leads to asymptotic expansions for gauge-invariant observables, in terms of
the gauge-invariant expansion parameter $x$. Order-by-order gauge invariance
then follows from the uniqueness of the coefficients of an asymptotic
expansion Bender  _et al._ (1999). With this perspective, the Nielsen-Fukuda-
Kugo identities Nielsen (1975); Fukuda and Kugo (1976) demonstrate the order-
by-order gauge invariance of the free energy density, which is equal to the
gauge-dependent effective action evaluated on solutions that extremize it.

As mentioned, the next-to-leading order (NLO) potential comes from integrating
out vector-bosons at two loops, while contributions from the scalar fields
appear first at NNLO. These should be computed within the effective
description for the light scalar fields. Thus, at LO the squared masses of the
Higgs and Goldstone fields are, respectively,

$\displaystyle m_{H}^{2}(\phi)$
$\displaystyle=\partial^{2}_{\phi}V_{\text{LO}},$ $\displaystyle
m_{G}^{2}(\phi)$ $\displaystyle=\phi^{-1}\partial_{\phi}V_{\text{LO}}.$ (8)

Utilising the LO potential rather than the tree-level potential here resums
the Higgs and Goldstone self-energies to LO in $x$; see figure 1c.

As mentioned, all quantities, including the minima, should be expanded in
powers of $x$. For example, the minimization condition is

$\displaystyle\partial_{\phi}V_{\text{eff}}(\phi)\left.\right|_{\phi=\phi_{\text{min}}}=0,$
(9)

where

$\displaystyle\phi_{\text{min}}=\phi_{\text{LO}}+x\phi_{\text{NLO}}+\mathellipsis,$
(10)

and $\phi_{\text{LO}}$ solves
$\partial_{\phi}V_{\text{LO}}(\phi)\left.\right|_{\phi=\phi_{\text{LO}}}=0$.
Higher-order terms of $\phi_{\text{min}}$ are found by using equation (10) in
a Taylor expansion of equation (9).

The effective potential evaluated at a minimum represents the free energy
density of that phase. And the difference in free energy density between
phases can be expressed as

$\displaystyle\Delta
V\equiv\left[V_{\text{eff}}(\phi_{\text{min}})-V_{\text{eff}}(0)\right]$ (11)

We say that a phase-transition occurs for some critical mass, or value of $y$,
defined by $\Delta V(y=y_{\text{c}})=0$. From this one can determine the
critical temperature $T_{\text{c}}$ by solving $y(T_{\text{c}})=y_{\text{c}}$,
and using the known temperature dependence of $y$ for a given 4d model. This
critical mass should also be found order-by-order in $x$, to wit
$y_{\text{c}}=y_{\text{LO}}+xy_{\text{NLO}}+\ldots$. To leading order the
critical mass is the solution of $\Delta
V_{\text{LO}}=\left[V_{\text{LO}}(\phi_{\text{LO}})-V_{\text{LO}}(0)\right]_{y=y_{\text{LO}}}=0$.
And the next-to-leading order critical mass is

$\displaystyle y_{\text{NLO}}=-\left.\frac{\Delta
V_{\text{NLO}}}{\partial_{y}\Delta
V_{\text{LO}}}\right|_{\phi=\phi_{\text{LO}},y=y_{\text{LO}}}.$ (12)

Consider now an observable, $F(\phi,y)=F_{\text{LO}}+xF_{\text{NLO}}+\ldots$,
evaluated at the critical mass. The expansion is of the form

$\displaystyle F(\phi_{\text{min}},y_{\text{c}})=F_{\text{LO}}$
$\displaystyle+F_{\text{NLO}}+y_{\text{NLO}}\partial_{y}F_{\text{LO}}+\bar{\phi}_{\text{NLO}}\partial_{\phi}F_{\text{LO}}+\dots,$
(13)

where $\bar{\phi}_{\text{NLO}}$ is given by
$\bar{\phi}_{\text{NLO}}=\phi_{\text{NLO}}+y_{\text{NLO}}\partial_{y}\phi_{\text{LO}}$,
and all terms are evaluated at $\phi_{\text{LO}}$, $y_{\text{LO}}$.

Figure 2: The renormalization-scale invariant quantity
$y_{\text{c}}-\beta_{y}\log\mu_{3}$, where $y_{\text{c}}$ is the critical
mass, $\mu_{3}$ is the renormalisation scale, and $\beta_{y}$ is the beta
function for $y$. Theoretical uncertainties for lattice data are shown as
(barely visible) error bars, while those for the perturbative calculations are
shown as coloured bands. The latter are estimated either by varying the
renormalization scale over $\mu_{3}/(2m_{A})\in[1/\sqrt{10},\sqrt{10}]$, or by
varying the $R_{\xi}$ gauge parameter over $\xi\in[0,5]$. The hashed region
shown on the right reflects the presence of a third unphysical phase in this
approximation. The orange data point shows the location of the endpoint of the
line of first-order phase transitions, as determined on the lattice.

Note the simplicity of using strict-perturbation theory as compared to a
numerical approach: there is no need to numerically minimize (the real parts
of) complicated potentials since everything is expressed in terms of the LO
vev and critical mass.

Let us now turn to the latent heat, which determines the strength of a first-
order phase transition. It is given by the temperature derivative of $\Delta
V$. This can be written in terms of derivatives with respect to the 3d
effective parameters $x$, $y$ and $g_{3}^{3}$, using the chain rule of
differentiation Farakos  _et al._ (1995). Dependence on the derivative with
respect to $g_{3}^{3}$ cancels at the critical temperature Kajantie  _et al._
(1996b). This leads us to calculate the following scalar condensates,

$\displaystyle\Delta\braket{\Phi^{\dagger}\Phi}$
$\displaystyle\equiv\frac{\partial}{\partial y}\Delta V,$
$\displaystyle\Delta\braket{(\Phi^{\dagger}\Phi)^{2}}$
$\displaystyle\equiv\frac{\partial}{\partial x}\Delta V,$ (14)

which determine the contribution to the latent heat from the infrared EFT
scale.

##  IV Results

Results at LO follow from the LO effective potential given in equation (6):
the sum of tree-level and the one-loop term in figure 1(a). The minima can be
found analytically, which yields simple expressions for thermodynamic
observables.

NLO corrections arise from two-loop diagrams with loop momenta of order the
vector boson mass. These diagrams, shown in figure 1b, appear through
integrating out the heavy vectors for the EFT of the light scalars.
Contributions to these diagrams from smaller loop momenta are subdominant, so
scalar masses can be set to zero within loop integrals. The total NLO
contribution to the effective potential is

$\displaystyle
V_{\text{NLO}}=\frac{\phi^{2}}{(4\pi)^{2}}\left(-\frac{51}{32}\log\frac{|\phi|}{\mu_{3}}-\frac{63}{32}\log\frac{3}{2}+\frac{33}{64}\right),$
(15)

where $\mu_{3}$ is the 3d renormalization scale Farakos  _et al._ (1994).

The NNLO potential consists of one-loop diagrams within the light scalar EFT.
Expanded strictly in $x$, this is equivalent to the diagrams in figure 1c.
These diagrams give

$\displaystyle
V_{\text{NNLO}}=-\frac{1}{12\pi}\left[\left(m_{H}^{2}(\phi)\right)^{3/2}+3\left(m_{G}^{2}(\phi)\right)^{3/2}\right].$
(16)

The resummation of vector petals follows from equation (8).

With the effective potential in hand, we can calculate desired observables.
Following the previous section we find

$\displaystyle y_{\text{c}}$
$\displaystyle=\frac{1-\frac{51}{2}x\log\tilde{\mu_{3}}-2\sqrt{2}x^{3/2}}{2(8\pi)^{2}x},$
(17) $\displaystyle\Delta\braket{\Phi^{\dagger}\Phi}_{\text{c}}$
$\displaystyle=\frac{1+\frac{51}{2}x+13\sqrt{2}x^{3/2}}{2(8\pi x)^{2}},$ (18)
$\displaystyle\Delta\braket{\left(\Phi^{\dagger}\Phi\right)^{2}}_{\text{c}}$
$\displaystyle=\frac{1+51x+14\sqrt{2}x^{3/2}}{4(8\pi x)^{4}}.$ (19)

where $\tilde{\mu}_{3}\equiv
e^{-\tfrac{11}{34}+\tfrac{42}{34}\log\tfrac{3}{2}}(8\pi x\mu_{3})\approx
1.19(8\pi x\mu_{3})$. These expressions are accurate up to $O(x^{2})$ in the
numerators.

Some information on the convergence of the perturbative expansion can be
gleaned by comparing the magnitudes of successive terms. One finds that NLO
terms in equations (17)-(19) can dominate over LO for $x\gtrsim 0.03$,
choosing a renormalization scale $\mu_{3}\sim 1/(8\pi x)$. However, in each
case the NNLO term remains smaller until significantly larger $x$. This
suggests that the NLO term is anomalously large, and the perturbative
expansion may be reliable at somewhat larger values of $x$, as indeed we find
in our comparison to lattice data.

The explicit logarithm of the renormalization scale in equation (17) matches
the implicit running of $y$, to ensure that $y-y_{\text{c}}$ is
renormalization scale invariant at this order; for the beta functions see
Refs. Farakos  _et al._ (1994); Gould  _et al._ (2022). By contrast, the
absence of logarithms in the scalar condensates reflects the renormalization
group invariance of these physical quantities. While the effective potential
is $\mu_{3}$ dependent at NLO, this dependence cancels exactly in the
observable quantity $\Delta V$.

The $x$-expansion is gauge-invariant and renormalization scale invariant
order-by-order. However, residual renormalization scale dependence can be
useful to estimate theoretical uncertainties. For $y_{\text{c}}$, this can be
carried out by including the exact renormalization group running, rather than
truncating at the order of $x$ calculated. Doing so includes running from
terms suppressed by higher powers of $x$ than we have calculated, which
therefore do not cancel. By contrast, the scalar condensates are simply
renormalization group invariant, both exactly and order-by-order, so there is
no way to introduce renormalization group running. This is nevertheless a
desirable feature of the $x$-expansion, as it gives unambiguous predictions at
each order.

Figures 2 and 3 compare perturbative and lattice results, the latter taken
from Refs. Kajantie  _et al._ (1996b); Gurtler  _et al._ (1997); Rummukainen
_et al._ (1998); Laine and Rummukainen (1999); Gould  _et al._ (2022). The
$x$-expansion for $y_{\text{c}}$ agrees well with the lattice over the entire
range of $x$ for which there is a first-order phase transition, and the
renormalization scale dependence gives a good estimate of the disagreement. By
choosing $\mu_{3}\sim 2m_{A}$, one avoids large logarithms. There is a
significant improvement from LO to NLO, and while the NNLO correction is
small, it does shift towards the lattice data. Figure 3 shows similarly good
agreement for the scalar condensates, though in this case the renormalization
scale invariance of the $x$-expansion means there are no uncertainty bands.
For the quadratic condensate, the NNLO correction improves agreement with the
lattice only for $x\lesssim 0.05$, suggesting that above this higher-order
terms become important. The quartic condensate shows the largest discrepancies
between the $x$-expansion and the lattice data, which is consistent with the
larger expansion coefficients for this quantity.

For comparison, in figure 2 we also show predictions based on a different
perturbative approach: numerically minimising the real part of the one-loop
and two-loop effective potentials. This approach is commonly taken in the
literature but, for radiatively-induced transitions, it is not an expansion in
any small parameter, because loop-level contributions are of the same size as
tree-level terms. Further, directly minimising the real part of the effective
potential cancels neither gauge nor renormalization-scale dependence at each
order. The large uncertainty at one-loop order in this approach is due to the
presence of a third unphysical phase for $R_{\xi}$ gauge parameters
$\xi\gtrsim 3$. Predictions for all the condensates in this approach gain an
imaginary part for $y<0$.

Figure 3: The jumps in the scalar condensates as functions of $x$, computed in
the $x$-expansion and on the lattice. Both condensates are manifestly gauge
and renormalization-scale invariant.

##  V Conclusions

We find that strong first-order phase transitions can be described rather well
by perturbation theory. We have developed an EFT approach to transitions with
radiatively-generated barriers, and performed calculations up to NNLO. The
main results of this paper are presented in equations (17)–(19) and in Figs. 2
and 3. For the SU(2) Higgs theory, we find good numerical agreement between
our perturbative expansion, in powers of $x$, and existing lattice data.

The existence of the $x$-expansion was recognised in the late 90s Kajantie
_et al._ (1997b, 1998); Moore and Rummukainen (2001), yet in the two
intervening decades it has not been used once, to our knowledge.222This
excludes the contemporary Ref. Hirvonen  _et al._ (2022), by one of the
present authors (J.L.). In this paper, we revive this forgotten idea, and make
the first major advancement since its introduction. We elucidate the
underlying structure of the $x$-expansion, and develop the tools necessary to
extend the expansion to higher orders. This allows us to push calculations to
NNLO for the first time —the lowest order at which resummations are necessary
within 3d.

The $x$-expansion is theoretically well-behaved in the sense that it does not
suffer from IR-divergences or gauge-dependence, and renormalization scale-
dependence cancels order-by-order. In addition, contrary to previous
calculations Arnold (1994); Arnold and Espinosa (1993), our results indicate
good convergence. Though the results are consistent with lattice calculations,
and better behaved than those of other methods, the expansion can still be
improved. As we have seen, going to NLO is very important to get a reasonably
accurate prediction. The next improvement to NNLO is numerically smaller, but
improves agreement with the lattice results for $x\lesssim 0.05$.

Our results can be contrasted with QCD, where there is a sizeable mismatch
between perturbative and lattice calculations of the free energy density Zhai
and Kastening (1995); Braaten and Nieto (1996). This can be partially ascribed
to the absence of a Higgs phase in QCD, within which perturbation theory is
comparatively well under control, but also to the relatively large magnitude
of the QCD gauge coupling, and its slow (logarithmic) approach to asymptotic
freedom.

The $x$-expansion, like all other perturbative approaches to non-Abelian gauge
fields at high temperature, breaks down at finite order due to the Linde
Problem Linde (1980). For the free energy density, nonperturbative infrared
physics contributes first at $O(g^{6}T^{4})$. For the quantities studied here,
this means that starting at N5LO, which is suppressed by $O(x^{3})$, the
$x$-expansion is fundamentally nonperturbative. However, the Linde Problem
cannot modify the leading five coefficients of the asymptotic expansion, and
hence it does not stop the $x$-expansion performing well for small enough $x$.

Future studies can utilize the approach taken here for other models. The
remaining higher-order terms in the expansion should also be calculated: N3LO
and N4LO, respectively $O(x^{2})$ and $O(x^{5/2})$. Doing so requires
calculation of three-loop diagrams. These two terms are the highest order
terms that are completely calculable in perturbation theory and therefore give
the final word on how well the expansion actually works.

Finally, given the importance of a strict expansion for equilibrium
quantities, we expect similar behaviour for near-equilibrium ones, like the
bubble-nucleation rate Ekstedt (2022). As such our results are an important
stepping-stone for accurate predictions of gravitational waves.

###### Acknowledgements.

We would like to thank Sinan Güyer and Kari Rummukainen for making available
their lattice results prior to publication of Ref. Gould  _et al._ (2022). The
work of A.E. has been supported by the Deutsche Forschungsgemeinschaft under
Germany’s Excellence Strategy - EXC $2121$ Quantum Universe - $390833306$; and
by the Swedish Research Council, project number VR:$2021$-$00363$. O.G. (ORCID
0000-0002-7815-3379) was supported by U.K. Science and Technology Facilities
Council (STFC) Consolidated grant ST/T000732/1.

## References

  * Kuzmin  _et al._ (1985) V. A. Kuzmin, V. A. Rubakov, and M. E. Shaposhnikov, Phys. Lett. B 155, 36 (1985).
  * Hindmarsh  _et al._ (2017) M. Hindmarsh, S. J. Huber, K. Rummukainen, and D. J. Weir, Phys. Rev. D 96, 103520 (2017), [Erratum: Phys.Rev.D 101, 089902 (2020)], arXiv:1704.05871 [astro-ph.CO] .
  * Arzoumanian  _et al._ (2020) Z. Arzoumanian _et al._ (NANOGrav), Astrophys. J. Lett. 905, L34 (2020), arXiv:2009.04496 [astro-ph.HE] .
  * Amaro-Seoane  _et al._ (2017) P. Amaro-Seoane _et al._ (LISA), (2017), arXiv:1702.00786 [astro-ph.IM] .
  * Kawamura  _et al._ (2011) S. Kawamura _et al._ , Class. Quant. Grav. 28, 094011 (2011).
  * Ruan  _et al._ (2020) W.-H. Ruan, Z.-K. Guo, R.-G. Cai, and Y.-Z. Zhang, Int. J. Mod. Phys. A 35, 2050075 (2020), arXiv:1807.09495 [gr-qc] .
  * El-Neaj  _et al._ (2020) Y. A. El-Neaj _et al._ (AEDGE), EPJ Quant. Technol. 7, 6 (2020), arXiv:1908.00802 [gr-qc] .
  * Schwaller (2015) P. Schwaller, Phys. Rev. Lett. 115, 181101 (2015), arXiv:1504.07263 [hep-ph] .
  * Croon  _et al._ (2018) D. Croon, V. Sanz, and G. White, JHEP 08, 203 (2018), arXiv:1806.02332 [hep-ph] .
  * Vaskonen (2017) V. Vaskonen, Phys. Rev. D 95, 123515 (2017), arXiv:1611.02073 [hep-ph] .
  * Addazi  _et al._ (2019) A. Addazi, A. Marcianò, and R. Pasechnik, MDPI Physics 1, 92 (2019), arXiv:1811.09074 [hep-ph] .
  * Croon  _et al._ (2019) D. Croon, T. E. Gonzalo, and G. White, JHEP 02, 083 (2019), arXiv:1812.02747 [hep-ph] .
  * Croon  _et al._ (2021) D. Croon, O. Gould, P. Schicho, T. V. I. Tenkanen, and G. White, JHEP 04, 055 (2021), arXiv:2009.10080 [hep-ph] .
  * Gould and Tenkanen (2021) O. Gould and T. V. I. Tenkanen, JHEP 06, 069 (2021), arXiv:2104.04399 [hep-ph] .
  * Buchmuller  _et al._ (1994) W. Buchmuller, Z. Fodor, and A. Hebecker, Phys. Lett. B 331, 131 (1994), arXiv:hep-ph/9403391 .
  * Hirvonen  _et al._ (2022) J. Hirvonen, J. Löfgren, M. J. Ramsey-Musolf, P. Schicho, and T. V. I. Tenkanen, JHEP 07, 135 (2022), arXiv:2112.08912 [hep-ph] .
  * Löfgren  _et al._ (2021) J. Löfgren, M. J. Ramsey-Musolf, P. Schicho, and T. V. I. Tenkanen, (2021), arXiv:2112.05472 [hep-ph] .
  * Linde (1980) A. D. Linde, Phys. Lett. B 96, 289 (1980).
  * Laine (1994) M. Laine, Phys. Lett. B 335, 173 (1994), arXiv:hep-ph/9406268 .
  * Laine (1995) M. Laine, Phys. Rev. D 51, 4525 (1995), arXiv:hep-ph/9411252 .
  * Farakos  _et al._ (1994) K. Farakos, K. Kajantie, K. Rummukainen, and M. E. Shaposhnikov, Nucl. Phys. B 425, 67 (1994), arXiv:hep-ph/9404201 .
  * Kajantie  _et al._ (1996a) K. Kajantie, M. Laine, K. Rummukainen, and M. E. Shaposhnikov, Nucl. Phys. B 458, 90 (1996a), arXiv:hep-ph/9508379 .
  * Braaten and Nieto (1995) E. Braaten and A. Nieto, Phys. Rev. D 51, 6990 (1995), arXiv:hep-ph/9501375 .
  * Braaten and Nieto (1996) E. Braaten and A. Nieto, Phys. Rev. D 53, 3421 (1996), arXiv:hep-ph/9510408 .
  * Kajantie  _et al._ (1996b) K. Kajantie, M. Laine, K. Rummukainen, and M. E. Shaposhnikov, Nucl. Phys. B 466, 189 (1996b), arXiv:hep-lat/9510020 .
  * Gurtler  _et al._ (1997) M. Gurtler, E.-M. Ilgenfritz, and A. Schiller, Phys. Rev. D 56, 3888 (1997), arXiv:hep-lat/9704013 .
  * Rummukainen  _et al._ (1998) K. Rummukainen, M. Tsypin, K. Kajantie, M. Laine, and M. E. Shaposhnikov, Nucl. Phys. B 532, 283 (1998), arXiv:hep-lat/9805013 .
  * Moore and Rummukainen (2001) G. D. Moore and K. Rummukainen, Phys. Rev. D 63, 045002 (2001), arXiv:hep-ph/0009132 .
  * Gould  _et al._ (2022) O. Gould, S. Güyer, and K. Rummukainen, (2022), arXiv:2205.07238 [hep-lat] .
  * Kajantie  _et al._ (1997a) K. Kajantie, M. Laine, K. Rummukainen, and M. E. Shaposhnikov, Nucl. Phys. B 493, 413 (1997a), arXiv:hep-lat/9612006 .
  * Laine (1996) M. Laine, Nucl. Phys. B 481, 43 (1996), [Erratum: Nucl.Phys.B 548, 637–638 (1999)], arXiv:hep-ph/9605283 .
  * Cline and Kainulainen (1996) J. M. Cline and K. Kainulainen, Nucl. Phys. B 482, 73 (1996), arXiv:hep-ph/9605235 .
  * Farrar and Losada (1997) G. R. Farrar and M. Losada, Phys. Lett. B 406, 60 (1997), arXiv:hep-ph/9612346 .
  * Cline and Kainulainen (1998) J. M. Cline and K. Kainulainen, Nucl. Phys. B 510, 88 (1998), arXiv:hep-ph/9705201 .
  * Losada (1997) M. Losada, Phys. Rev. D 56, 2893 (1997), arXiv:hep-ph/9605266 .
  * Andersen (1999) J. O. Andersen, Eur. Phys. J. C 11, 563 (1999), arXiv:hep-ph/9804280 .
  * Andersen  _et al._ (2018) J. O. Andersen, T. Gorda, A. Helset, L. Niemi, T. V. I. Tenkanen, A. Tranberg, A. Vuorinen, and D. J. Weir, Phys. Rev. Lett. 121, 191802 (2018), arXiv:1711.09849 [hep-ph] .
  * Gorda  _et al._ (2019) T. Gorda, A. Helset, L. Niemi, T. V. I. Tenkanen, and D. J. Weir, JHEP 02, 081 (2019), arXiv:1802.05056 [hep-ph] .
  * Brauner  _et al._ (2017) T. Brauner, T. V. I. Tenkanen, A. Tranberg, A. Vuorinen, and D. J. Weir, JHEP 03, 007 (2017), arXiv:1609.06230 [hep-ph] .
  * Gould  _et al._ (2019) O. Gould, J. Kozaczuk, L. Niemi, M. J. Ramsey-Musolf, T. V. I. Tenkanen, and D. J. Weir, Phys. Rev. D 100, 115024 (2019), arXiv:1903.11604 [hep-ph] .
  * Niemi  _et al._ (2019) L. Niemi, H. H. Patel, M. J. Ramsey-Musolf, T. V. I. Tenkanen, and D. J. Weir, Phys. Rev. D 100, 035002 (2019), arXiv:1802.10500 [hep-ph] .
  * Fujikawa  _et al._ (1972) K. Fujikawa, B. W. Lee, and A. I. Sanda, Phys. Rev. D 6, 2923 (1972).
  * Fukuda and Kugo (1976) R. Fukuda and T. Kugo, Phys. Rev. D 13, 3469 (1976).
  * Peskin and Schroeder (1995) M. E. Peskin and D. V. Schroeder,  _An Introduction to quantum field theory_ (Addison-Wesley, Reading, USA, 1995).
  * Arnold and Espinosa (1993) P. B. Arnold and O. Espinosa, Phys. Rev. D 47, 3546 (1993), [Erratum: Phys.Rev.D 50, 6662 (1994)], arXiv:hep-ph/9212235 .
  * Ekstedt and Löfgren (2020) A. Ekstedt and J. Löfgren, JHEP 12, 136 (2020), arXiv:2006.12614 [hep-ph] .
  * Gould and Hirvonen (2021) O. Gould and J. Hirvonen, Phys. Rev. D 104, 096015 (2021), arXiv:2108.04377 [hep-ph] .
  * Andreassen  _et al._ (2017) A. Andreassen, D. Farhi, W. Frost, and M. D. Schwartz, Phys. Rev. D 95, 085011 (2017), arXiv:1604.06090 [hep-th] .
  * Patel and Radovcic (2017) H. H. Patel and B. Radovcic, Phys. Lett. B 773, 527 (2017), arXiv:1704.00775 [hep-ph] .
  * Nielsen (1975) N. K. Nielsen, Nucl. Phys. B 101, 173 (1975).
  * Bender  _et al._ (1999) C. M. Bender, S. Orszag, and S. A. Orszag,  _Advanced mathematical methods for scientists and engineers I: Asymptotic methods and perturbation theory_ , Vol. 1 (Springer Science & Business Media, 1999).
  * Farakos  _et al._ (1995) K. Farakos, K. Kajantie, K. Rummukainen, and M. E. Shaposhnikov, Nucl. Phys. B 442, 317 (1995), arXiv:hep-lat/9412091 .
  * Laine and Rummukainen (1999) M. Laine and K. Rummukainen, Nucl. Phys. B Proc. Suppl. 73, 180 (1999), arXiv:hep-lat/9809045 .
  * Kajantie  _et al._ (1997b) K. Kajantie, M. Laine, K. Rummukainen, and M. E. Shaposhnikov, Nucl. Phys. B 503, 357 (1997b), arXiv:hep-ph/9704416 .
  * Kajantie  _et al._ (1998) K. Kajantie, M. Karjalainen, M. Laine, and J. Peisa, Nucl. Phys. B 520, 345 (1998), arXiv:hep-lat/9711048 .
  * Arnold (1994) P. B. Arnold, in  _8th International Seminar on High-energy Physics_ (1994) arXiv:hep-ph/9410294 .
  * Zhai and Kastening (1995) C.-x. Zhai and B. M. Kastening, Phys. Rev. D 52, 7232 (1995), arXiv:hep-ph/9507380 .
  * Ekstedt (2022) A. Ekstedt, (2022), arXiv:2205.05145 [hep-ph] .


# On Direction Preserving Discretizations for Computing Phase-Space Densities

David J. Chappell(a), Martin Richter(a,b) and Gregor Tanner(b)
(a) School of Science and Technology,
Nottingham Trent University,
Clifton Campus, Clifton Lane,
Nottingham NG11 8NS, UK
(b) School of Mathematical Sciences,
University of Nottingham,
University Park,
Nottingham NG7 2RD, UK

###### Abstract

Ray flow methods provide efficient tools for modelling wave energy transport
in complex systems at high-frequencies. We compare two Petrov-Galerkin
discretizations of a phase-space boundary integral model for stationary wave
energy densities in two-dimensional domains. The directional dependence is
approximated using a finite set of directions oriented into the domain from
the boundary. The propagation direction can be preserved across multi-
component domains when the directions within the local set for a given region
of the boundary are taken as a subset of a global direction set. In this work
we compare the use of piecewise constant and piecewise linear test functions,
which physically corresponds to the interpolation scheme used when the
transport is in a direction not belonging to the finite global set.

##  1 INTRODUCTION

Dynamical energy analysis (DEA) is an approach for modelling wave energy
densities at high-frequencies that was first proposed just over ten years ago
[1, 2, 3]. DEA is based on a linear integral operator description of phase-
space density transport along ray trajectories between positions on the
boundary of a domain or sub-domain. Recent developments have seen the
capability of DEA extended to industrial applications [1, 2], as well as
stochastic propagation through uncertain structures [4, 5, 6].

A Petrov-Galerkin method that is efficient for modelling densities with strong
directional dependence was recently proposed [7]. This class of problems has
usually proved problematic for the DEA method [1]. In [7], a basis
approximation using Dirac delta distributions was proposed to approximate the
directional dependence. The direction of transport can be preserved throughout
multi-domains provided that the Dirac delta specified direction set local to
any part of the boundary is inherited from a common global set of directions;
the propagation is then able to continue along rays with the same global
direction through different sub-domains. A beneficial consequence is that the
proposed methodology can be applied directly on complicated domains formed
from a potentially large number of simpler sub-domains as is typically the
case with the finite element type meshes used for industrial applications.

In this study we investigate the effect of modifying the choice of test
functions from the piecewise constant indicator functions employed in the
Petrov-Galerkin scheme proposed in [7]. For classical potential problems, it
is known that the convergence of Petrov-Galerkin schemes based on Dirac delta
distribution basis approximations together with splines as test functions
crucially depends on this choice [8, 9, 10]. In particular, higher order test
functions were shown to lower the regularity requirements on the boundary
data. In comparison to more standard projection methods for integral
equations, such as the collocation and Galerkin methods, these Petrov-Galerkin
schemes combine the efficiency and ease of implementation of collocation based
schemes with potentially even lower regularity requirements for convergence
than the Galerkin method [9]. These features are important in this work since
ray tracing solutions can often have low regularity and the added
dimensionality of working in phase-space means that a faster implementation
with fewer integrals to compute is desirable.

##  2 PHASE-SPACE BOUNDARY INTEGRAL MODEL

We are concerned with transporting densities along ray trajectories through
multi-domains $\Omega=\cup_{j=1}^{K}\Omega_{j}$. The ray dynamics in
$\Omega_{j}$ is defined by the Hamiltonian
$H_{j}(\mathbf{r},\mathbf{p})=|\mathbf{p}|/\eta(\mathbf{r})\equiv 1$ for
$j=1,\ldots,K$. The phase-space coordinates $(\mathbf{r},\mathbf{p})$ denote
the position $\mathbf{r}$ and momentum $\mathbf{p}$ vectors, respectively.
Each $\Omega_{j}$, $j=1,\ldots,K$ is assumed to be a convex polygon containing
a homogeneous medium. As a consequence we may write
$\eta(\mathbf{r})=\eta_{j}$ when $\mathbf{r}\in\Omega_{j}$ for $j=1,\ldots,K$,
where $\eta_{j}$, $j=1,\ldots,K$ are constants. When $\eta_{j}=c_{j}^{-1}$ is
taken to be the inverse of the phase speed in $\Omega_{j}$, then $H_{j}$
defines the ray trajectories obtained via leading order high-frequency
asymptotics for the Helmholtz equation

$c_{j}^{2}\>\Delta u+\omega^{2}u=0,$ (1)

at angular frequency $\omega$. The choice of $\eta$ can easily be modified for
flexural waves with nonlinear dispersion - see [7].

Figure 1: The boundary mapping
$\varphi_{i,j}(s_{j},p_{j})=(s^{\prime}_{i},p^{\prime}_{i})\in\Gamma_{i}\times(-\eta_{i},\eta_{i})$
of the phase-space boundary coordinates
$X_{j}=(s_{j},p_{j})\in\Gamma_{j}\times(-\eta_{j},\eta_{j})$, which maps to
the next intersection with a boundary edge followed by either a specular
reflection or transmission.

In order to express our model in boundary integral form it is convenient to
introduce the phase-space boundary coordinates $Y_{j}=(s_{j},p_{j})$ on
$\Gamma_{j}=\partial\Omega_{j}$, $j=1,\ldots,K$, where $s_{j}$ is an arclength
parameter for $\Gamma_{j}$ and $p_{j}=\eta_{j}\sin(\theta_{j})$ is the
tangential component of $\mathbf{p}$ at $s_{j}$. Here $\theta_{j}$ is the
angle formed between the ray (oriented into $\Omega_{j}$) and the inward
normal vector at $s_{j}$ as shown in Fig. 1. We next introduce a local
boundary flow map
$\varphi_{i,j}(s_{j},p_{j})=(s^{\prime}_{i}(s_{j},p_{j}),p^{\prime}_{i}(s_{j},p_{j}))$,
which describes the discrete evolution of the rays at times coinciding with
boundary intersections. As written here, $\varphi_{i,j}$ maps the boundary
phase-space coordinate $(s_{j},p_{j})$ on the boundary of $\Omega_{j}$ to
$(s^{\prime}_{i}(s_{j},p_{j}),p^{\prime}_{i}(s_{j},p_{j}))$ on the boundary of
$\Omega_{i}$. To write $\varphi_{i,j}$ in this form we have implicitly assumed
that either $i=j$, or $\Gamma_{i}$ and $\Gamma_{j}$ share a common edge
through which the ray can transmit. As before,
$p^{\prime}_{i}(s_{j},p_{j})=\eta_{i}\sin(\theta^{\prime}_{i}(s_{j},p_{j}))$
denotes the tangential slowness at $s_{i}^{\prime}$, and $\theta^{\prime}_{i}$
is the angle formed between the outgoing ray and the inward normal vector to
$\Gamma_{i}$ at $s^{\prime}_{i}$. In the simplest cases, $\theta^{\prime}_{i}$
is obtained from either a specular reflection when $i=j$, or if $i\neq j$ and
$\eta_{i}=\eta_{j}$ then the ray will continue in the same direction into
$\Omega_{i}$.

Phase-space densities are transported throughout $\Omega$ using a local
boundary operator $\mathcal{B}_{j}$, which transports a density $f$ along the
boundary flow $\varphi_{i,j}$ as follows [1, 7]

$\mathcal{B}_{j}[f](X_{i}):=\int
e^{-\mu_{j}D(X_{i},Y_{j})}w_{i,j}(Y_{j})\delta(X_{i}-\varphi_{i,j}(Y_{j}))f(Y_{j})\,\mathrm{d}Y_{j}.$
(2)

Here $X_{i}\in\Gamma_{i}\times(-\eta_{i},\eta_{i})$ for some $i=1,2,\ldots,K$
and $w_{i,j}$ has been introduced to incorporate reflection/transmission
coefficients. An exponential damping term with coefficient $\mu_{j}>0$ has
also been introduced, where $D(X_{i},Y_{j})$ denotes the distance between
$s_{j}$ and the solution point. The global boundary operator $\mathcal{B}$ is
then given by $\mathcal{B}=\sum_{j}\mathcal{B}_{j}$, where the sum is taken
over each $\Omega_{j}$ that shares an edge with $\Omega_{i}$, including
$\Omega_{i}$ itself.

The stationary boundary density $\rho$ can be expressed in terms of a Neumann
series

$\rho=\sum_{n=0}^{\infty}\mathcal{B}^{n}[\rho_{0}]=(I-\mathcal{B})^{-1}[\rho_{0}],$
(3)

where $\rho_{0}$ is a given initial boundary density and $\mathcal{B}^{n}$
represents $n$ iterates of the operator $\mathcal{B}$. Once $\rho$ has been
evaluated using (3), the interior density $\rho_{\Omega}$ can be calculated by
projecting onto a prescribed solution point $\mathbf{r}\in\Omega_{j}$ using
[2]

$\rho_{\Omega}(\mathbf{r})=\eta_{j}^{2}\int_{0}^{2\pi}e^{-\mu_{j}D(\mathbf{r},s_{j})}\rho(s_{j}(\mathbf{r},\Theta),p_{j}(\mathbf{r},\Theta))\,\mathrm{d}\Theta.$
(4)

Here, $\Theta\in[0,2\pi)$ is the polar angle parametrising trajectories
approaching $\mathbf{r}$ from $s_{j}(\mathbf{r},\Theta)\in\Gamma_{j}$ and $D$
is used to represent the length of the trajectory between
$\mathbf{r}\in\Omega_{j}$ and $s_{j}\in\Gamma_{j}$.

##  3 PETROV-GALERKIN DISCRETIZATION

Figure 2: The local direction coordinates
$\phi_{n}(s_{j})\in(-\frac{\pi}{2},\frac{\pi}{2})$, $n=1,2,3$ for $L=8$ global
directions given by $\Phi_{l}=2\pi(l-1)/8$ for $l=1,2,\ldots,8$. The local and
global directions are connected via $\phi_{n}=\gamma-\Phi_{n+1}$ for
$n=1,2,3$, where $\gamma\in[0,2\pi)$ is the global direction of the inward
normal vector at $s_{j}$ (here $\gamma=\pi/2$).

In this section we introduce two direction preserving discretizations of the
boundary operator (2) using Petrov-Galerkin projections in order to
numerically solve for $\rho$ by solving a discretised form of equation (3). We
first split $\Gamma_{j}$ into boundary elements $E^{j}_{m}$ for
$m=1,2,\ldots,M_{j}$ and define the global ray directions
$\Phi_{l}\in[0,2\pi)$ for $l=1,2,\ldots L$. These global direction are defined
anti-clockwise relative to the positive $x_{1}$-axis as specified in [7]. Here
we take $\Phi_{l}=2\pi(l-1)/L$, but note that this choice is flexible and can
easily be amended to include dominant transmission paths, where known. Let
$\phi_{n}(s_{j})\in(-\pi/2,\pi/2)$, $n=1,2,\ldots,N_{m}$ denote the local ray
directions at $s_{j}\in\Gamma_{j}$. The local directions at $s_{j}$ are simply
a subset of the global directions taken as all directions which satisfy the
property of being directed into $\Omega_{j}$ from $s_{j}$. The local
directions are also given a local numbering based on their direction relative
to the interior normal at $s_{j}$ \- see Fig. 2.

We apply an approximation of the form

$\rho(s_{j},p_{j})\approx\sum_{m=1}^{M_{j}}\sum_{n=1}^{N_{m}}\rho_{(j,m,n)}b_{m}(s_{j})\delta(p_{j}-\tilde{p}_{n}(s_{j})),\qquad
j=1,\dots,K,$ (5)

where $\tilde{p}_{n}(s_{j})=\eta_{j}\sin(\phi_{n}(s_{j}))$ and
$b_{m}(s_{j})=|E^{j}_{m}|^{-1/2}:=(\mathrm{diam}({E^{j}}_{m}))^{-1/2}$ for
$s_{j}\in{E^{j}}_{m}$ or $b_{m}(s_{j})=0$ otherwise. We impose a standard
(Bubnov) Galerkin projection in $s_{j}$ with the orthonormal basis $b_{m}$,
$m=1,2,\ldots,M_{j}$. It is therefore only in the momentum variable $p_{j}$
that we apply a Petrov-Galerkin projection and compare two possible choices of
test functions that are orthonormal in the $L^{2}$ inner product to
$\delta(p_{j}-\tilde{p}_{n}(s_{j}))$ for $n=1,2,\ldots,N_{m}$. The first set
of test functions that we consider are the indicator functions

$\chi_{n}(p_{j})=\tilde{\chi}_{n}(\arcsin(p_{j}/\eta_{j}))$

introduced in [7], where $\tilde{\chi}_{n}(\theta_{j})=1$ if
$\theta_{j}\in((\phi_{n-1}+\phi_{n})/2,(\phi_{n}+\phi_{n+1})/2)$ and zero
otherwise. The second set of test functions that we propose are the piecewise
linear hat functions $B_{n}(p_{j})=\tilde{B}_{n}(\arcsin(p_{j}/\eta_{j}))$,
where

$\tilde{B}_{n}(\theta_{j})=\left\\{\begin{array}[]{ccc}\displaystyle
1+\frac{\theta_{j}-\phi_{n}}{\phi_{n}-\phi_{n-1}}&\mathrm{if}&\theta_{j}\in(\phi_{n-1},\phi_{n}],\vspace{2mm}\\\
\displaystyle
1+\frac{\phi_{n}-\theta_{j}}{\phi_{n+1}-\phi_{n}}&\mathrm{if}&\theta_{j}\in(\phi_{n},\phi_{n+1}],\vspace{2mm}\\\
0&&\mathrm{otherwise.}\end{array}\right.$

In both cases we note that

$\langle\delta(\cdot-\tilde{p}_{n}(s_{j})),\chi_{n^{\prime}}\rangle_{L^{2}(-\eta_{j},\eta_{j})}=\langle\delta(\cdot-\tilde{p}_{n}(s_{j})),B_{n^{\prime}}\rangle_{L^{2}(-\eta_{j},\eta_{j})}=0$
(6)

for all $n\neq n^{\prime}$, and if $n=n^{\prime}$ then the right hand side of
(6) is instead one.

We now apply the various Galerkin projections described above to the operator
$\mathcal{B}$. Using piecewise linear test functions in direction leads to a
matrix representation $B$ of $\mathcal{B}$ as follows:

$\displaystyle B_{I,J}$ $\displaystyle=$
$\displaystyle\int_{\Gamma_{j}\times(-\eta_{j},\eta_{j})}e^{-\mu_{j}D(\varphi_{i,j}(Y_{j}),Y_{j})}w_{i,j}(Y_{j})b_{m}(s_{j})b_{m^{\prime}}(s^{\prime}_{i}(Y_{j}))\delta(p_{j}-\tilde{p}_{n}(s_{j}))B_{n^{\prime}}(p^{\prime}_{i}(Y_{j}))\,\mathrm{d}Y_{j}\vspace{2mm}$
$\displaystyle=$
$\displaystyle\int_{\Gamma_{j}}e^{-\mu_{j}D_{i}(s_{j})}w_{i,j}(s_{j},\tilde{p}_{n}(s_{j}))b_{m}(s_{j})b_{m^{\prime}}(s^{\prime}_{i}(s_{j},\tilde{p}_{n}(s_{j})))B_{n^{\prime}}(p^{\prime}_{i}(s_{j},\tilde{p}_{n}(s_{j})))\,\mathrm{d}s_{j}$
$\displaystyle=$
$\displaystyle\frac{w_{i,j}(\tilde{p}_{n})}{|E^{j}_{m}|^{1/2}}\int_{E^{j}_{m}}e^{-\mu_{j}D_{i}(s_{j})}b_{m^{\prime}}(s^{\prime}_{i}(s_{j},\tilde{p}_{n}(s_{j})))B_{n^{\prime}}(p^{\prime}_{i}(s_{j},\tilde{p}_{n}(s_{j})))\,\mathrm{d}s_{j},$

where $I=(i,m^{\prime},n^{\prime})$ and $J=(j,m,n)$ are multi-indices. To
implement the piecewise constant test functions $\chi_{n^{\prime}}$ as in [7]
instead, we simply exchange $B_{n^{\prime}}$ for $\chi_{n^{\prime}}$ in (3).
In addition, the notation $D_{i}(s_{j})$ represents the length of the
trajectory between $s_{j}\in\Gamma_{j}$ and
$s_{i}^{\prime}(s_{j},\tilde{p}_{n}(s_{j}))\in\Gamma_{i}$. We remark that the
four integrals required for the phase-space Galerkin projection have been
simplified to just one integral over $E^{j}_{m}\subset\Gamma_{j}$. Note that
in the third line of (3) we have made the additional assumption that the
reflection/transmission coefficients
$w_{i,j}(s_{j},\tilde{p}_{n}(s_{j}))=w_{i,j}(\tilde{p}_{n})$ depend only on
the arrival direction of the ray along an interface and not on the position
along the interface. For homogeneous polygonal sub-domains, the Euclidean
distance function $D_{i}(s_{j})$ is linear in $s_{j}\in E^{j}_{m}$ and hence
the one integral remaining in the third line of (3) can be evaluated
relatively simply. One can also apply analytic spatial integration for higher
order (Legendre polynomial) spatial basis functions as detailed in [11].

The expansion coefficients in (5) can now be obtained by solving a linear
system

${\bm{\rho}}=(I-B)^{-1}{\bm{\rho}}_{0}$

corresponding to the discretized operator equation (3). Here, the vectors
$\bm{\rho}_{0}$ and $\bm{\rho}$ contain the expansion coefficients for
$\rho_{0}$ and $\rho$, respectively. The source vector $\bm{\rho}_{0}$ can be
evaluated from the prescribed initial density $\rho_{0}$ by making use of the
orthonormality properties of the bases in both position and momentum as
follows

$[\bm{\rho}_{0}]_{J}=\int_{\Gamma_{j}\times(-\eta_{j},\eta_{j})}\rho_{0}(s_{j},p_{j})b_{m}(s_{j})B_{n}(p_{j})\mathrm{d}Y_{j}=\frac{\eta_{j}}{|E^{j}_{m}|^{1/2}}\int_{E^{j}_{m}}\int_{\phi_{n-1}}^{\phi_{n+1}}\rho_{0}(s_{j},p_{j}(\theta_{j}))\tilde{B}_{n}(\theta_{j})\cos(\theta_{j})\mathrm{d}\theta_{j}\mathrm{d}s_{j}.$

The corresponding expression when using the test functions $\chi_{n}$ instead
of $B_{n}$ is slightly simpler since the piecewise constants are independent
of the integration variable and hence

$[\bm{\rho}_{0}]_{J}=\int_{\Gamma_{j}\times(-\eta_{j},\eta_{j})}\rho_{0}(s_{j},p_{j})b_{m}(s_{j})\chi_{n}(p_{j})\mathrm{d}Y_{j}=\frac{\eta_{j}}{|E^{j}_{m}|^{1/2}}\int_{E^{j}_{m}}\int_{(\phi_{n-1}+\phi_{n})/2}^{(\phi_{n}+\phi_{n+1})/2}\rho_{0}(s_{j},p_{j}(\theta_{j}))\cos(\theta_{j})\mathrm{d}\theta_{j}\mathrm{d}s_{j}.$

Numerical results comparing these discretization schemes will be presented
during the conference and included in an extended version of this paper.

##  4 ACKNOWLEDGMENTS

Support from the EPSRC (grant no. EP/R012008/1) is gratefully acknowledged.

## References

  * Chappell  _et al._ [2013] D. Chappell, G. Tanner, D. Löchel, and N. Søndergaard, P. Roy. Soc. Lond. A Mat 469 (2013).
  * Hartmann  _et al._ [2019] T. Hartmann, S. Morita, G. Tanner, and D. Chappell, Wave Motion 87, 132 – 150 (2019).
  * Tanner [2009] G. Tanner, J. Sound Vib. 320, 1023 – 1038 (2009).
  * Chappell and Tanner [2014] D. Chappell and G. Tanner, Chaos 24, p. 043137 (2014).
  * Bajars and Chappell [2018] J. Bajars and D. J. Chappell, Journal of Computational Physics 373, 130–147 (2018).
  * Bajars and Chappell [2020] J. Bajars and D. J. Chappell, Commun Nonlinear Sci 80, p. 104973 (2020).
  * Chappell  _et al._ [2021] D. Chappell, J. Crofts, M. Richter, and G. Tanner, SIAM J. Sci. Comput. 43, B884–B906 (2021).
  * Lusikka, Ruotsalainen, and Saranen [1986] I. Lusikka, K. Ruotsalainen, and J. Saranen, Engineering Analysis 3, 144–153 (1986).
  * Ruotsalainen and Saranen [1987] K. Ruotsalainen and J. Saranen, SIAM J. Numer. Anal. 24, 816–827 (1987).
  * Ruotsalainen and Saranen [1988] K. Ruotsalainen and J. Saranen, Math. Meth. Appl. Sci. 10, 439–445 (1988).
  * Bajars  _et al._ [2017] J. Bajars, D. Chappell, T. Hartmann, and G. Tanner, J. Sci Comput 72, 1290–1312 (2017).


# Comparative analysis of the yields of dicentrics and chromosomal
translocations

Dorota Młynarczyk  Pedro Puig  Departament de Matemàtiques, Universitat
Autònoma de Barcelona, Bellaterra (Spain)  Centre de Recerca Matemàtica,
Bellaterra (Spain)  Joan F. Barquinero  Departament de Biologia Animal,
Biologia Vegetal i Ecologia, Universitat Autònoma de Barcelona, Bellaterra
(Spain)  Carmen Armero  Departament d’Estadística i Investigació Operativa,
Universitat de València, València (Spain)  Virgilio Gómez-Rubio  Department of
Mathematics, School of Industrial Engineering, Universidad de Castilla-La
Mancha, Albacete (Spain)

## Abstract

### Purpose

Chromosomal dicentrics and translocations are commonly employed as biomarkers
to estimate radiation doses. The main goal of this article is to perform a
comparative analysis of yields of both types of aberrations. The objective is
to determine if there are relevant distinctions between both yields, allowing
for a comprehensive assessment of their respective suitability and accuracy in
the estimation of radiation doses.

### Materials and Methods

The analysis involved data from a partial-radiation simulation study with the
calibration data obtained through two scoring methods: conventional and PAINT
modified. Subsequently, a Bayesian bivariate zero-inflated Poisson model was
employed to compare the posterior marginal density of the mean of dicentrics
and translocations and assess the differences between them.

### Results

When employing the conventional method of scoring, the findings indicate that
there is no notable disparity between the yield of observed translocations and
dicentrics. However, when utilizing the PAINT modified method, a notable
discrepancy is observed for higher doses, indicating a relevant difference in
the mean number of the two types aberrations.

### Conclusions

The choice of scoring method significantly influences the analysis of
radiation-induced aberrations, especially when distinguishing between complex
and simple chromosomal formations. Further research and analysis are necessary
to gain a deeper understanding of the factors and mechanisms impacting the
formation of dicentrics and translocations.

### Keywords

bivariate zero-inflated Poisson model, FISH-painting, PAINT and conventional
nomenclatures, radiation-induced chromosome aberrations

##  1 Introduction

The quantitative assessment of dicentrics and translocations, two well-studied
chromosomal aberrations induced by exposure to ionizing radiation, plays a
pivotal role in radiation biology. Sometimes it is assumed that yields for
these two types of aberrations are very similar, but concerns may arise
regarding their potential inequality, highlighting the complex mechanisms that
influence these chromosomal changes in response to varying radiation doses.
The current radiation biodosimetry manual ([IAEA 2011]) recommends that each
laboratory should develop its own fitting dose-response curves for
translocations, but, in practice, in situations where a laboratory lacks a
specific translocations curve, it is occasionally acceptable to use the
calibration curve for dicentrics ([Barquinero et al. 2017]). Furthermore, the
analysis of aberrations may also be influenced by the detection technique
employed and the terminology used to characterize these alterations
([Barquinero et al. 1999]). Thus, the objective of this article is to compare
yields for dicentrics and translocations obtained using different scoring
systems, with the aim of evaluating the technique’s influence on the
disparities between these yields. While some authors have previously tackled
this issue ([Barquinero et al. 1999]; [Lindholm et al. 1998]), our current
article emphasizes a new mathematical approach.

Radiation-induced chromosome exchange aberrations produced in lymphocytes in
G0 (quiescent) stage have been classified as symmetrical or asymmetrical (i.e
reciprocal translocations or dicentrics) ([Savage and Papworth 1982]). For
biological dosimetry purposes and in solid stained metaphases, symmetrical
exchanges are difficult to detect, unless the resulting chromosomes are
markedly different from the normal karyotype. For this reason, initially
biological dosimetry was based on the detection of dicentrics or dicentrics
plus rings ([IAEA 2001]). From the visualization of a solid stained dicentric
plus its corresponding acentric fragment it was logical to infer that this
exchange resulted to the interaction of two broken chromosomes, and by
association the same for symmetric exchanges like translocations. The
introduction of fluorescence in situ hybridization techniques (FISH) to detect
whole chromosomes (chromosome-painting) allowed an easy detection of
translocations ([Lucas et al. 1992]). However, with the same technique, it was
evident that radiation-induced chromosome exchanges could be formed by the
interaction of more than two breaks ([Cremer 1990]; [Schmid et al. 1992]) and
exchanges were classified between simple and complex. The former involve two
breaks and the latter involving three or more breaks in two or more
chromosomes ([Savage and Simpson 1994]).

The first calibration curves produced by chromosome-painting used the terms
translocation and dicentric similarly to that was done when analyzing solid
stained metaphases ([Lucas et al. 1992, Bauchinger et al. 1993, Fernández et
al. 1995]), i.e. distinguishing aberrations: translocation(t), dicentric(dic),
ring(r), insertion(ins), and additional acentric(ace) fragment ([ISCN 1985]).
This system, hereafter referred to as the conventional system, may also
identify complete or incomplete dicentrics and translocations, along with
centric or acentric rings. However, the presence of complex exchanges mainly
at higher doses led the development of specific nomenclatures to describe the
radiation-induced exchanges. Initially, two highly different nomenclatures
were proposed: the nomenclature with the acronym PAINT (from Protocol for
Aberration Identification and Nomenclature Terminology; [Tucker et al. 1995])
was purely descriptive of each abnormal piece present in the metaphase,
without cross-reference to other aberrant pieces in the same metaphase; and
the CAB terminology (from Chromosomes-Arms-Breaks; [Cornforth 2001]) or the
so-called S&S system ([Savage and Simpson 1994]) that proposed a code with
numerals and letters that refer to the number of abnormal pieces observed and
how common the exchange was expected to be. Understanding complex aberrations
is greatly aided by the CAB terminology, which focused on the mechanical
aspects of exchange formation. However, in daily practice the use of the CAB
vocabulary was not easy to handle, and currently the most widely used method
to score the chromosome exchanges using FISH is to describe each abnormal
metaphase as a unit using the PAINT nomenclature in a slightly modified way
that allows to infer the mechanistic aspects of exchange formation and allows
classifying between simple and complex aberrations ([Knehr et al. 1998]).
Specifically, within the PAINT modified system, terms like ’apparently simple
dicentric/translocation’ (ASD/AST) are employed to describe aberrations
involving bicolored chromosomes with single-color junctions each ([Tucker et
al. 1995]), which may be the result of complex, undetected aberrations.

Figure 1:  An example of a metaphase illustrating FISH-based chromosome
‘painting’. In red painted chromosomes #1, #4, and #11, in blue counterstained
chromosomes.

Despite the existence of these nomenclatures, some biological dosimetry
laboratories still rely on conventional terminology to describe chromosomal
alterations. Therefore, the question is how different scoring systems
(conventional vs. PAINT modified) influence the analysis of yields of
aberrations, and whether these can lead to misleading interpretations of the
results. To illustrate the problem, let’s consider an example from Figure 1.
The metaphase contains four aberrant pieces with painted material involved.
The metaphase would be scored as: dic(BA), t(Ba), ace(ab), t(Ab). These
aberrations could be scored by combining aberrations that imply the minimum
number of breaks, i.e., dic(BA) with ace(ab), which is an ASD or a dicentric,
and t(Ba) with t(Ab), which is an AST or a translocation. However, one could
also combine dic(BA) with t(Ab), resulting in a complex aberration (neither
ASD nor AST), and it will be counted as one dicentric plus one translocation
in the conventional nomenclature. Similarly, combining t(Ba) with ace(ab)
results in a complex aberration, which will be counted as a translocation.

For the purpose of this research, we decided to utilize the data obtained from
the study conducted by [Duran et al.] (2002). The study involved a blood
sample taken from a healthy 32-year-old male that were exposed to X-rays at
radiation doses of 2, 3, 4, and 5 Gy within a laboratory environment.
Subsequent to the exposure, the samples underwent the FISH protocol with
chromosomes 1, 4, and 11 being painted. The scoring of chromosomal aberrations
was carried out using both the conventional and PAINT modified methods. This
implies that we have counts of dicentrics and translocations via the
conventional scoring approach, as well as counts of ASD and AST using the
PAINT modified approach for each cell, what enables a comparison between the
results obtained from both methods.

[Duran et al.] (2002) studied specifically partial body radiation exposures.
This type of exposure occurs when only a specific region or part of the body
is subjected to radiation, such as during targeted radiation therapy or
accidental exposure to a localized radiation source. To replicate this
situation experimentally, various proportions of irradiated and non-irradiated
blood were mixed together: 0.875, 0.75, 0.5, 0.25, and 0.125. Furthermore, the
analysis also considered the impact of total irradiation by including samples
with a proportion of 1 (without the addition of non-irradiated blood to the
samples). Subsequently, these mixed samples were carefully scored using both
considered terminologies.

| dose (Gy) | 2 | 2 | 2 | 2 | 2 | 2 | 3 | 3 | 3 | 3 | 3 | 3
---|---|---|---|---|---|---|---|---|---|---|---|---|---
| proportion of irradiated cells $p$ | 1 | 0.875 | 0.75 | 0.50 | 0.25 | 0.125 | 1 | 0.875 | 0.75 | 0.50 | 0.25 | 0.125
| number of cells | 525 | 974 | 1551 | 1322 | 1516 | 1322 | 394 | 509 | 824 | 1009 | 1070 | 1096
conventional | mean dicentrics | 0.12 | 0.07 | 0.06 | 0.04 | 0.02 | 0.01 | 0.21 | 0.22 | 0.10 | 0.08 | 0.03 | 0.02
| mean translocations | 0.11 | 0.08 | 0.06 | 0.05 | 0.02 | 0.01 | 0.27 | 0.21 | 0.10 | 0.10 | 0.03 | 0.04
PAINT mod. | mean ASD | 0.10 | 0.07 | 0.06 | 0.04 | 0.02 | 0.01 | 0.16 | 0.21 | 0.10 | 0.07 | 0.03 | 0.02
| mean AST | 0.08 | 0.07 | 0.05 | 0.04 | 0.02 | 0.01 | 0.20 | 0.17 | 0.08 | 0.07 | 0.02 | 0.03
| dose (Gy) | 4 | 4 | 4 | 4 | 4 | 4 | 5 | 5 | 5 | 5 | 5 | 5
| proportion of irradiated cells $p$ | 1 | 0.875 | 0.75 | 0.50 | 0.25 | 0.125 | 1 | 0.875 | 0.75 | 0.50 | 0.25 | 0.125
| number of cells | 250 | 463 | 504 | 775 | 1035 | 1040 | 133 | 374 | 553 | 523 | 890 | 1260
conventional | mean dicentrics | 0.48 | 0.23 | 0.16 | 0.10 | 0.04 | 0.02 | 0.47 | 0.35 | 0.20 | 0.11 | 0.06 | 0.02
| mean translocations | 0.50 | 0.20 | 0.19 | 0.12 | 0.04 | 0.02 | 0.59 | 0.30 | 0.18 | 0.13 | 0.06 | 0.02
PAINT mod. | mean ASD | 0.38 | 0.20 | 0.14 | 0.08 | 0.04 | 0.01 | 0.33 | 0.27 | 0.18 | 0.08 | 0.05 | 0.02
| mean AST | 0.33 | 0.16 | 0.14 | 0.08 | 0.03 | 0.02 | 0.42 | 0.20 | 0.13 | 0.09 | 0.04 | 0.01

Table 1: Mean number of aberrations for each dose and dilution $p$, utilizing
both the conventional and PAINT mod. scoring systems.

Table 1 displays the mean number of aberrations for each dose and dilution $p$
(proportion of irradiated cells), utilizing both the conventional and PAINT
mod. nomenclature. In total, 20,332 cells were scored, with the lowest count
of 3,733 at a dose of 5Gy and of 1,302 for proportion $p=1$. This table
highlights relevant differences between the two scoring systems under
consideration as evident in the variations in the mean numbers of
dicentrics/ASD and translocations/AST. The raw data can be found in the
supplementary material of this article.

##  2 Materials and methods

From a statistical standpoint, the bivariate Poisson model seems to be a
suitable choice when working with both dicentrics and translocations, as it
allows for joint analysis and inference about the underlying processes. This
model assumes that both count variables are generated by two different Poisson
processes, with a specific correlation term between them. Moreover, the
bivariate Poisson model can be extended to cover scenarios with partial body
exposures. The experimental design, which combines irradiated and non-
irradiated samples, results in a significant number of zero counts in the
observed data, which can be effectively handled by a bivariate zero-inflated
Poisson model. While bivariate Poisson models have been widely employed for
analyzing sports and healthcare data ([Karlis and Ntzoufras 2003]), their
application in the field of biodosimetry has been largely unexplored.
Furthermore, the bivariate approach utilized in this study for analyzing
dicentrics and translocations can also be applied to analyze dicentrics and
rings or dicentrics and foci. However, it would be more difficult in the case
of foci due to their significant dependence on the amount of time that has
passed following irradiation ([Młynarczyk et al. 2022]).

To begin the probabilistic formulation of the sampling model, let’s clarify
that we are defining a single model for both scoring techniques. Therefore, in
this section, we will simply refer to aberrations as dicentrics and
translocations, without mentioning ASD or AST. Assume that $N$ cells were
examined. We denote the number of dicentrics $X_{j}$ and of translocations
$Y_{j}$ found in cell $j$, $j=1,\ldots,N$. Then, it is assumed that $X_{j}$
and $Y_{j}$ jointly follow a conditional zero-inflated bivariate Poisson
distribution, with probability function

$f_{ZIBP}(x_{j},y_{j}\mid\omega,\lambda_{1j},\lambda_{2j},\lambda_{3j})=\begin{cases}\omega
f_{BP}(x_{j},y_{j}\mid\lambda_{1j},\lambda_{2j},\lambda_{3j})&(x_{j},y_{j})\neq(0,0)\\\
(1-\omega)+\omega
f_{BP}(x_{j},y_{j}\mid\lambda_{1j},\lambda_{2j},\lambda_{3j})&(x_{j},y_{j})=(0,0)\\\
\end{cases}$ (1)

where $(1-\omega)$ is the proportion of structural zeros (proportion of not-
irradiated cells), and
$f_{BP}(x_{j},y_{j}\mid\lambda_{1j},\lambda_{2j},\lambda_{3j})$ is the
conditional probability function of the bivariate Poisson distribution given
by

$f_{BP}(x_{j},y_{j}\mid\lambda_{1j},\lambda_{2j},\lambda_{3j})=\exp(-(\lambda_{1j}+\lambda_{2j}+\lambda_{3j}))\frac{\lambda_{1j}^{x_{j}}}{x_{j}!}\frac{\lambda_{2j}^{y_{j}}}{y_{j}!}\sum_{k=0}^{\min(x_{j},y_{j})}\binom{x_{j}}{k}\binom{y_{j}}{k}k!(\frac{\lambda_{3j}}{\lambda_{1j}\lambda_{2j}})^{k}.$

The effect on the irradiated cells is determined by the parameters
$\lambda_{1j},\lambda_{2j}$, and $\lambda_{3j}$ in our model. In more detail,
$\lambda_{1j}+\lambda_{3j}$ represents the marginal mean of the number of
dicentrics, $\lambda_{2j}+\lambda_{3j}$ the marginal mean of translocations,
and $\lambda_{3j}$ the covariance term between the number of dicentrics and
translocations (all them for the irradiated cells). These parameters, all them
greater than zero, can be modelled as functions of the radiation dose using a
regression approach. The mean number of aberrations is frequently estimated in
biodosimetric research using a linear or quadratic function of dose; we chose
to concentrate on the quadratic model because the linear model is also
included in it. We specify them as,

$\begin{aligned} \text{dicentrics}\\\ \text{translocations}\\\
\text{covariance}\\\ \end{aligned}\quad\begin{aligned}
\lambda_{1j}&=\beta_{11}+\beta_{12}\cdot d_{j}+\beta_{13}\cdot d^{2}_{j},\\\
\lambda_{2j}&=\beta_{21}+\beta_{22}\cdot d_{j}+\beta_{23}\cdot d^{2}_{j},\\\
\lambda_{3j}&=\beta_{31},\end{aligned}$ (2)

where $d_{j}$ denotes the radiation dose received by $j$-th cell, and
$\beta_{\bullet\bullet}$ denotes the corresponding regression coefficients.
The expected number of dicentrics, including irradiated and not-irradiated
cells, is
$\mu_{1j}=E(X_{j}\mid\beta_{11},\beta_{12},\beta_{13})=\omega(\lambda_{1j}+\lambda_{3j})$.
Similarly for translocations
$\mu_{2j}=E(Y_{j}\mid\beta_{21},\beta_{22},\beta_{23})=\omega(\lambda_{2j}+\lambda_{3j})$.

It has been studied that some cells exposed to radiation fail to survive until
the metaphase stage of the cell cycle, when dicentrics and translocations are
analyzed ([Hall and Giaccia 2006]). As a result, the actual proportion of
irradiated cells present in the sample at the time of the analysis is lower
than the initially assumed proportion. This discrepancy between both
proportions poses a challenge in accurately assessing the radiation-induced
chromosomal aberrations. Therefore, it is essential to account for this
limitation, so researchers employ a statistical correction to estimate the
true proportion of irradiated cells ([Pujol et al. 2016]). The survival rate,
the proportion of the irradiated cells which survive until the metaphase, is
described as a decreasing exponential function of the dose,

$s(d)=\text{exp}(-\gamma d),$

where $\gamma$ is the so-called survival index, $\gamma\in[0,1]$. Therefore
the actual proportion $\omega$ of scored cells irradiated at dose $d$ is given
by (see supplementary material for details)

$\omega(d,p)=\frac{1}{1+\text{exp}(\gamma d)(1-p)/p},$ (3)

where $p$ is the initial proportion of the irradiated cells. Note that in the
given experiment $d$ and $p$ are known, but the survival index $\gamma$ needs
to be estimated.

Let denote by $\mathbf{x}=(x_{1},\ldots,x_{N})$ the counts of dicentrics in
cells $1,\ldots,N$, and respectively the counts of translocations by
$\mathbf{y}=(y_{1},\ldots,y_{N})$. The vector
$\boldsymbol{d}=(d_{1},\ldots,d_{N})$ stands for the radiation doses to which
each cell was exposed and $\boldsymbol{p}=(p_{1},\ldots,p_{N})$ is the vector
of initial proportion of irradiated cells in the sample from which the
observation comes. Thus the likelihood of the model is given by

$\displaystyle
L(\mathbf{x},\mathbf{y},\boldsymbol{d},\boldsymbol{p}\mid\boldsymbol{\beta},\gamma)$
$\displaystyle=\prod_{j=1}^{N}\Biggr{(}1_{(0,0)}(x_{j},y_{j})\left(1-\frac{1}{1+\text{exp}(\gamma
d_{j})(1-p_{j})/p_{j}}\right)$ $\displaystyle+\frac{1}{1+\text{exp}(\gamma
d_{j})(1-p_{j})/p_{j}}\
\text{exp}(-(\lambda_{1j}+\lambda_{2j}+\lambda_{3j}))\frac{\lambda_{1j}^{x_{j}}}{x_{j}!}\frac{\lambda_{2j}^{y_{j}}}{y_{j}!}$
$\displaystyle\times\sum_{k=0}^{\min(x_{j},y_{j})}\binom{x_{j}}{k}\binom{y_{j}}{k}k!(\frac{\lambda_{3j}}{\lambda_{1j}\lambda_{2j}})^{k}\Biggr{)},$

where $\lambda_{ij}$ are defined in (2), and $1_{(0,0)}(x_{j},y_{j})$ is an
indicator function that is 1 when $x_{j}=y_{j}=0$ and 0 otherwise. Within the
Bayesian framework, the main interest is in computing the posterior
distribution
$\pi(\boldsymbol{\beta},\gamma\mid\mathbf{x},\mathbf{y},\boldsymbol{d},\boldsymbol{p})$
of the parameters of the model, i.e. the regression coefficients
$\boldsymbol{\beta}=(\beta_{11},\ldots,\beta_{31})$ and the survival index
$\gamma$, through the Bayes´ theorem

$\pi(\boldsymbol{\beta},\gamma\mid\mathbf{x},\mathbf{y},\boldsymbol{d},\boldsymbol{p})\propto
L(\mathbf{x},\mathbf{y},\boldsymbol{d},\boldsymbol{p}\mid\boldsymbol{\beta},\gamma)\pi(\boldsymbol{\beta},\gamma),$

where $\pi(\boldsymbol{\beta},\gamma)$ is the prior distribution for
$(\boldsymbol{\beta},\gamma)$. We assume prior independence between
$\boldsymbol{\beta}$ and $\gamma$, so
$\pi(\boldsymbol{\beta},\gamma)=\pi(\boldsymbol{\beta})\pi(\gamma)$. Given the
fact that $\lambda_{1j}$ and $\lambda_{2j}$ should be non-negative, the prior
distribution $\pi(\boldsymbol{\beta})$ was chosen to be a non-informative
uniform distribution on the interval $(0,10)$. The prior distribution for
$\gamma$ was chosen to be a uniform distribution on the interval $[0,1]$
because, according to the definition, the survival index belongs to the
interval $[0,1]$. Approximated samples from the posterior distribution
$\pi(\boldsymbol{\beta},\gamma\mid\mathbf{x},\mathbf{y},\boldsymbol{d},\boldsymbol{p})$
were generated using Markov chain Monte Carlo methods. The JAGS software
([Plummer 2003]), in particular, was utilized with specific configurations,
including the use of 2 chains, conducting 25,000 simulations, burn-in of 1000
iterations, and thinning every 25 iterations (the program itself is available
in the supplementary material).

Note that $\lambda_{1j}$ and $\lambda_{2j}$ depend on the dose received by the
cell $j$ and are determined by equations (2). By fixing the dose and utilizing
samples from the posterior distribution, it becomes straightforward to derive
the posterior distribution for the expected number of aberrations. From this
point onward, we will focus only on the irradiated cells (with the actual
proportion $\omega(d,p)=1$), since we regard this as the most interesting
case. For a fixed dose $d$, we can omit the subindex $j$ in equations (2), and
denote the mean number of dicentrics of the irradiated cells for dose $d$ as
$\mu_{1d}$ and of translocations as $\mu_{2d}$. Our main goal is to discuss
whether these mean values differ from each other. To do so, we can consider
the posterior distribution of the difference between means
$\pi(\mu_{1d}-\mu_{2d}\mid\mathbf{x},\mathbf{y},\boldsymbol{d},\boldsymbol{p})$.
It can be easily approximated through simulation, enabling the calculation of
the posterior mean of the difference, credible intervals, and conducting
graphical analysis by plotting posterior densities for each dose. We will
report Highest Density Intervals (HDI), which are credible intervals designed
to have a higher probability density for all values inside the interval
compared to those outside ([Kruschke 2014]).

We can additionally compare the models by assessing the probability that the
mean number of dicentrics of the irradiated cells, $\mu_{1d}$, is lower than
the mean number of translocations, $\mu_{2d}$, separately calculated for each
dose $d$. Let’s define this as hypothesis $H_{1}:\mu_{1d}-\mu_{2d}<0$.
Conversely, the alternative scenario will be denoted as $H_{2}$, with
$H_{2}:\mu_{1d}-\mu_{2d}\geq 0$. Evaluating the ratio between the posterior
probability of $H_{1}$ and that of $H_{2}$, i.e.

$\frac{\pi(\mu_{1d}-\mu_{2d}<0\mid\mathbf{x},\mathbf{y},\boldsymbol{d},\boldsymbol{p})}{\pi(\mu_{1d}-\mu_{2d}\geq
0\mid\mathbf{x},\mathbf{y},\boldsymbol{d},\boldsymbol{p})}$ (4)

may assist us in determining which situation is more likely. Given our prior
assumption that both hypotheses are equally probable a priori, i.e.,
$P(H_{1})=P(H_{2})=0.5$, this resulting odds can be interpreted as the Bayes
factor ([Kass and Raftery 1995]) (details can be found in supplementary
material). According to the Jeffreys scale ([Jeffreys 1961]), a Bayes factor
close to 1 suggests only ’bare mention’ evidence of a difference between
$\mu_{1d}$ and $\mu_{2d}$. A number greater than 3 but less than 10 is
considered ’substantial’ evidence. Results between 10 and 30 indicate ’strong’
evidence, and between 30 and 100 are categorized as ’very strong.’ A Bayes
factor exceeding 100 is considered ’decisive’.

##  3 Results

Our objective is to examine the feasibility of employing the same calibration
curve for both dicentrics and translocations, which entails comparing the
parameters $\mu_{1d}$ and $\mu_{2d}$ defined in Section 2. Figure 2 depicts
the posterior distribution of the means of the number of dicentrics $\mu_{1d}$
and of translocations $\mu_{1d}$ for conventional scoring method and for ASD
and AST for PAINT modified technique, across different doses. Notably, this
information reveals that the posterior mean values of dicentrics and
translocations are relatively similar when using conventional nomenclature.
The density of dicentrics of the irradiated cells is skewed to the left for
each dose compared to translocations. Conversely, in the case of the PAINT
modified scoring, the posterior distribution for AST is shifted towards the
left. This difference becomes more pronounced as the dose increases,
indicating a relevant difference between them, what can also be noticed
looking at the Figure 3. This figure represents the posterior density of the
difference between estimated means of dicentrics and translocations (or ASD
and AST),
$\pi(\mu_{1d}-\mu_{2d}\mid\mathbf{x},\mathbf{y},\boldsymbol{d},\boldsymbol{p})$
for the irradiated cells at each dose and for both scoring techniques.
Posterior mean and $95\%$ HDI intervals of the differences for each case can
be found in the Table 2. All of these findings indicate that there are some
distinctions in the means, especially in the case of a dose of 4Gy for the
PAINT mod. method. Nevertheless, it’s important to mention that more
noticeable differences can be detected at higher doses.

Figure 2:  Posterior densities of the mean number of dicentrics,
$\pi(\mu_{1d}\mid\mathbf{x},\mathbf{y},\boldsymbol{d},\boldsymbol{p})$, and of
the mean number of translocations,
$\pi(\mu_{2d}\mid\mathbf{x},\mathbf{y},\boldsymbol{d},\boldsymbol{p})$, given
by the model for two scoring techniques. Figure 3:  Posterior distribution of
the difference between means of dicentrics and translocations
$\pi(\mu_{1d}-\mu_{2d}\mid\mathbf{x},\mathbf{y},\boldsymbol{d},\boldsymbol{p})$
given by the model for two scoring techniques. The red line is a reference
line when $\mu_{1d}-\mu_{2d}=0$.

| conventional | PAINT mod.
---|---|---
dose | 2 | 3 | 4 | 5 | 2 | 3 | 4 | 5
mean | -0.006 | -0.008 | -0.01 | -0.011 | 0.001 | 0.009 | 0.022 | 0.041
HDI | [-0.019,0.008] | [-0.022,0.008] | [-0.035,0.017] | [-0.06,0.039] | [-0.013,0.015] | [-0.004,0.023] | [0.001,0.044] | [-0.001,0.085]

Table 2: Posterior mean of the difference between means of dicentrics and
translocations and $95\%$ HDI for each dose and scoring technique.

The results presented in Table 3 showcase the posterior odds, defined
previously (4), as the ratio of probabilities between instances where
$\mu_{1d}$ is lower than $\mu_{2d}$ and vice versa. However, for a more robust
interpretation of the results, for the PAINT mod. method, the Bayes factor was
computed as the reciprocal of this fraction. According to the interpretation
provided by ([Jeffreys 1961]), the results offer substantial evidence
supporting the hypothesis that $\mu_{1d}$ is lower than $\mu_{2d}$ for the
conventional technique, a condition generally observed, as depicted in Figure
2. For the PAINT technique, when the mean of ASD exceeds the mean of AST,
strong evidence is found indicating that $\mu_{1d}$ is greater than or equal
to $\mu_{2d}$ for doses of 4Gy and 5Gy.

dose (Gy) | 2 | 3 | 4 | 5
---|---|---|---|---
conventional | 3.940 | 6.369 | 3.446 | 1.941
PAINT mod. | 1.301 | 8.340 | 46.729 | 33.333
Table 3: Results of Bayes factor as defined by ratio (4) for the conventional
method of scoring and computed as the reciprocal of the fraction for the PAINT
mod. method.  Figure 4:  Posterior density of the survival index, $\gamma$,
for the two methods of scoring.

The current information about the survival index $\gamma$ of the irradiated
cells is represented by the posterior distribution in Figure 4. It is
important to note that this survival rate, as defined in Section 2, is dose-
dependent. Higher doses of radiation generally correspond to lower survival
proportions, while lower doses tend to result in higher survival proportions.
Additionally, the actual proportion of irradiated cells is also influenced by
the initial proportion of irradiated cells as can be seen by (3). Table 4
illustrates the relationship between the initial proportion of cells p, the
received dose of radiation d, and the resulting actual proportion of surviving
cells $\omega$. These values are determined based on the mean value of the
posterior density of the survival index $\gamma$. The estimated mean value of
$\gamma$ for the PAINT mod. method was 0.19, while for the conventional method
it was 0.22. These results indicate a significant decrease in the proportion
of irradiated cells that survived until metaphase. For example, when using the
conventional method with a dose of 5Gy and an initial proportion of 0.5, the
estimated actual proportion of irradiated cells in the analyzed sample is
approximately half of the initial proportion.

|  | Initial proportion of irradiated cells $p$
---|---|---
Scoring method | dose $d$ (Gy) | 1 | 0.875 | 0.75 | 0.5 | 0.25 | 0.125
conventional | 2 | 1 | 0.818 | 0.658 | 0.39 | 0.176 | 0.084
| 3 | 1 | 0.782 | 0.606 | 0.339 | 0.146 | 0.068
| 4 | 1 | 0.742 | 0.552 | 0.291 | 0.120 | 0.055
| 5 | 1 | 0.697 | 0.496 | 0.247 | 0.099 | 0.045
PAINT modified | 2 | 1 | 0.827 | 0.672 | 0.406 | 0.185 | 0.089
| 3 | 1 | 0.798 | 0.629 | 0.361 | 0.158 | 0.075
| 4 | 1 | 0.766 | 0.583 | 0.318 | 0.135 | 0.062
| 5 | 1 | 0.730 | 0.536 | 0.278 | 0.114 | 0.052
Table 4: Actual proportion of irradiated cells $\omega$ estimated by the
posterior mean of the survival $\gamma$ given the dose of radiation $d$ (in
Gy) and the initial proportion $p$.

##  4 Discussion

As expected, the results of this study suggest that the choice of the scoring
method significantly influences the analysis of the radiation-induced
aberrations. However, the differences become less evident when the focus is on
yields of aberrations using the same scoring technique. When using the
conventional scoring method, the results suggest that there is a very small
difference between the mean numbers of translocations and dicentrics observed.
However, with the PAINT mod. method, a more pronounced variation in the
results was found, but it’s important to emphasize that our focus was solely
on apparently simple aberrations. We believe that these results are primarily
attributable to the fact that the conventional method does not differentiate
between complex and simple chromosomal formations, whereas the PAINT mod.
method does.

The prevailing understanding in past studies suggests that dicentrics and
translocations have equal probability of appearing after chromosomal breaks
induced by ionizing radiation ([Lucas 1996]). However, it is important to note
that there is no unanimous consensus of an equal formation ratio (1:1) and
this theoretical assumption is occasionally questioned ([Barquinero et al.
1998]). Based on the findings from the PAINT mod. method, our results indicate
that, when complex aberrations are excluded, there is a noteworthy difference
in the mean number of dicentrics and translocations, particularly at higher
radiation doses. This result suggests that the 1:1 ratio may not be valid; if
this is the case, the results may imply that translocations would occur more
frequently than dicentrics in complex aberrations.

However, another explanation of the results for PAINT mod. method, when the
yield of AST is shifted left comparing to ASD, may be the limit of detection
of exchanged chromosome fragments. It has been described that using FISH
techniques the minimum detectable size is approximately 11 and 14 Mb for
painted and unpainted material, respectively ([Kodama et al. 1997]). ASTs in
which small fragments are exchanged may go unnoticed. On the contrary, if two
chromosomes break at their terminal part and form a dicentric chromosome, it
will be clearly visible.

Another possibility is that the probability to form an ADS or an AST changes
as the dose increase. Exchange type aberrations have a fast kinetics formation
([Darroudi et al. 1998]), and using fussion PCC techniques to evaluate the
formation of dicentrics it was described that after 2 Gy irradiation the
dicentrics were mostly formed during the first two hours post irradiation,
while the resting unsolved damage seemed to reconstitute the original
chromosomes ([Pujol et al. 2020]), it can be hypothetized that at higher doses
the probability to form an ASD increases with respect the probability to form
a AST. Consequently, further investigation and analysis are necessary to gain
a better understanding of the underlying mechanisms and factors that influence
the formation of dicentrics and translocations. The results observed in the
present study have an implication in biological dosimetry studies. For past
dose assessment, it is important to choose a specific biomarker, whether using
the conventional nomenclature or PAINT modified one, and that the results
obtained at high doses are not always interchangeable.

## Ethical approval

## Disclosure statement

No potential conflict of interest was reported by the authors.

## Author contributions

Conceptualization, D.M., P.P., JF.B.; Methodology, D.M., P.P., C.A. and
V.G.R.; Software, D.M.; Formal Analysis, D.M.; Resources, JF.B. ; Data
Curation, JF.B.; Writing–Original Draft, D.M. and JF.B.;
Writing–Review&Editing, D.M., P.P., C.A., V.G.R, JF.B.

## Data and code availability

All data and computer code used for analysis during this study are included in
supplementary material.

## Funding

This work was supported by the Consejería de Educación, Cultura y Deportes
(Junta de Comunidades de Castilla-La Mancha, Spain) [the Project MECESBayes,
SBPLY/17/180501/000491]; Ministerio de Ciencia e Innovación (Spain) [research
grants PID2019-106341GB-I00, PID2022-137414OB-I00]; the Spanish Consejo de
Seguridad Nuclear [BOE-A-2019-311]; and the Spanish State Research Agency [the
Severo Ochoa and Marıa de Maeztu Program for Centers and units of Excellence
in R&D (CEX2020-001084-M)].

## Notes on contributors

Dorota Mlynarczyk, PhD student, Mathematician, PhD student at the Department
of Mathematics, Faculty of Sciences, Universitat Autònoma de Barcelona,
Bellaterra (Catalonia), Spain.

Pedro Puig, PhD, Mathematician, Professor at the Department of Mathematics and
member of the Center de Recerca Matemàtica, Faculty of Sciences, Universitat
Autònoma de Barcelona, Bellaterra (Catalonia), Spain.

Joan-Francesc Barquinero, PhD, Biologist, University Professor at the
Department of Animal Biology, Plant Biology and Ecology, Faculty of
Biosciencies, Universitat Autònoma de Barcelona, Bellaterra (Catalonia),
Spain.

Carmen Armero, PhD, Mathematician, Professor at the Department of Statistics
and Operations Research, Faculty of Mathematics, Universitat de València,
Burjassot (Valencian Community), Spain.

Virgilio Gómez-Rubio, PhD, Mathematician, Professor at the Department of
Mathematics, School of Industrial Engineering, Universidad de Castilla-La
Mancha, Albacete (Castilla-La Mancha), Spain.

## References

  * [Bauchinger et al. 1993] Bauchinger M, Schmid E, Zitzelsberger H, Braselmann H, Nahrstedt U. 1993. Radiation-induced chromosome aberrations analysed by two-colour fluorescence in situ hybridization with composite whole chromosome-specific DNA probes and a pancentromeric DNA probe. Int J Radiat Biol. 64(2):179-84. doi: 10.1080/09553009314551271.
  * [Barquinero et al. 1998] Barquinero JF, Knehr S, Braselmann H, Figel M, Bauchinger M. 1998. DNA-proportional distribution of radiation-induced chromosome aberrations analysed by fluorescence in situ hybridization painting of all chromosomes of a human female karyotype. Int J Radiat Biol. 74(3):315-23. doi: 10.1080/095530098141456.
  * [Barquinero et al. 1999] Barquinero JF, Cigarrán S, Caballín MR, Braselmann H, Ribas M, Egozcue J, Barrios L. 1999. Comparison of X-ray dose-response curves obtained by chromosome painting using conventional and PAINT nomenclatures. Int J Radiat Biol. 75(12):1557-66. doi: 10.1080/095530099139160.
  * [Barquinero et al. 2017] Barquinero JF, Beinke C, Borràs M, Buraczewska I, Darroudi F, Gregoire E, Hristova R, Kulka U, Lindholm C, Moreno M, et al. 2017. RENEB biodosimetry intercomparison analyzing translocations by FISH. Int J Radiat Biol. 93(1):30-35. doi: 10.1080/09553002.2016.1222092.
  * [Cornforth 2001] Cornforth MN. Analyzing radiation-induced complex chromosome rearrangements by combinatorial painting. 2001. Radiat Res. 155(5):643-59. doi: 10.1667/0033-7587(2001)155[0643:ARICCR]2.0.CO;2.
  * [Cremer 1990] Cremer T, Popp S, Emmerich P, Lichter P, Cremer C. 1990. Rapid metaphase and interphase detection of radiation-induced chromosome aberrations in human lymphocytes by chromosomal suppression in situ hybridization. Cytometry. 11(1):110-8. doi: 10.1002/cyto.990110113.
  * [Darroudi et al. 1998] Darroudi F, Fomina J, Meijers M, Natarajan AT. 1998. Kinetics of the formation of chromosome aberrations in X-irradiated human lymphocytes, using PCC and FISH. Mutat Res. 404(1-2):55-65. doi: 10.1016/s0027-5107(98)00095-5.
  * [Duran et al.] Duran A, Barquinero JF, Caballín MR, Ribas M, Puig P, Egozcue J, Barrios L. 2002.Suitability of FISH painting techniques for the detection of partial-body irradiations for biological dosimetry. Radiat Res. 157(4):461-8. doi: 10.1667/0033-7587(2002)157[0461:sofptf]2.0.co;2.
  * [Fernández et al. 1995] Fernández JL, Campos A, Goyanes V, Losada C, Veiras C, Edwards AA. 1995. X-ray biological dosimetry performed by selective painting of human chromosomes 1 and 2. Int J Radiat Biol. 67(3):295-302. doi: 10.1080/09553009514550351.
  * [Hall and Giaccia 2006] Hall EJ, Giaccia AJ. 2006. Radiobiology for the Radiologist. 6th ed. Philadelphia (PA): Lippincott Williams and Wilkins.
  * [IAEA 2001] [IAEA] International Atomic Energy Agency. 2001. Cytogenetic Analysis for Radiation Dose Assessment. Austria: IAEA in Austria.
  * [IAEA 2011] [IAEA] International Atomic Energy Agency. 2011. Cytogenetic Dosimetry: Applications in Preparedness for and Response to Radiation Emergencies. Austria: IAEA in Austria.
  * [ISCN 1985] [ISCN] An International System for Human Cytogenetic Nomenclature. 1985. Report of the Standing Committee on Human Cytogenetic Nomenclature. Birth Defects Orig Artic Ser. 21(1):1-117.
  * [Plummer 2003] Plummer M. 2003. JAGS: A Program for Analysis of Bayesian Graphical Models Using Gibbs Sampling. Proceedings of the 3rd International Workshop on Distributed Statistical Computing (DSC 2003); Mar 20-22; Vienna, Austria.
  * [Jeffreys 1961] Jeffreys H. 1961. Theory of Probability. 3rd ed. Oxford (UK): Clarendon Press.
  * [Karlis and Ntzoufras 2003] Karlis D, Ntzoufras I. 2003. Analysis of sports data by using bivariate Poisson models. J R Stat Soc Series D. 52: 381-393. doi: 10.1111/1467-9884.00366.
  * [Kass and Raftery 1995] Kass RE, Raftery AE. 1995. Bayes Factors. J Am Stat Assoc. 90(430): 773-795. doi: 10.1080/01621459.1995.10476572.
  * [Knehr et al. 1998] Knehr S, Zitzelsberger H, Bauchinger M. 1998. FISH-based analysis of radiation-induced chromosomal aberrations using different nomenclature systems. Int J Radiat Biol. 73(2):135-41. doi: 10.1080/095530098142509
  * [Kodama et al. 1997] Kodama Y, Nakano M, Ohtaki K, Delongchamp R, Awa AA, Nakamura N. 1997. Estimation of minimal size of translocated chromosome segments detectable by fluorescence in situ hybridization. Int J Radiat Biol. 71(1), 35-39. doi: 10.1080/095530097144391.
  * [Kruschke 2014] Kruschke J. 2014. Doing Bayesian data analysis: A tutorial with R, JAGS, and Stan. 2nd ed. Boston (MA): Academic Press / Elsevier.
  * [Lindholm et al. 1998] Lindholm C, Luomahaara S, Koivistoinen A, Ilus T, Edwards AA, Salomaa S. 1998. Comparison of dose-response curves for chromosomal aberrations established by chromosome painting and conventional analysis. Int J Radiat Biol. 74(1):27-34. doi: 10.1080/095530098141690.
  * [Lucas et al. 1992] Lucas JN, Awa A, Straume T, Poggensee M, Kodama Y, Nakano M, Ohtaki K, Weier HU, Pinkel D, Gray J, et al. 1992. Rapid translocation frequency analysis in humans decades after exposure to ionizing radiation. Int J Radiat Biol. 62(1):53-63. doi: 10.1080/09553009214551821.
  * [Lucas 1996] Lucas JN. Theoretical predictions on the equality of radiation-produced dicentrics and translocations detected by chromosome painting. Int J Radiat Biol. 1996 Feb;69(2):145-53. doi: 10.1080/095530096145977.
  * [Młynarczyk et al. 2022] Młynarczyk D, Puig P, Armero C, Gómez-Rubio V, Barquinero JF, Pujol-Canadell M. 2022. Radiation dose estimation with time-since-exposure uncertainty using the $\gamma$-H2AX biomarker. Sci Rep. 12(1):19877. doi:10.1038/s41598-022-24331-1.
  * [Pujol et al. 2016] Pujol M, Barrios L, Puig P, Caballín MR, Barquinero JF. 2016. A New Model for Biological Dose Assessment in Cases of Heterogeneous Exposures to Ionizing Radiation. Radiat Res. 185(2):151-62. doi: 10.1667/RR14145.1.
  * [Pujol et al. 2020] Pujol-Canadell M, Puig R, Armengol G, Barrios L, Barquinero JF. 2020. Chromosomal aberration dynamics through the cell cycle. DNA Repair (Amst). 89:102838. doi: 10.1016/j.dnarep.2020.102838.
  * [Savage and Papworth 1982] Savage JR, Papworth DG. 1982. Frequency and distribution studies of asymmetrical versus symmetrical chromosome aberrations. Mutat Res. 95(1):7-18. doi: 10.1016/0027-5107(82)90062-8.
  * [Savage and Simpson 1994] Savage JR, Simpson P. 1994. On the scoring of FISH-”painted” chromosome-type exchange aberrations. Mutat Res. 307(1):345-53. doi: 10.1016/0027-5107(94)90308-5.
  * [Schmid et al. 1992] Schmid E, Zitzelsberger H, Braselmann H, Gray JW, Bauchinger M. 1992. Radiation-induced chromosome aberrations analysed by fluorescence in situ hybridization with a triple combination of composite whole chromosome-specific DNA probes. Int J Radiat Biol. 62(6):673-8. doi: 10.1080/09553009214552621.
  * [Tucker et al. 1995] Tucker JD, Morgan WF, Awa AA, Bauchinger M, Blakey D, Cornforth MN, Littlefield LG, Natarajan AT, Shasserre C. 1995. PAINT: a proposed nomenclature for structural aberrations detected by whole chromosome painting. Mutat Res. 347(1):21-4. doi: 10.1016/0165-7992(95)90028-4.