aainstitutetext: Institute for Theoretical Physics, Karlsruhe Institute of
Technology (KIT), 76131 Karlsruhe, Germanybbinstitutetext: Rudolf Peierls
Centre for Theoretical Physics, Parks Road, Oxford OX1 3PU, UK

# SMEFT predictions for $gg\to hh$ at full NLO QCD and truncation
uncertainties

Gudrun Heinrich a Jannis Lang b Ludovic Scyboz<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We present a calculation of the NLO QCD corrections for Higgs-boson pair
production in gluon fusion including effects of anomalous couplings within
Standard Model Effective Field Theory (SMEFT). We study effects of different
truncation options of the EFT expansion in $1/\Lambda$ and of double operator
insertions, both at total cross-section level as well as for the distribution
of the invariant mass of the Higgs-boson pair, at $\sqrt{s}=13$ TeV. The NLO
corrections are implemented in the generator ggHH_SMEFT in the Powheg-Box-V2
framework.

###### Keywords:

LHC, Higgs-boson couplings, NLO, EFT

††preprint: KA-TP-14-2022
. OUTP-22-05P
. P3H-22-045

##  1 Introduction

The importance of Higgs-boson pair production as a process allowing us to shed
more light on the Higgs potential is undisputed. Deviations from the Standard
Model (SM) form, manifesting themselves in anomalous Higgs-boson self-
couplings, would be a clear sign of new physics. If the Higgs-boson trilinear
coupling was found to be different from the SM value, small deviations in
other Higgs couplings could also be expected. In order to expose such slight
deviations, it is crucial to control the uncertainties of the theory
predictions, including the description of anomalous couplings within an
effective field theory (EFT) framework. The theory uncertainties have various
sources, the dominant ones in the SM currently being uncertainties related to
the treatment of the top-quark mass in different renormalisation schemes. NLO
QCD corrections including the full top-quark mass dependence are available
Borowka:2016ehy ; Borowka:2016ypz ; Baglio:2018lrj ; Baglio:2020ini  and have
been included in calculations where higher orders have been performed in the
heavy-top limit Grazzini:2018bsd ; Chen:2019lzz ; Chen:2019fhs , thus reducing
the scale uncertainties and the uncertainties due to missing top-quark mass
effects, while the top-mass scheme uncertainties remain an issue
Baglio:2020wgt .

Going beyond the SM description of the process $gg\to hh$, effects of
anomalous couplings have been studied at NLO in the Born-improved heavy-top
limit (HTL) for both CP-conserving Grober:2015cwa  as well as CP-violating
Grober:2017gut  operators. NLO corrections with full $m_{t}$-dependence have
been incorporated within a non-linear EFT parametrisation (also called Higgs
Effective Field Theory, HEFT) in Ref. Buchalla:2018yce . In Ref.
deFlorian:2017qfk  the combination of NNLO corrections in the HTL has been
performed within the HEFT framework. Finally, in Ref. deFlorian:2021azd , the
full NLO corrections of Ref. Buchalla:2018yce  have been combined with the
NNLO corrections of Ref. deFlorian:2017qfk  to provide approximate NNLO
predictions, dubbed NNLO′, which include the full top-quark mass dependence up
to NLO and higher order corrections up to NNLO in the HTL, combined with
operators related to the five most relevant anomalous couplings for the
process $gg\to hh$.

By including the EFT parametrisation of new physics effects into the
predictions for Higgs-boson pair production, new uncertainties arise, related
to the truncation of the EFT expansion, which relies on an assessment of the
relevance of operator contributions to the Lagrangian in a certain well-
defined counting scheme. Furthermore, at amplitude squared level, there are
several possibilities to truncate the expansion in the canonical dimension,
related to the inclusion of squared dimension six terms and double operator
insertions. The discussion of truncation uncertainties recently gained
considerable attention in the literature Biekotter:2016ecg ; Brivio:2022pyi ;
Dawson:2021xei ; Lang:2021hnd ; Ethier:2021bye ; Trott:2021vqa ;
Martin:2021cvs ; Battaglia:2021nys ; Aoude:2022imd , showing that the
uncertainties have to be assessed on a case-by-case basis.

In the following we will present results within SMEFT for Higgs-boson pair
production in gluon fusion, including full NLO QCD corrections as well as
options to include linear or quadratic terms in the EFT expansion and double
operator insertions. This allows us to investigate various scenarios of
truncation of the EFT expansion and to assess the related uncertainties
compared to the size of the scale uncertainties at NLO QCD. We will also point
out differences between a SMEFT and a HEFT description in this context. In
section 2, we recap the definition of the SMEFT and HEFT effective field-
theory frameworks, specifically of those operators that are of interest for
Higgs-pair production, and we outline different ways of defining the
truncation at the level of the squared amplitude. In section 3 we give a brief
summary of the implementation of the SMEFT predictions at full NLO QCD into
the Powheg-Box-V2 Nason:2004rx ; Frixione:2007vw ; Alioli:2010xd  event
generator. We discuss our results in section 4, and conclude in section 5.

##  2 Description of the EFT framework and the calculation

###  2.1 EFT descriptions of Higgs-boson pair production

In this section, we introduce our conventions and contrast the SMEFT and HEFT
descriptions at Lagrangian level.

In Standard Model Effective Field Theory (SMEFT) Buchmuller:1985jz ;
Grzadkowski:2010es ; Brivio:2017vri , a low energy description of unknown
interactions at a new physics scale $\Lambda$ is constructed as an expansion
in inverse powers of $\Lambda$, with operators $\mathcal{O}_{i}$ of canonical
dimension larger than four and corresponding Wilson coefficients $C_{i}$,

$\mathcal{L}_{\text{SMEFT}}=\mathcal{L}_{\text{SM}}+\sum_{i}\frac{C_{i}^{(6)}}{\Lambda^{2}}\mathcal{O}_{i}^{\rm{dim6}}+{\cal
O}(\frac{1}{\Lambda^{3}})\;.$ (1)

In SMEFT it is assumed that the physical Higgs boson is part of a doublet
transforming linearly under $SU(2)_{L}\times U(1)$. The SMEFT Lagrangian is
usually given in the so-called Warsaw basis Grzadkowski:2010es , where the
terms relevant to the process $gg\to hh$ read

$\begin{split}\Delta\mathcal{L}_{\text{Warsaw}}&=\frac{C_{H,\Box}}{\Lambda^{2}}(\phi^{\dagger}\phi)\Box(\phi^{\dagger}\phi)+\frac{C_{HD}}{\Lambda^{2}}(\phi^{\dagger}D_{\mu}\phi)^{*}(\phi^{\dagger}D^{\mu}\phi)+\frac{C_{H}}{\Lambda^{2}}(\phi^{\dagger}\phi)^{3}\\\
&+\left(\frac{C_{uH}}{\Lambda^{2}}\phi^{\dagger}{\phi}\bar{q}_{L}\phi^{c}t_{R}+h.c.\right)+\frac{C_{HG}}{\Lambda^{2}}\phi^{\dagger}\phi
G_{\mu\nu}^{a}G^{\mu\nu,a}\;.\end{split}$ (2)

The dipole operator $\bar{{\cal O}}_{uG}$ is not included here because it can
be shown that it carries an extra loop suppression factor $1/(4\pi)^{2}$
relative to the other contributions if weak coupling to the heavy sector is
assumed Buchalla:2022vjp ; Buchalla:2018yce ; Arzt:1994gp . In the case of a
UV completion where the coupling to the heavy sector is strong, SMEFT would
not be the appropriate description of the full theory at low energies anyway.

Higgs Effective Field Theory (HEFT) Feruglio:1992wf ; Burgess:1999ha ;
Grinstein:2007iv ; Contino:2010mh ; Alonso:2012px ; Buchalla:2013rka  instead
is based on an expansion in terms of loop orders, which also can be formulated
in terms of chiral dimension ($d_{\chi}$) counting Weinberg:1978kz ;
Buchalla:2013eza ; Krause:2016uhw . The expansion parameter is given by
$f^{2}/\Lambda^{2}\simeq\frac{1}{16\pi^{2}}$, where $f$ is a typical energy
scale at which the EFT expansion is valid (such as the pion decay constant in
chiral perturbation theory),

$\displaystyle{\cal L}_{d_{\chi}}={\cal
L}_{(d_{\chi}=2)}+\sum_{L=1}^{\infty}\sum_{i}\left(\frac{1}{16\pi^{2}}\right)^{L}c_{i}^{(L)}O^{(L)}_{i}\;.$
(3)

The HEFT Lagrangian relevant to Higgs-boson pair production in gluon fusion
can be parametrised by five a priori independent anomalous couplings as
follows Buchalla:2018yce

$\displaystyle\Delta{\cal
L}_{\text{HEFT}}=-m_{t}\left(c_{t}\frac{h}{v}+c_{tt}\frac{h^{2}}{v^{2}}\right)\,\bar{t}\,t-c_{hhh}\frac{m_{h}^{2}}{2v}h^{3}+\frac{\alpha_{s}}{8\pi}\left(c_{ggh}\frac{h}{v}+c_{gghh}\frac{h^{2}}{v^{2}}\right)\,G^{a}_{\mu\nu}G^{a,\mu\nu}\;.$
(4)

Expanding the Higgs doublet in eq. (2) around its vacuum expectation value and
applying a field redefinition to the physical Higgs boson

$\displaystyle h\to
h+v^{2}\frac{C_{H,\textrm{kin}}}{\Lambda^{2}}\left(h+\frac{h^{2}}{v}+\frac{h^{3}}{3v^{2}}\right)\;,$
(5)

with

$C_{H,\textrm{kin}}:=C_{H,\Box}-\frac{1}{4}\,C_{HD}\;,$

the Higgs kinetic term acquires its canonical form (up to ${\cal
O}\left(\Lambda^{-4}\right)$ terms). After that, the couplings can be related
through a comparison of the coefficients of the corresponding terms in the
Lagrangian, which leads to the expressions given in Table 1. Note that in the
Warsaw basis $C_{HG}$ implicitly contains a factor of $\alpha_{s}$ relative to
$c_{ggh}$ and $c_{gghh}$ and therefore the translation becomes scale-dependent
even if no (electroweak) RGE running of the Wilson coefficients is included.

HEFT | Warsaw
---|---
$c_{hhh}$ | $1-2\frac{v^{2}}{\Lambda^{2}}\frac{v^{2}}{m_{h}^{2}}\,C_{H}+3\frac{v^{2}}{\Lambda^{2}}\,C_{H,\textrm{kin}}$
$c_{t}$ | $1+\frac{v^{2}}{\Lambda^{2}}\,C_{H,\textrm{kin}}-\frac{v^{2}}{\Lambda^{2}}\frac{v}{\sqrt{2}m_{t}}\,C_{uH}$
$c_{tt}$ | $-\frac{v^{2}}{\Lambda^{2}}\frac{3v}{2\sqrt{2}m_{t}}\,C_{uH}+\frac{v^{2}}{\Lambda^{2}}\,C_{H,\textrm{kin}}$
$c_{ggh}$ | $\frac{v^{2}}{\Lambda^{2}}\frac{8\pi}{\alpha_{s}}\,C_{HG}$
$c_{gghh}$ | $\frac{v^{2}}{\Lambda^{2}}\frac{4\pi}{\alpha_{s}}\,C_{HG}$
Table 1: Translation at Lagrangian level between different operator basis
choices.

The translation given in Table 1 suggests that there is no explicit dependence
on the scale $\Lambda$ in the HEFT Wilson coefficients. As mentioned above,
the effective HEFT Lagrangian is expanded in powers of
$f^{2}/\Lambda^{2}\simeq 1/(4\pi)^{2}$, with $\Lambda\simeq 4\pi f$ the scale
of new physics and $f$ a reference scale for energies where the EFT expansion
is valid; for the case of strongly coupled UV completions $f$ corresponds to
the scale of dynamical symmetry breaking. In section 4, we will still use the
translation for a specific value of the scale $\Lambda$ to compare SMEFT and
HEFT results.

However, it should be pointed out that a translation between the coefficients
at Lagrangian level must be applied with care. The EFT parametrisations have a
validity range limited by unitarity constraints and the assumption that
$C_{i}/\Lambda^{2}$ in SMEFT is a small quantity. Furthermore, due to
different assumptions about the transformation of the Higgs field under the EW
symmetry transformations, there are relations between certain coefficients in
SMEFT, which are not present in HEFT. Therefore a naive translation from HEFT
(which is, in this regard, the more general theory) to SMEFT can lead out of
the validity range of SMEFT for certain points in the coupling parameter
space, even though they are perfectly valid points in HEFT.

###  2.2 Operator insertions at amplitude-squared level

An EFT description is based on an expansion in a parameter encoding the scale
hierarchies which underlie the EFT description. Therefore an uncertainty
arises due to the truncation of the EFT expansion. In particular, there is the
question whether to truncate the SMEFT expansion at amplitude-squared level
strictly beyond dimension-6, or to include squared dimension-6 operators,
which are of order $1/\Lambda^{4}$ and therefore formally suppressed at the
same order as dimension-8 operators. Furthermore, double operator insertions
into a single diagram are usually neglected as they form a subset of operators
at order $1/\Lambda^{4}$ which is not uniquely defined (for example, they can
be related to a different set of operators through the equations of motion).
Such issues have been discussed recently in Refs. Brivio:2022pyi ;
Dawson:2021xei ; Battaglia:2021nys ; Trott:2021vqa ; Martin:2021cvs .

In the next section we present a Monte Carlo program, which includes the full
$m_{t}$-dependent NLO QCD corrections, and allows the systematic study of
truncation effects for the case of Higgs-boson pair production in gluon
fusion. In order to construct the different truncation options we divide the
amplitude into three parts: the pure SM contribution (SM), single dimension-6
operator insertions ($\rm{dim6}$) and double dimension-6 operator insertions
($\rm{dim6}^{2}$):

$\displaystyle{\cal M}=$ $\displaystyle\
\vbox{\hbox{\includegraphics[page=1,scale={0.9}]{figures/gghh_diagrams}}}+\vbox{\hbox{\includegraphics[page=2,scale={0.9}]{figures/gghh_diagrams}}}+\vbox{\hbox{\includegraphics[page=3,scale={0.9}]{figures/gghh_diagrams}}}$
$\displaystyle\
+\vbox{\hbox{\includegraphics[page=4,scale={0.9}]{figures/gghh_diagrams}}}+\vbox{\hbox{\includegraphics[page=5,scale={0.9}]{figures/gghh_diagrams}}}+\dots$
$\displaystyle=$ $\displaystyle\ {\cal M}_{\text{SM}}+{\cal
M}_{\rm{dim6}}+{\cal M}_{\rm{dim6}^{2}}\;,$ (6)

where $C^{\prime}$ denotes the corresponding coupling combination listed in
Table 1. We consider four possibilities to choose which parts of $|{\cal
M}|^{2}$ from eq. (6) may enter in the squared amplitude forming the cross
section:

$\displaystyle\sigma\simeq\left\\{\begin{aligned} &\
\sigma_{\text{SM}}+\sigma_{\text{SM}\times\rm{dim6}}&\textrm{(a)}\\\ &\
\sigma_{\left(\text{SM}+\rm{dim6}\right)\times\left(\text{SM}+\rm{dim6}\right)}&\textrm{(b)}\\\
&\
\sigma_{\left(\text{SM}+\rm{dim6}\right)\times\left(\text{SM}+\rm{dim6}\right)}+\sigma_{\text{SM}\times\rm{dim6}^{2}}&\textrm{(c)}\\\
&\
\sigma_{\left(\text{SM}+\rm{dim6}+\rm{dim6}^{2}\right)\times\left(\text{SM}+\rm{dim6}+\rm{dim6}^{2}\right)}&\textrm{(d)}\end{aligned}\right.$
(7)

Case (a) denotes the first order of an expansion of $\sigma\sim|{\cal M}|^{2}$
in $\Lambda^{-2}$, (b) is the first order of an expansion of ${\cal M}$ in
$\Lambda^{-2}$. Case (c) includes all terms of ${\cal
O}\left(\Lambda^{-4}\right)$ coming from single and double dimension-6
operator insertions, however it lacks the contribution at the same order from
dimension-8 operators and ${\cal O}\left(\Lambda^{-4}\right)$ terms following
from the field redefinition of eq. (5). Case (d) would correspond to HEFT upon
using the translation of the parameters as in Table 1, except for differences
due to the running in $\alpha_{s}$, because no linearisation whatsoever in
$1/\Lambda$ is present there.

##  3 Implementation and usage of the code within the Powheg-Box

###  3.1 Implementation of the NLO QCD corrections

Our implementation is similar to the NLO HEFT code Powheg-Box-V2/ggHH, which
was presented in Refs. Heinrich:2019bkc ; Heinrich:2020ckp  and is publicly
available. For the virtual two-loop corrections, the NLO QCD corrections in
the SM as calculated in Refs. Borowka:2016ehy ; Borowka:2016ypz  have been
used and extended to the SMEFT framework. The operator insertions have been
included in a modular way, such that the different options described in the
previous section can be calculated.

For the real emission part, the one-loop $2\to 3$ matrix elements have been
produced with GoSam Cullen:2011ac ; Cullen:2014yla , based on a model which we
generated in UFO format Degrande:2011ua , specifying the anomalous couplings
such that GoSam is able to calculate the different contributions according to
the chosen truncation option. The existing interface Luisoni:2013cuh  to
Powheg Nason:2004rx ; Frixione:2007vw ; Alioli:2010xd  has been modified to
hand over event parameters to GoSam in such a way that the factor
$\alpha_{s}$, which is included in the definition of the Higgs-gluon couplings
in SMEFT, is evaluated at the correct energy scale. The renormalisation of the
top-quark mass is performed in the on-shell scheme.

The code is built such that it splits the amplitude evaluation according to
eq. (6), and the squared amplitude is calculated with truncation option (a),
(b), (c) or (d) as defined by the user in the input card.

###  3.2 Usage within the Powheg-Box-V2

Usage of the program ggHH_SMEFT is similar to the HEFT public code ggHH
Heinrich:2020ckp , and both are provided within the POWHEG-BOX-V2
Alioli:2010xd  under User-Processes-V2. The user defines the EFT parameters
directly in the input card, of which we provide an example in
testrun/powheg.input-save. The EFT mass scale is set by:

Lambda=1.0 :

the input value of the SMEFT heavy mass scale $\Lambda$ (in TeV).

The value of the SMEFT coefficients in the Warsaw basis is set by the
following keywords:

CHbox :

the Higgs kinetic term coefficient $C_{H,\Box}$,

CHD :

the Higgs kinetic term coefficient $C_{HD}$,

CH :

the Higgs trilinear coupling term $C_{H}$,

CuH :

the Higgs Yukawa coupling to up-type quarks term $C_{uH}$,

CHG :

the effective coupling of gluons to Higgs bosons $C_{HG}$.

The input keyword usesmeft can be set to 0 or 1. When usesmeft=1, the above
parameters for CHbox, CHD, CH, CuH, CHG are taken from the input card and
translated internally to be used in the computation of the amplitude. If
usesmeft=0, the values of the parameters chhh, ct, ctt, cggh, cgghh are used
instead. If usesmeft=1 is set, the values for the parameters chhh, ct, ctt,
cggh, cgghh will be ignored.

Finally, the different SMEFT truncation options can be selected via the
keyword multiple-insertion, where the options (a)–(d) in eq. (7) correspond to
the values 0–3 of this flag. Apart from the above, the usage of the code is as
described in Heinrich:2020ckp  and in the Docs folder of the code. In
particular, we remind the reader that the 2-loop grids for the virtual
contributions have been generated with a fixed value of the Higgs and the top-
quark masses, set respectively to $m_{h}=125$ GeV and $m_{t}=173$ GeV, and
that these should not be changed (unless the user would like to calculate the
leading order only). For the generation of full-fledged Monte-Carlo events an
interface to Pythia 8 Sjostrand:2004ef ; Sjostrand:2014zea  and Herwig 7
Marchesini:1987cf ; Bellm:2015jjp  is provided that is identical to the one in
Powheg-Box-V2/ggHH.

##  4 Results

Our results were produced for a centre-of-mass energy of $\sqrt{s}=13$ TeV
using the PDF4LHC15_nlo_30_pdfas Butterworth:2015oua  parton distribution
functions interfaced to our code via LHAPDF Buckley:2014ana , along with the
corresponding value for $\alpha_{s}$. The masses of the Higgs boson and the
top quark have been fixed, as in the virtual amplitude, to $m_{h}=125$ GeV,
$m_{t}=173$ GeV and their widths have been set to zero. Jets are clustered
with the anti-$k_{T}$ algorithm Cacciari:2008gp  as implemented in the FastJet
package Cacciari:2005hq ; Cacciari:2011ma , with jet radius $R=0.4$ and a
minimum transverse momentum $p_{T,\mathrm{min}}^{\rm{jet}}=20$ GeV. We set the
renormalisation and factorisation scales to $\mu_{R}=\mu_{F}=m_{hh}/2$.

###  4.1 Total cross sections

We consider three benchmark points, given in Table 2. Each of the benchmark
points belongs to a characteristic $m_{hh}$ shape cluster at NLO, as derived
in Ref. Capozi:2019xsi .111Here, we have updated a subset of the benchmark
points identified in Ref. Capozi:2019xsi , performing the same clustering into
seven characteristic shapes of the $m_{hh}$ distribution with the help of
unsupervised machine learning applied to HEFT results. However, we extracted
central benchmarks under slightly tighter constraints to reflect more recent
experimental measurements CMS:2020gsy ; ATLAS:2021vrm , i.e. we used $0.83\leq
c_{t}\leq 1.17$ for all benchmarks, and additionally $|c_{tt}|<0.05$ in the
case of benchmark 1. The values for $c_{gghh}$ have been modified in order to
fulfil the SMEFT relation $c_{ggh}=2c_{gghh}$. Here we consider three out of
the seven shape types:

  * •

benchmark 1: enhanced low $m_{hh}$ region,

  * •

benchmark 3: enhanced low $m_{hh}$ and second local maximum above
$m_{hh}\simeq 2m_{t}$,

  * •

benchmark 6: close-by double peaks or shoulder left.

| benchmark
---
(${}^{\ast}=$ modified)
$c_{hhh}$ | $c_{t}$ | $c_{tt}$ | $c_{ggh}$ | $c_{gghh}$ | $C_{H,\textrm{kin}}$ | $C_{H}$ | $C_{uH}$ | $C_{HG}$ | $\Lambda$
SM | $1$ | $1$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$ | $0$ |  $1\;$TeV
$1^{\ast}$ | $5.105$ | $1.1$ | $0$ | $0$ | $0$ | $4.95$ | $-6.81$ | $3.28$ | $0$ |  $1\;$TeV
$3^{\ast}$ | $2.21$ | $1.05$ | $-\frac{1}{3}$ | $0.5$ | $0.25^{\ast}$ | $13.5$ | $2.64$ | $12.6$ | $0.0387$ |  $1\;$TeV
$6^{\ast}$ | $-0.684$ | $0.9$ | $-\frac{1}{6}$ | $0.5$ | $0.25$ | $0.561$ | $3.80$ | $2.20$ | $0.0387$ |  $1\;$TeV
Table 2: Benchmark points used for the total cross sections and the
distributions of the invariant mass of the Higgs-boson pair, cf. Table 3 and
Figs. 2–4. The value of $C_{HG}$ is determined using
$\alpha_{s}(m_{Z})=0.118$.

In Table 3 we list total cross sections at 13 TeV for these three benchmark
points for option (b) with $\Lambda=1$ TeV and $\Lambda=2$ TeV, as well as the
ratio to the SM cross section. The result for the total cross section with
option (a) and the result for HEFT are also shown.

benchmark |  $\sigma_{\rm{NLO}}$[fb] | K-factor | ratio to SM |  $\sigma_{\rm{NLO}}$[fb] |  $\sigma_{\rm{NLO}}$[fb]
---|---|---|---|---|---
| option (b) | option (b) | option (b) | option (a) | HEFT
SM | 27.94${}^{+13.7\%}_{-12.8\%}$ | 1.67 | 1 | - | -
$\Lambda=1$ TeV
1 | 74.29${}^{+19.8\%}_{-15.6\%}$ | 2.13 | 2.66 | -61.17 | 94.32
3 | 69.20${}^{+11.7\%}_{-10.3\%}$ | 1.82 | 2.47 | 29.64 | 72.43
6 | 72.51${}^{+20.6\%}_{-16.4\%}$ | 1.90 | 2.60 | 52.89 | 91.40
$\Lambda=2$ TeV
1 | 14.03${}^{+12.0\%}_{-11.9\%}$ | 1.56 | 0.502 | 5.58 | -
3 | 30.81${}^{+16.0\%}_{-14.4\%}$ | 1.71 | 1.10 | 28.35 | -
6 | 35.39${}^{+17.5\%}_{-15.2\%}$ | 1.76 | 1.27 | 34.18 | -
Table 3: Total cross sections for Higgs-boson pair production at full NLO QCD
for three benchmark points and truncation option (b). The total cross sections
for truncation option (a) are also given, in order to highlight the difference
to the linearised case, as well as the values for HEFT. One can clearly see
that truncation option (a) is not a valid option for $\Lambda=1$ TeV, leading
to an unphysical cross section. The uncertainties are scale uncertainties
based on 3-point scale variations.

The squared modulus of the amplitude can be parametrised in terms of all
possible coupling combinations appearing at fixed order, where the $a_{i}$
have been defined as the coefficients of coupling combinations in HEFT at NLO
QCD Buchalla:2018yce ,

$\displaystyle|{\cal M}_{\textrm{BSM}}|^{2}=$ $\displaystyle\,a_{1}\cdot
c_{t}^{4}+a_{2}\cdot c_{tt}^{2}+a_{3}\cdot c_{t}^{2}c_{hhh}^{2}+a_{4}\cdot
c_{ggh}^{2}c_{hhh}^{2}+a_{5}\cdot c_{gghh}^{2}+a_{6}\cdot
c_{tt}c_{t}^{2}+a_{7}\cdot c_{t}^{3}c_{hhh}$ $\displaystyle+$
$\displaystyle\,a_{8}\cdot c_{tt}c_{t}c_{hhh}+a_{9}\cdot
c_{tt}c_{ggh}c_{hhh}+a_{10}\cdot c_{tt}c_{gghh}+a_{11}\cdot
c_{t}^{2}c_{ggh}c_{hhh}+a_{12}\cdot c_{t}^{2}c_{gghh}$ $\displaystyle+$
$\displaystyle\,a_{13}\cdot c_{t}c_{hhh}^{2}c_{ggh}+a_{14}\cdot
c_{t}c_{hhh}c_{gghh}+a_{15}\cdot c_{ggh}c_{hhh}c_{gghh}+a_{16}\cdot
c_{t}^{3}c_{ggh}$ $\displaystyle+$ $\displaystyle\,a_{17}\cdot
c_{t}c_{tt}c_{ggh}+a_{18}\cdot c_{t}c_{ggh}^{2}c_{hhh}+a_{19}\cdot
c_{t}c_{ggh}c_{gghh}+a_{20}\cdot c_{t}^{2}c_{ggh}^{2}$ $\displaystyle+$
$\displaystyle\,a_{21}\cdot c_{tt}c_{ggh}^{2}+a_{22}\cdot
c_{ggh}^{3}c_{hhh}+a_{23}\cdot c_{ggh}^{2}c_{gghh}\;.$

Truncation options (a) and (c) of eq. (7) are expansions at cross-section
level, while (d) is the direct translation from HEFT to SMEFT. Therefore, for
those cases the application of the translation of Table 1, including all terms
at the desired order in $1/\Lambda$, is sufficient. For the truncation option
(b), there are combinations which cannot be reconstructed from HEFT and
therefore they have been calculated explicitly and implemented in analytic
form.

Figure 1:  Heat maps showing the dependence of the cross section on the
couplings $C_{H}$, $C_{uH}$ (left) and $C_{H}$, $C_{H,\textrm{kin}}$ (right)
with $\Lambda=1$ TeV for different truncation options. Top: option (a) (linear
dim-6), middle: option (b) (quadratic dim-6), bottom: option (d) (no
linearisation in $1/\Lambda$). The white areas denote regions in parameter
space where the corresponding cross section would be negative.

In Fig. 1 we show heat maps illustrating the dependence of the cross section
on the couplings $C_{H}$, $C_{uH}$ (left) and $C_{H}$, $C_{H,\textrm{kin}}$
(right) with $\Lambda=1$ TeV for the truncation options (a) (linear dim-6,
upper plots), (b) (quadratic dim-6, middle plots) and (d) (no linearisation,
lower plots). We observe that the results for the total cross sections
(normalised to the SM case) are substantially different between these options.
The purely linear dim-6 contributions lead to negative cross sections over
large parts of the parameter space, as manifested by the white areas in the
top row of plots. This feature is not present at all when we consider the
quadratic dim-6 truncation option. Furthermore, in the linear dim-6 case, we
find completely flat directions in the considered parameter range for a
combined variation of the respective Wilson coefficients. In the quadratic
dim-6 case, option (b), iso-contours have an elliptic shape due to the
quadratic terms in the cross section, while for option (d) the elliptic iso-
contours are distorted due to higher terms in the polynomials in the Wilson
coefficients.

###  4.2 Distribution of the Higgs-pair invariant mass $m_{hh}$: truncation
effects

We now consider differential results and show the effect of the different
truncation options on the distribution of the invariant mass $m_{hh}$ of the
Higgs-boson pair. We present results for the three benchmark points given in
Table 2. We also show scale uncertainty bands for option (b) and the SM case,
resulting from 3-point scale variations around the central scale
$\mu_{R}=\mu_{F}=c\cdot m_{hh}/2$, with $c\in\\{\frac{1}{2},1,2\\}$. For the
SM case as well as for benchmark point 1 we have verified that 7-point scale
variations lie within the 3-point scale uncertainty band.

The SMEFT results for benchmark points 1, 3 and 6 are shown in Figs. 2, 3 and
4, respectively, at $\Lambda=1$ TeV (upper panels), $\Lambda=2$ TeV (middle
panels) and $\Lambda=4$ TeV (lower panels), for the different truncation
options (a) through (d), where the left (right) row shows LO (NLO) results. In
the upper panels, we also display the HEFT results for the same benchmark
point, with the translation given numerically in Table 2 for a value of
$\Lambda=1$ TeV. The dark blue curve corresponds to the linear dim-6 case,
i.e. truncation option (a). As mentioned above, the negative differential
cross-section values in the linear dim-6 case indicate that points in the
coupling parameter space which are valid in HEFT can lead, upon naive
translation of the Wilson coefficients, to parameter points for which the
SMEFT expansion is not valid. The orange curve corresponds to the truncation
option (b), where squared dim-6 contributions are taken into account. This
choice, as well as option (c), leads to positive differential cross
sections.222Nevertheless, one can question whether including an incomplete
subset of $1/\Lambda^{4}$ contributions is theoretically well-defined (in
particular whether such truncation options leave the theory gauge-invariant
Trott:2021vqa ). We take the view that these options can at least be useful in
determining the magnitude of truncation uncertainties in $gg\to hh$.
Furthermore, field redefinition ambiguities are expected to be suppressed in
the $gg\to hh$ case as they enter via interference with a loop-induced SM
amplitude Martin:2021cvs .

Figure 2:  Differential cross sections for the invariant mass $m_{hh}$ of the
Higgs-boson pair for benchmark point 1 of Table 2. Top row: $\Lambda=1$ TeV,
middle row: $\Lambda=2$ TeV, bottom row: $\Lambda=4$ TeV. Left: LO, right:
NLO.

For benchmark point 1, the characteristic shape in HEFT features an enhanced
low-$m_{hh}$ region. From Fig. 2 we see that this shape is not preserved in
SMEFT as $\Lambda$ increases. Disregarding the fully linearised option, which
is simply not a valid option for this benchmark point, we see that the
distribution develops a dip for $\Lambda=2$ TeV; as the other couplings are
almost equal to the SM case, increasing $\Lambda$ translates to decreasing the
effective trilinear coupling, which becomes numerically close to the value of
maximal destructive interference between box- and triangle-contributions. As
the heavy scale is increased further to $\Lambda=4$ TeV, the distribution
approaches the SM shape. Note that for benchmark point 1, where $c_{ggh}$ and
$c_{gghh}$ are zero, curves from option (d) (green) and HEFT (cyan) are
identical at $\Lambda=1$ TeV.

For benchmark point 3, the differences between the truncation options are very
pronounced, see Fig. 3. The value of $c_{hhh}$ is close to the value of
maximal destructive interference between box- and triangle-type contributions
when considering the cross section as a function of $c_{hhh}$ alone, therefore
delicate cancellations are likely to take place. Again, the $m_{hh}$-shape
changes as $\Lambda$ increases from 1 TeV to 2 TeV.

Figure 3:  Differential cross sections for the invariant mass $m_{hh}$ of the
Higgs-boson pair for benchmark point 3 of Table 2. Top row: $\Lambda=1$ TeV,
middle row: $\Lambda=2$ TeV, bottom row: $\Lambda=4$ TeV. Left: LO, right:
NLO.

Furthermore, we observe that the contribution from the interference of double
dimension-6 operator insertions with the SM, which appeared to be sub-dominant
even for $\Lambda=1$ TeV in the case of benchmark point 1, modifies the
interference pattern considerably in the case of benchmark point 3, as can be
seen by comparing truncation option (b) (orange) with option (c) (red), where
the latter includes the double operator insertions interfered with the SM
amplitude. While for $\Lambda=1$ TeV the differences between the truncation
options are large, benchmark point 3 shows a faster convergence to the SM
shape as $\Lambda$ increases. For $\Lambda=4$ TeV, the scale uncertainty bands
largely overlap with the SM uncertainty bands, except at very low and very
high $m_{hh}$.

Figure 4:  Differential cross sections for the invariant mass $m_{hh}$ of the
Higgs-boson pair for benchmark point 6 of Table 2. Top row: $\Lambda=1$ TeV,
middle row: $\Lambda=2$ TeV, bottom row: $\Lambda=4$ TeV. Left: LO, right:
NLO.

For benchmark point 6, the pattern of destructive interference between
different parts of the amplitude (e.g. box- and triangle-type diagrams) in
HEFT is similar to the one in the SM case. However, in SMEFT (taking the
squared dim-6 level – option (b) – as reference), this interference pattern is
modified, leading to a smaller cross section than in HEFT. Furthermore, the
characteristic shape (see Fig. 4) is not preserved for any of the considered
$\Lambda$ values: in HEFT, the characteristic feature of benchmark 6 is a
shoulder left. In SMEFT, this shoulder is absent (except for option (d) and
$\Lambda=1$ TeV, which corresponds to HEFT apart from the running of
$\alpha_{s}$).

Looking at the explicit values of the SMEFT coupling parameters in Table 2,
stemming from the naive translation at $\Lambda=1$ TeV between HEFT and SMEFT,
it becomes clear that the parameters are too large for the SMEFT expansion to
be valid. Therefore the large differences seen in the results at $\Lambda=1$
TeV cannot be regarded as a truncation uncertainty. However, for $\Lambda=2$
TeV these values are divided by a factor of 4 and for $\Lambda=4$ TeV by a
factor of 16, the latter case leading to perfectly valid SMEFT points for all
three benchmarks. Except for benchmark 1, at $\Lambda=4$ TeV, the scale
uncertainty bands overlap with the SM case for almost all the $m_{hh}$ range,
so it seems that the parameter space where the SMEFT expansion is valid but
still clearly distinguishable from the SM case (within NLO uncertainties) is
rather small. This implies that the variety of characteristic $m_{hh}$-shapes
is diminished in SMEFT compared to HEFT.

##  5 Conclusions

We have presented the full NLO QCD corrections to Higgs-boson pair production
in gluon fusion within the SMEFT framework. The corresponding matrix elements
have been implemented in the Powheg-Box-V2333The code can be downloaded from
https://powhegbox.mib.infn.it/ under User-Processes-V2/ggHH_SMEFT. in a
flexible way, and allow the user to investigate the impact of different
truncation options of the series in inverse powers of the new physics scale
$\Lambda$. In particular, we have compared truncations in the inclusion of
dimension-6 operators both at amplitude- as well as at cross-section level,
and we have studied the effects of double operator insertions into single
Feynman diagrams. We have also compared results from non-linear effective
field theory (HEFT) to the different SMEFT truncation options for three
benchmark points. While some of these truncation options only involve a subset
of operators to be included in a consistent counting scheme, truncating the
EFT expansion at the order $1/\Lambda^{2}$ either at cross-section level, or
at amplitude level (where the latter option leads to $1/\Lambda^{4}$-terms at
cross-section level) are both options which are typically employed in fits to
experimental data and which are useful for an estimation of truncation
uncertainties.

It is well known that for values of the Wilson coefficients
$C_{i}\,E^{2}/\Lambda^{2}\ll 1$, the differences between the considered
truncation options should be small for energies typical for LHC processes,
i.e. $E\lesssim 1$ TeV. However, Higgs-boson pair production in gluon fusion
is a process with delicate cancellations between different parts of the
amplitude, such that small differences in the treatment of the Wilson
coefficients can have a rather large effect, especially in differential
distributions.

Furthermore, we have shown that a naive translation between HEFT and SMEFT
Wilson coefficients, based on a comparison of terms at Lagrangian level, has
to be done with great care (if at all). The two counting schemes rely on
different assumptions, and for the SMEFT expansion to be valid, the values of
$C_{i}\,E^{2}/\Lambda^{2}$ should be small. However, starting from a point in
the space of anomalous couplings in HEFT that fulfils all the current
constraints from measurements, after naive translation one easily ends up at
values of $C_{i}$ where the SMEFT expansion in $C_{i}\,E^{2}/\Lambda^{2}$ is
not valid for $E\gtrsim m_{hh}$ and $\Lambda\simeq 1$ TeV, the latter being a
value often used in the literature for this purpose.

We have illustrated the effects of such translations and of the different
truncation options both at total cross-section level as well as for the
$m_{hh}$ distribution, at three benchmark points which are characteristic for
a certain $m_{hh}$-shape in HEFT. In SMEFT, it is clear that the shape must
change for fixed $C_{i}$ as a function of $\Lambda$, as the SM shape is always
approached for large values of $\Lambda$. At $\Lambda=4$ TeV, two of the three
benchmark points lead to a $m_{hh}$ distribution that is so close to the SM
shape that it is indistinguishable from the SM within the NLO scale
uncertainties.

## Acknowledgements

We would like to thank Stephen Jones and Matthias Kerner for collaboration
related to the $ggHH$@NLO project and Gerhard Buchalla and Michael Trott for
useful discussions. This research was supported by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under grant 396021762
- TRR 257. LS is supported by the Royal Society under grant number
RP\R1\180112 and by Somerville College.

## References

  * (1) S. Borowka, N. Greiner, G. Heinrich, S. P. Jones, M. Kerner, J. Schlenk et al., _Higgs Boson Pair Production in Gluon Fusion at Next-to-Leading Order with Full Top-Quark Mass Dependence_ , _Phys. Rev. Lett._ 117 (2016) 012001 [1604.06447].
  * (2) S. Borowka, N. Greiner, G. Heinrich, S. P. Jones, M. Kerner, J. Schlenk et al., _Full top quark mass dependence in Higgs boson pair production at NLO_ , _JHEP_ 10 (2016) 107 [1608.04798].
  * (3) J. Baglio, F. Campanario, S. Glaus, M. Mühlleitner, M. Spira and J. Streicher, _Gluon fusion into Higgs pairs at NLO QCD and the top mass scheme_ , _Eur. Phys. J. C_ 79 (2019) 459 [1811.05692].
  * (4) J. Baglio, F. Campanario, S. Glaus, M. Mühlleitner, J. Ronca, M. Spira et al., _Higgs-Pair Production via Gluon Fusion at Hadron Colliders: NLO QCD Corrections_ , _JHEP_ 04 (2020) 181 [2003.03227].
  * (5) M. Grazzini, G. Heinrich, S. Jones, S. Kallweit, M. Kerner, J. M. Lindert et al., _Higgs boson pair production at NNLO with top quark mass effects_ , _JHEP_ 05 (2018) 059 [1803.02463].
  * (6) L.-B. Chen, H. T. Li, H.-S. Shao and J. Wang, _Higgs boson pair production via gluon fusion at N 3LO in QCD_, _Phys. Lett. B_ 803 (2020) 135292 [1909.06808].
  * (7) L.-B. Chen, H. T. Li, H.-S. Shao and J. Wang, _The gluon-fusion production of Higgs boson pair: N 3LO QCD corrections and top-quark mass effects_, _JHEP_ 03 (2020) 072 [1912.13001].
  * (8) J. Baglio, F. Campanario, S. Glaus, M. Mühlleitner, J. Ronca and M. Spira, _$gg\to HH$ : Combined uncertainties_, _Phys. Rev. D_ 103 (2021) 056002 [2008.11626].
  * (9) R. Gröber, M. Mühlleitner, M. Spira and J. Streicher, _NLO QCD Corrections to Higgs Pair Production including Dimension-6 Operators_ , _JHEP_ 09 (2015) 092 [1504.06577].
  * (10) R. Gröber, M. Mühlleitner and M. Spira, _Higgs Pair Production at NLO QCD for CP-violating Higgs Sectors_ , _Nucl. Phys. B_ 925 (2017) 1 [1705.05314].
  * (11) G. Buchalla, M. Capozi, A. Celis, G. Heinrich and L. Scyboz, _Higgs boson pair production in non-linear Effective Field Theory with full $m_{t}$-dependence at NLO QCD_, _JHEP_ 09 (2018) 057 [1806.05162].
  * (12) D. de Florian, I. Fabre and J. Mazzitelli, _Higgs boson pair production at NNLO in QCD including dimension 6 operators_ , _JHEP_ 10 (2017) 215 [1704.05700].
  * (13) D. de Florian, I. Fabre, G. Heinrich, J. Mazzitelli and L. Scyboz, _Anomalous couplings in Higgs-boson pair production at approximate NNLO QCD_ , _JHEP_ 09 (2021) 161 [2106.14050].
  * (14) A. Biekötter, J. Brehmer and T. Plehn, _Extending the limits of Higgs effective theory_ , _Phys. Rev. D_ 94 (2016) 055032 [1602.05202].
  * (15) I. Brivio et al., _Truncation, validity, uncertainties_ , 2201.04974.
  * (16) S. Dawson, S. Homiller and M. Sullivan, _Impact of dimension-eight SMEFT contributions: A case study_ , _Phys. Rev. D_ 104 (2021) 115013 [2110.06929].
  * (17) J. Lang, S. Liebler, H. Schäfer-Siebert and D. Zeppenfeld, _Effective field theory versus UV-complete model: vector boson scattering as a case study_ , _Eur. Phys. J. C_ 81 (2021) 659 [2103.16517].
  * (18) SMEFiT collaboration, J. J. Ethier, G. Magni, F. Maltoni, L. Mantani, E. R. Nocera, J. Rojo et al., _Combined SMEFT interpretation of Higgs, diboson, and top quark data from the LHC_ , _JHEP_ 11 (2021) 089 [2105.00006].
  * (19) M. Trott, _Methodology for theory uncertainties in the standard model effective field theory_ , _Phys. Rev. D_ 104 (2021) 095023 [2106.13794].
  * (20) A. Martin and M. Trott, _$ggh$ variations_, _Phys. Rev. D_ 105 (2022) 076004 [2109.05595].
  * (21) M. Battaglia, M. Grazzini, M. Spira and M. Wiesemann, _Sensitivity to BSM effects in the Higgs p T spectrum within SMEFT_, _JHEP_ 11 (2021) 173 [2109.02987].
  * (22) R. Aoude, E. Madge, F. Maltoni and L. Mantani, _Quantum SMEFT tomography: top quark pair production at the LHC_ , 2203.05619.
  * (23) P. Nason, _A New method for combining NLO QCD with shower Monte Carlo algorithms_ , _JHEP_ 11 (2004) 040 [hep-ph/0409146].
  * (24) S. Frixione, P. Nason and C. Oleari, _Matching NLO QCD computations with Parton Shower simulations: the POWHEG method_ , _JHEP_ 11 (2007) 070 [0709.2092].
  * (25) S. Alioli, P. Nason, C. Oleari and E. Re, _A general framework for implementing NLO calculations in shower Monte Carlo programs: the POWHEG BOX_ , _JHEP_ 06 (2010) 043 [1002.2581].
  * (26) W. Buchmüller and D. Wyler, _Effective Lagrangian Analysis of New Interactions and Flavor Conservation_ , _Nucl. Phys. B_ 268 (1986) 621.
  * (27) B. Grzadkowski, M. Iskrzynski, M. Misiak and J. Rosiek, _Dimension-Six Terms in the Standard Model Lagrangian_ , _JHEP_ 10 (2010) 085 [1008.4884].
  * (28) I. Brivio and M. Trott, _The Standard Model as an Effective Field Theory_ , _Phys. Rept._ 793 (2019) 1 [1706.08945].
  * (29) G. Buchalla, G. Heinrich, C. Müller-Salditt and F. Pandler, _Loop counting matters in SMEFT_ , 2204.11808.
  * (30) C. Arzt, M. B. Einhorn and J. Wudka, _Patterns of deviation from the standard model_ , _Nucl. Phys. B_ 433 (1995) 41 [hep-ph/9405214].
  * (31) F. Feruglio, _The Chiral approach to the electroweak interactions_ , _Int. J. Mod. Phys. A_ 8 (1993) 4937 [hep-ph/9301281].
  * (32) C. P. Burgess, J. Matias and M. Pospelov, _A Higgs or not a Higgs? What to do if you discover a new scalar particle_ , _Int. J. Mod. Phys. A_ 17 (2002) 1841 [hep-ph/9912459].
  * (33) B. Grinstein and M. Trott, _A Higgs-Higgs bound state due to new physics at a TeV_ , _Phys. Rev. D_ 76 (2007) 073002 [0704.1505].
  * (34) R. Contino, C. Grojean, M. Moretti, F. Piccinini and R. Rattazzi, _Strong Double Higgs Production at the LHC_ , _JHEP_ 05 (2010) 089 [1002.1011].
  * (35) R. Alonso, M. B. Gavela, L. Merlo, S. Rigolin and J. Yepes, _The Effective Chiral Lagrangian for a Light Dynamical ”Higgs Particle”_ , _Phys. Lett. B_ 722 (2013) 330 [1212.3305].
  * (36) G. Buchalla, O. Catà and C. Krause, _Complete Electroweak Chiral Lagrangian with a Light Higgs at NLO_ , _Nucl. Phys. B_ 880 (2014) 552 [1307.5017].
  * (37) S. Weinberg, _Phenomenological Lagrangians_ , _Physica A_ 96 (1979) 327.
  * (38) G. Buchalla, O. Catà and C. Krause, _On the Power Counting in Effective Field Theories_ , _Phys. Lett. B_ 731 (2014) 80 [1312.5624].
  * (39) C. G. Krause, _Higgs Effective Field Theories - Systematics and Applications_ , Ph.D. thesis, Munich U., 2016.  1610.08537.  10.5282/edoc.19873.
  * (40) G. Heinrich, S. P. Jones, M. Kerner, G. Luisoni and L. Scyboz, _Probing the trilinear Higgs boson coupling in di-Higgs production at NLO QCD including parton shower effects_ , _JHEP_ 06 (2019) 066 [1903.08137].
  * (41) G. Heinrich, S. P. Jones, M. Kerner and L. Scyboz, _A non-linear EFT description of $gg\to HH$ at NLO interfaced to POWHEG_, _JHEP_ 10 (2020) 021 [2006.16877].
  * (42) G. Cullen, N. Greiner, G. Heinrich, G. Luisoni, P. Mastrolia, G. Ossola et al., _Automated One-Loop Calculations with GoSam_ , _Eur. Phys. J. C_ 72 (2012) 1889 [1111.2034].
  * (43) G. Cullen et al., _G $\scriptsize{O}$S$\scriptsize{AM}$-2.0: a tool for automated one-loop calculations within the Standard Model and beyond_, _Eur. Phys. J. C_ 74 (2014) 3001 [1404.7096].
  * (44) C. Degrande, C. Duhr, B. Fuks, D. Grellscheid, O. Mattelaer and T. Reiter, _UFO - The Universal FeynRules Output_ , _Comput. Phys. Commun._ 183 (2012) 1201 [1108.2040].
  * (45) G. Luisoni, P. Nason, C. Oleari and F. Tramontano, _$HW^{\pm}$ /HZ + 0 and 1 jet at NLO with the POWHEG BOX interfaced to GoSam and their merging within MiNLO_, _JHEP_ 10 (2013) 083 [1306.2542].
  * (46) T. Sjostrand and P. Z. Skands, _Transverse-momentum-ordered showers and interleaved multiple interactions_ , _Eur. Phys. J._ C39 (2005) 129 [hep-ph/0408302].
  * (47) T. Sjöstrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten et al., _An Introduction to PYTHIA 8.2_ , _Comput. Phys. Commun._ 191 (2015) 159 [1410.3012].
  * (48) G. Marchesini and B. R. Webber, _Monte Carlo Simulation of General Hard Processes with Coherent QCD Radiation_ , _Nucl. Phys._ B310 (1988) 461.
  * (49) J. Bellm et al., _Herwig 7.0/Herwig++ 3.0 release note_ , _Eur. Phys. J. C_ 76 (2016) 196 [1512.01178].
  * (50) J. Butterworth et al., _PDF4LHC recommendations for LHC Run II_ , _J. Phys._ G43 (2016) 023001 [1510.03865].
  * (51) A. Buckley, J. Ferrando, S. Lloyd, K. Nordström, B. Page, M. Rüfenacht et al., _LHAPDF6: parton density access in the LHC precision era_ , _Eur. Phys. J._ C75 (2015) 132 [1412.7420].
  * (52) M. Cacciari, G. P. Salam and G. Soyez, _The anti- $k_{T}$ jet clustering algorithm_, _JHEP_ 04 (2008) 063 [0802.1189].
  * (53) M. Cacciari and G. P. Salam, _Dispelling the $N^{3}$ myth for the $k_{t}$ jet-finder_, _Phys.Lett._ B641 (2006) 57 [hep-ph/0512210].
  * (54) M. Cacciari, G. P. Salam and G. Soyez, _FastJet User Manual_ , _Eur.Phys.J._ C72 (2012) 1896 [1111.6097].
  * (55) M. Capozi and G. Heinrich, _Exploring anomalous couplings in Higgs boson pair production through shape analysis_ , _JHEP_ 03 (2020) 091 [1908.08923].
  * (56) CMS collaboration, _Combined Higgs boson production and decay measurements with up to 137 fb -1 of proton-proton collision data at $\sqrt{s}$ = 13 TeV_, tech. rep., 2020.
  * (57) ATLAS collaboration, _Combined measurements of Higgs boson production and decay using up to $139$ fb-1 of proton-proton collision data at $\sqrt{s}=13$ TeV collected with the ATLAS experiment_, tech. rep., 2021.


Machine learning topological defects in confluent tissues

Andrew Killeen1*, Thibault Bertrand2†, Chiu Fan Lee1‡,

1 Department of Bioengineering, Imperial College London, South Kensington
Campus, London SW7 2AZ, U.K.

2 Department of Mathematics, Imperial College London, South Kensington Campus,
London SW7 2AZ, U.K.

*<EMAIL_ADDRESS>
†<EMAIL_ADDRESS>
‡<EMAIL_ADDRESS>
## Abstract

Active nematics is an emerging paradigm for characterising biological systems.
One aspect of particularly intense focus is the role active nematic defects
play in these systems, as they have been found to mediate a growing number of
biological processes. Accurately detecting and classifying these defects in
biological systems is, therefore, of vital importance to improving our
understanding of such processes. While robust methods for defect detection
exist for systems of elongated constituents, other systems, such as epithelial
layers, are not well suited to such methods. Here, we address this problem by
developing a convolutional neural network to detect and classify nematic
defects in confluent cell layers. Crucially, our method is readily
implementable on experimental images of cell layers and is specifically
designed to be suitable for cells that are not rod-shaped. We demonstrate that
our machine learning model outperforms current defect detection techniques and
that this manifests itself in our method requiring less data to accurately
capture defect properties. This could drastically improve the accuracy of
experimental data interpretation whilst also reducing costs, advancing the
study of nematic defects in biological systems.

## Author summary

Defects in the local alignment of cells have been found to play a functional
role in homeostatic and morphogenetic processes in many different biological
systems. Detecting these defects is, therefore, very important for improving
our understanding of these processes. However, current defect detection
techniques are not well suited to cell layers in which cells are not elongated
in shape, and so the direction of cell orientation can be poorly defined, even
though defects mediate important homeostatic processes in these layers. Here,
we address this problem and develop a machine learning method to detect and
classify defects which is specifically designed for systems for which existing
methods are not appropriate. We show that our method outperforms current
techniques, detecting defects in confluent cell layers more accurately. We
then demonstrate that this improved performance means that properties of these
defects, often the target of experimental studies, can be characterised more
accurately with less data. We anticipate this could drastically improve
experiments investigating defects in these systems, improving our knowledge of
important biological processes.

## Introduction

Tissue dynamics underpins a wide variety of biological processes such as wound
healing [1], cancer metastasis [2] and morphogenesis [3]. Many of these
processes concern confluent tissues, such as epithelial and endothelial cell
layers, making suitable descriptions of the dynamics of these systems a
prerequisite for our understanding of these processes. Unlike constituents in
a passive material, cells within a confluent tissue can generate forces and
exert stresses on their neighbours and underlying substrate. As such, active
matter physics provides a natural framework for describing confluent tissues
and has provided numerous insights into these systems [4]. Active matter is an
emergent field of physics concerned with describing many-body systems far from
equilibrium, where the system is driven from equilibrium by energy expended by
individual constituents [5].

A fruitful connection between active matter and biology is the widely accepted
use of active nematic theory to model the interplay between cell shapes and
tissue dynamics [6, 7]. Nematic systems consist of elongated constituents that
exhibit orientational order with no preferred direction within this
orientation, i.e. they are head-tail symmetric. As such, the orientation of a
cells long axis is a nematic object, and the average local direction of cell
elongation can be thought of as a nematic field. The nematic field is coupled
to the velocity field, with the energy expenditure of individual cells driving
a rich variety of out of equilibrium behavior [8, 9]. A particularly fertile
line of study within confluent tissues is the formation, dynamics and
properties of topological defects in the nematic field [10, 11], as they have
been found to mediate important homeostatic and morphogenetic processes [7].

Topological defects are singularities in the nematic field, points where its
orientation does not vary smoothly but is discontinuous. In active nematic
systems, two types of defects are typically found: comet-shaped singularities,
known as $+1/2$ defects (Fig. 1a), and trefoil-shaped singularities, known as
$-1/2$ defects (Fig. 1b). It is these nematic defects that are being
highlighted as having a functional role in an increasing number of biological
processes. Comet-shaped $+1/2$ defects have been found to trigger cell
extrusion in epithelial layers [12], control the collective dynamics of
confluent layers of neural progenitor cells [13], and have been highlighted as
organisation centres during Hydra morphogenesis [14]. These $+1/2$ defects
also mediate processes in densely packed bacterial systems, triggering the
formation of fruiting bodies in Myxococcus Xanthus colonies [15], as well as
facilitating collective motion in Pseudomonas aeruginosa [16] and E. Coli
colonies [17]. On the other hand, $-1/2$ defects have been associated with
controlling areas of cell depletion in bacterial colonies [15]. We have also
recently shown that active nematic defects can arise in confluent cell layers
with no inherently nematic active forces [18]. Due to the prevalence and
functional role of $\pm 1/2$ in many systems, the efficient detection and
characterisation of topological defects in confluent tissues, and cellular
systems generally, is of fundamental interest to both biology and physics.

Fig 1: Topological defects in confluent tissues. Examples of (a) comet-shaped
$+1/2$ and (b) trefoil-shaped $-1/2$ defects in a confluent cell layer with
the orientation of the long axis of each cell plotted in red.

While defect detection algorithms exist, their application to imaging data
often requires a sophisticated understanding of the underlying physics.
Current algorithms entail locating degenerate points in the nematic field
followed by inspecting how the orientation of the nematic field changes around
this point [19], usually by calculating a quantity known as the winding number
at this point. This can be effective in systems where the nematic field is
well defined across the domain [20, 21] including tissues where cells are
elongated or rod-like, such as spindle shaped fibroblasts [22]. However, this
method is not suited to systems where the nematic field is not well-defined
everywhere, as degenerate points in the field could just be areas of low order
and do not necessarily indicate the existence of a defect. This is often the
case in confluent tissues such as epithelial layers, where the cells are not
rod-like and can be nearly isotropic in shape at times. Previous work studying
defects in these systems searched for defects by calculating the winding
number on a pre-defined grid of points in the nematic field [18, 12, 23, 24].
This method required thousands of defects to be detected to adequately discern
the average properties of defects, such as tissue stress and velocity fields,
which are often the target of experimental studies. The necessity of such
large amounts of data suggests the method of defect detection is inefficient
and imprecise, which begs the question as to whether better methods of
detection are possible for these systems.

One possibility is to utilize machine learning to improve defect detection.
Machine learning methods are being exploited in an increasing variety of
applications in active nematic systems [25, 26, 27]. They have also been used
to detect topological defects in various systems [28, 29, 30], including
cellular systems [31]. This previous study identified degenerate points in the
nematic field of a cellular layer and then used a feed-forward, fully-
connected neural network to perform a binary classification by labelling the
points as either $+1/2$ or $-1/2$ defects. However, as previously discussed,
this method is less applicable to experimental cellular systems where the
nematic field is not well-defined everywhere, and points with low nematic
order do not necessarily indicate the existence of a defect. Additionally,
this study did not demonstrate that its machine learning method outperformed
existing techniques for detection. As such, there is still a need to develop a
machine learning method that can outperform current techniques and be readily
utilized in an experimental setting.

Here, we address this problem by developing a methodology to detect nematic
defects in confluent tissues using a convolutional neural network (CNN). We
design the method such that it is well suited for use in systems that
currently lack effective detection techniques, is user-friendly and readily
implementable on experimental data. In contrast to previous work, we show that
it outperforms current detection techniques and further demonstrate its
efficacy by finding the mean velocity field around $+1/2$ defects and
comparing this to defects detected using the winding number method,
highlighting the improvement in capturing properties of topological defects
with limited data.

## Materials and methods

### Acquisition of test and training data

For our method to be useful, it needs to be suitable for use on experimental
data. For this reason, it takes as its input the $x$ and $y$ coordinates of
each cells centre of mass and the orientation of the long axis of each cell,
both of which can be readily acquired using standard segmentation software
[32]. However, a large amount of data is required to adequately train and test
the model. Here, to obtain sufficient amounts of data, we train and test our
model using data from a numerical model of a confluent cell layer: the active
vertex model (AVM) [33].

The AVM have been used extensively to study epithelial tissue dynamics [34]
and has been found to accurately replicate phenomena observed experimentally
[35, 36]. Moreover, data from an AVM is appropriate for developing our method
as it represents the cell layer in a manner very similar to how they are
represented once experimental images have been segmented (Fig. 2a). We
implement the AVM in the same manner as our previous work [18]. Briefly, we
represent the tissue as a confluent tiling of polygons, the degrees of freedom
being the cell vertices. In the overdamped limit, these vertices move
according to two types of forces: passive mechanical interactions between
cells which arise due to gradients in an effective tissue energy function, and
polar self-propulsive forces that model the motility of each cell. The
effective tissue energy for a tissue containing $N$ cell is

$E=\Sigma^{N}_{a=1}(A_{a}-A_{0})^{2}+(P_{a}-P_{0})^{2}\ ,$ (1)

where $A_{a}$ and $P_{a}$ are the areas and perimeters of cell $a$
respectively, with $A_{0}$ and $P_{0}$ being the target area and perimeter for
each cell. The first term encodes the incompressibility of each cell and the
cell layers resistance to height fluctuations. The second term in the energy
function encodes the competition between cortical tension and cell-cell
adhesion. The force on vertex $i$ due to mechanical interactions is then
$\mathbf{F}_{i}=-\nabla_{i}E$. Self-propulsion is modeled by each cell
generating a polar force of magnitude $f_{0}$, that acts along polarity vector
$\hat{\mathbf{n}}_{a}=(\cos{\theta_{a}},\sin{\theta_{a}})$. The self-
propulsion force on each vertex is then the average self-propulsion of the
three cells that neighbour vertex $i$,
$\mathbf{f}_{i}=\frac{f_{0}}{3}\sum_{a\in\mathcal{N}(i)}\hat{\mathbf{n}}_{a}$,
where $\mathcal{N}(i)$ denotes the list of cells that share vertex $i$. Each
vertex then moves according to

$\frac{d\mathbf{r}_{i}}{dt}=\frac{1}{\zeta}(\mathbf{F}_{i}+\mathbf{f}_{i})\ ,$
(2)

where $\zeta$ is the damping coefficient. The polarity vector of each cell
undergoes rotational diffusion according to

$\frac{d\theta_{a}}{dt}=\sqrt{2D_{r}}\xi_{a}(t)\ ,$ (3)

where $\xi_{a}(t)$ is a white noise process with zero mean and unit variance
and $D_{r}$ is the rotational diffusion coefficient. For a complete
description of the AVM implementation please see [18]. Using the positions of
the cell vertices, we find the centre of mass and long axis orientation of
each cell, which we input into our model. Details of parameter values used can
be found in S1 Appendix.

Fig 2: Defect identification and classification procedure. (a) Example of an
active vertex model configuration. We find the $x$ and $y$ coordinates of each
cells centre of mass, the orientation of each cells long axis is then plotted
at these points. (b) This information is then interpolated to a finer grid to
form the nematic field of the system, where the average local scalar nematic
order parameter $S$ can be calculated at each grid point using a sliding
window. Areas of low order ($S_{th}<0.15$) are identified as possible defect
regions and the centres of mass of these regions are identified (blue dots).
The nematic field around these points is then cropped to form a region of
interest (ROI) (blue box). (c) These ROIs are then input into a machine
learning model which classifies them as a $+1/2$ defect, a $-1/2$ defect or
not a defect.

### Identifying inputs to machine learning model

To process our input such that it is in a form that a machine learning model
can use to classify $\pm 1/2$ defects, we first identify ‘regions of interest’
(ROIs) within the domain. As defects by definition occur in regions with low
nematic order, we identify these areas as ROIs. To do this, we interpolate our
cell orientation data to a fine grid and average each point over a sliding
window to smooth out the data and create a nematic field (Fig. 2b). At each
grid point we then calculate the scalar nematic order parameter $S$, which is
defined as the largest eigenvalue of the nematic tensor $\mathbf{Q}=\langle
2\hat{u}_{m}\hat{u}_{n}-\delta_{mn}\rangle$, where
$\mathbf{\hat{u}}=(\cos{\theta},\sin{\theta})$ with $\theta$ being the
orientation of the field, $\delta_{mn}$ is the kronecker delta and
$\langle\cdot\rangle$ represents a spatial average over a sliding window. $S$
takes a value of one if the local nematic field is perfectly aligned and 0 if
the local field is isotropic. As we seek areas of low order, we identify
contiguous areas in our domain where $S$ is below a threshold value
$S_{th}=0.15$. We then take the centres of mass of these areas to be the
centres of our ROIs, cropping the field in a square around these points (Fig.
2b). We choose the value of $S_{th}$ such that it is high enough to capture
all disordered regions that may contain defects, but low enough such that
these regions are distinct and do not coalesce. We size our ROIs such that
they contain 5-7 cells, large enough so as to capture the core of the defect
but small enough to isolate the defects and avoid capturing multiple defects
in a single ROI. We then use our ROIs as inputs into a machine learning model
that classifies them as containing a $+1/2$ defect, $-1/2$ defect or neither
(Fig. 2c). Details of parameter values used for pre-processing the data can be
found in S1 Appendix. As the position of potential defect locations (the
centres of our ROIs) can be located anywhere in the system domain, our method
is effectively off-lattice in its detection, although they still lie on a fine
grid. This contrasts with previous work detecting defects in epithelial cell
layers [12, 23], which can only detect defects at predefined locations on a
coarse-grained lattice.

### Model architecture and training

We use a CNN to classify our ROIs. A schematic of the architecture can be seen
in Fig. 3c. We use two convolutional layers, each detecting 32 features. Due
to the size of our ROIs, we do not use any max pooling layers after these
convolutional layers. We then follow these convolutional layers with an
additional fully-connected layer of 100 artificial neurons before our output
layer of three neurons, representing our three possible outputs, or classes,
of a $+1/2$ defect, $-1/2$ defect or no defect. Having a third output of no
defect is key here and what makes our method particularly well suited to
epithelial tissues. As the cells in our tissue do not have a well-defined long
axis, neighbouring cells are not always nematically aligned and there can be
regions with low nematic order that do not necessarily contain nematic
defects. Including an option for our CNN to classify an ROI as having no
defect accounts for this possibility. The output of our convolutional and
fully-connected layers are rectified linear units, while the output layer is
softmax [37].

To generate training and testing data we manually classify 5000 ROIs, using
4500 to train our model and saving 500 for testing. To enlarge our training
data, we generate three new copies of each training ROI by rotating each one
by angles $-\pi/2$, $\pi$ and $\pi/2$. Also, as the type of defect is
invariant under reflections, we double this enlarged training data by
reflecting each ROI about its centreline, leading to 36000 training inputs. We
do not enlarge our testing data set.

We train our model by minimising the cross-entropy cost function, defined as

$C=-\sum_{i=1}^{N}\sum^{3}_{c=1}y_{i,c}\log{p_{i,c}}\ ,$ (4)

where $N$ is the number of items in each batch of training data, $y_{i,c}$ is
the correct label (0 or 1) for class $c$ of the $i^{th}$ ROI and $p_{i,c}$ is
the probability calculated by the model that the $i^{th}$ ROI belongs to class
$c$. We minimise $C$ using a stochastic gradient descent algorithm with a
batch size of $N=64$ ROIs. We train over 20 epochs with a learning rate of
0.05 for the first ten epochs and 0.005 for the following ten, initializing
our weights using a Glorot normal distribution [38]. For each epoch, a random
10% of our training images are held back for validation and used at the end of
the epoch to assess the accuracy of our model. Training for more than 20
epochs did not lead to any appreciable improvements in validation accuracy. A
complete list of parameter values used can be found in S1 Appendix.

### Winding number calculation and comparison

Manually labelling the ROIs allows a direct comparison between our method and
the current standard technique used in defect detection, calculation of the
winding number, to be drawn, as we can also classify each ROI by calculating
its winding number. We can then find the accuracy of both methods when
compared against our manually labelled ROIs, our ‘ground truth’. Previous work
on applying machine learning to detect nematic defects in tissues has used the
winding number as the ground truth [31], thereby making it impossible to
determine if the machine learning method is superior to current techniques.

The winding number is the amount the nematic field rotates as a closed loop is
traversed around the centre of the defect [39]. The $\pm 1/2$ defects found in
nematic systems are so called because the nematic field rotates by half a full
rotation, or $\pi$ radians, around the loop (S1 Fig.). The sign of the defects
depends on whether the rotation of the nematic field is in the same direction
as the direction in which the loop is being traversed. If the nematic field
rotates clockwise as the loop is traversed in a clockwise direction, the
defect is positive; if it rotates anti-clockwise, it is negative. We classify
each ROI by finding the winding number on the fine grid around the edge of the
ROI.

## Results

### Machine learning model outperforms winding number classification

The mean performance of our CNN model with each training epoch can be seen in
Fig. 3a. After training, our model clearly outperforms the winding number for
overall classification accuracy on the training data set, defined as the
percentage of correct predictions. However, these are ROIs that our model is
being trained on, meaning it has ‘seen’ them before in previous training
epochs. The real utility of our method depends on its ability to classify ROIs
it has not seen before, which we test using the 500 ROIs in our test data set.
Here our model is again more accurate than the winding number, achieving an
accuracy of 84.0% compared to the winding numbers 76.6%, demonstrating that
our method outperforms the current most widely used technique.

Defect detection techniques can often be sensitive to the window size used to
detect them. If our trained model is to be readily usable on experimental
data, it should achieve accurate results over a range of window sizes. To
investigate this, we assess the accuracy of our trained model and the winding
number method in classifying the test data at different grid sizes (See S1
Appendix). As our model takes as input a $9\times 9$ grid of points, changing
the grid size is akin to changing the ROI size. We find that, over a range of
grid sizes, our method outperforms the winding number method, demonstrating
its robustness in classifying defects even when the ROI size is not well tuned
to the size of defects in the system (See Fig 1 in S1 Appendix).

Fig 3: Machine learning model outperforms winding number classifier. (a) The
loss and validation accuracy of the neural network as it is trained. The
dashed line represents the accuracy of the winding number classification on
the training data (0.812). Error bars represent the standard error in the mean
over 50 realisations. (b) An example domain with defects detected using each
method: our ground truth (GT), neural network (NN) and winding number (W).

As an example of the defects detected using each method, we look at an example
domain from our AVM containing ROIs from our test data set (Fig. 3b). In line
with Fig. 3a, both techniques show good agreement with manually labelled
defects although the winding number appears to detect more false positives
than the neural network. To assess this further and properly delineate the
efficacy of both methods, we break down their performance for each class in
Table 1. We calculate the precision, sensitivity and F1 score of each method.
The precision is defined as

$\textrm{P}=\frac{\textrm{True Positive}}{\textrm{True Positive + False
Positive}}\ ,$ (5)

and determines what proportion of detections are correct. It quantifies how
detrimental false positives are to performance. Sensitivity, on the other
hand, examines the role of false negatives in performance. It is defined as

$\textrm{S}=\frac{\textrm{True Positive}}{\textrm{True Positive + False
Negative}}\ ,$ (6)

and establishes what proportion of true defects are actually identified. The
F1 score, defined as

$\textrm{F1}=\frac{2\textrm{PS}}{\textrm{P + S}}\ ,$ (7)

is a weighted average of the precision and sensitivity and so is a broader
metric of performance.

Table 1:  Performance of defect detection methods on test data. Precision (P), sensitivity (S) and F1 score for each class, as well as an average over all classes, weighted by the size of that class. |  $\mathbf{+1/2}$ | No defect  |  $\mathbf{-1/2}$ | Total
---|---|---|---|---
| P | S | F1  | P | S | F1  | P | S | F1  | P | S | F1
Neural Network  | 0.786 | 0.967 | 0.867  | 0.932 | 0.663 | 0.775  | 0.818 | 0.964 | 0.885  | 0.856 | 0.840 | 0.834
Winding Number  | 0.729 | 0.974 | 0.834  | 0.960 | 0.457 | 0.619  | 0.707 | 1.000 | 0.828  | 0.819 | 0.766 | 0.743

Both methods display a similar pattern of having a higher sensitivity than
precision for both defect categories, but a higher precision than sensitivity
when no defect is present. Additionally, both methods exhibit a lower F1 score
when no defect is present, a reflection of the larger differential between
precision and sensitivity scores. Errors in both methods, therefore, primarily
come from falsely detecting defects, as opposed to missing defects that should
be detected. This information is lost when looking just at the weighted
average values across all classes, which give more comparable precision and
sensitivity scores.

Where the two methods differ is in our CNN model having consistently higher F1
for each category. This is driven by its higher precision in each defect class
and higher sensitivity when no defect is present. The winding number, however,
is slightly more sensitive to detecting defects when they are present. Taken
together, these results show that the improved performance of our neural
network compared to the winding number primarily stems from it detecting fewer
false positive defects. We point out here that it could be argued that
precision and sensitivity should not be weighted equally, as they are in the
F1 score, and that there may be scenarios where ensuring detecting as many
defects as possible is more important than minimising detecting defects that
are not there. However, we now show that, while the winding number may detect
a slightly higher proportion of defects, the higher overall performance of our
model can manifest itself in a wider improvement to experimental results.

### Superior performance leads to improved capturing of defect properties

While results thus far point to the effectiveness of our model, to show that
this realizes itself in tangible improvements to wider results, we look at the
ability of our model, and the winding number, to ascertain the properties of
defects. Experimental studies often seek not only to detect defects but
examine tissue properties around them [12, 23]. To this end, we calculate the
average velocity field around $+1/2$ defects detected using each method. This
observable is particularly pertinent as one can infer global system properties
from the velocity of $+1/2$ defects. The velocity direction indicates whether
the system is behaving as an extensile (the net force on cells is pushing out
along its long axis) or contractile (pulling in along its long axis) nematic,
with tail-to-head motion indicating extensile forces and head-to-tail motion
indicating contractile [11]. Epithelial layers have been shown to exhibit both
forward and backward motion in experiments [23]. It is therefore valuable to
be able to distinguish defect motion accurately and efficiently.

Previous work using this AVM has determined that $+1/2$ defects move in a
tail-to-head direction, indicating extensile behavior, in this system [18].
However, obtaining the characteristic extensile flow field required averaging
the velocity field over many simulations and several thousand defects. While
this was achievable in a numerical model, time and cost constraints could make
the requirements of such vast amounts of data to understand the properties of
these defects prohibitive in an experimental setting. Due to the difficulty in
collecting experimental data, it is therefore crucial that defect properties
can be discerned using a minimal amount of information.

The average velocity fields for manually labelled, winding number detected and
neural network detected defects, using 150 $+1/2$ defects detected from the
test data set, can be seen in Fig. 4. 150 defects were used as this was number
of $+1/2$ defects manually labelled in the test data set and so the largest
number we could use to compare the different methods. The manually labelled
defects demonstrate the clearest tail-to-head, vortical flow fields
characteristic of extensile systems [11]. The flow field around CNN detected
defects clearly show better agreement with the manually labelled flow field
than the flow field found using the winding number, reflected in the higher
correlation between the two fields. The improvement between the two methods is
even more stark when looking at the difference in velocity magnitudes between
each method and the manually labelled flow field (Fig. 5). This illustrates
the impact of the reduced performance, particularly the reduced precision, of
the winding number and confirms the primacy of our model in detecting ‘better’
defects, as the anticipated mean-field behavior is clearer.

Additionally, we compared the velocity field around defects detected using our
model and defects detected on the same data using the ‘on-lattice’ winding
number method used previously [18] (S2 Fig.). The difference between the two
is even starker than between our model and the off-lattice winding number
method used in the present study. As well as highlighting the improvement
using an off-lattice method can bring, it further underlines the benefits of
our model compared to techniques used currently.

Fig 4: Less data is needed to characterise defect properties. Average velocity
fields around $+1/2$ defects for (a) manually labelled defects, (b) defects
detected using the neural network and (c) defects detected using the winding
number. The single point correlation function
($\langle\mathbf{v}_{GT}\cdot\mathbf{v}_{\Box}\rangle$) between the ground
truth (GT) field and the neural network (NN) and winding number (W) fields is
also shown. Fig 5: Winding number velocity field exhibits larger error than
machine learning flow field. Difference in velocity magnitudes between average
velocity fields for manually labelled defects and (a) defects detected using
the neural network and (b) defects detected using the winding number.

## Discussion

In this study, we have developed a new method for detecting nematic defects in
confluent tissues which, crucially, is readily implementable on experimental
data. Our model can therefore aid in the expanding study of characterising
cellular layers as active nematic systems, as active nematic defects are
increasingly found to play functional roles in these systems [6]. Importantly,
we demonstrate that our method displays superior performance to the current
standard use of the winding number in detecting defects and in capturing the
mean-field properties of these defects. This reduces the amount of data
required to obtain these properties, potentially improving experimental data
interpretation.

Interestingly, although the overall performance of our model is better, the
winding number is slightly more sensitive to detecting defects. This means
there could be applications where using the winding number would be more
suitable, if the cost of missing a defect in the domain is very high. However,
the improved performance in finding mean defect flow fields demonstrate that,
in practice, the increase in overall performance of our model makes it more
advantageous. This improved performance is likely due to the winding number
only using information around the edge of the ROI, where as our CNN can take
advantage of spatial information and correlations across the whole region.

In contrast to previous studies on using machine learning to detect nematic
defects [31], our method is specifically designed for noisy experimental
systems where the nematic field may not be well defined everywhere, and
consequently low nematic order may not guarantee the presence of a defect.
However, we anticipate our method will work well with any system whose nematic
field can be easily interpolated to a 2D grid. Indeed, applying our method to
active nematic systems, such as microtubule systems [40], would be an
interesting future application of our technique. Another interesting future
avenue of research includes extending the model to detect integer $+1$
defects, such as spiral or aster shaped singularities, as these have been
engineered to arise in cellular systems [41, 42], and have also been linked to
morphogenetic processes [14, 43].

The trained CNN model along with training data and Python scripts to detect
defects can be found at https://github.com/KilleenA/ML_DefectDetection.

## Supporting information

#### S1 Appendix

Parameter values used in AVM, hyperparameters used in machine learning model
and details of the grid size study.

#### S1 Fig.

Calculating a defects winding number. (a) Example of a trefoil-shaped $-1/2$
defects in a confluent cell layer with the orientation of the long axis of
each cell plotted in red. (b) Characterising this defect by its winding
number. As a closed loop is traversed around the $-1/2$ defect, the
orientation of the cells rotate by $\pi$ radians (half a full rotation), hence
the defect is half-integer. The sign of the defect is negative as the cells
rotate in the opposite direction to the direction of travel around the loop

#### S2 Fig.

Flow fields around $+1/2$ defects. Mean tissue velocity fields around 150
$+1/2$ defects detected using (a) our CNN model and (b) the winding number at
predefined points in the domain, used previously [18].

## Acknowledgments

We acknowledge the High Throughput Computing service provided by Imperial
College Research Computing Service. AK was supported by the EPSRC Centre for
Doctoral Training in Fluid Dynamics Across Scales (Grant EP/L016230/1).

## References

  * 1. Poujade M, Grasland-Mongrain E, Hertzog A, Jouanneau J, Chavrier P, Ladoux B, et al.  Collective migration of an epithelial monolayer in response to a model wound.  Proceedings of the National Academy of Sciences. 2007;104(41):15988–15993.  doi:10.1073/pnas.0705062104.
  * 2. Friedl P, Gilmour D.  Collective cell migration in morphogenesis, regeneration and cancer.  Nature Reviews Molecular Cell Biology. 2009;10(7):445–457.  doi:10.1038/nrm2720.
  * 3. Weijer CJ.  Collective cell migration in development.  Journal of Cell Science. 2009;122(18):3215–3223.  doi:10.1242/jcs.036517.
  * 4. Alert R, Trepat X.  Physical Models of Collective Cell Migration.  Annual Review of Condensed Matter Physics. 2020;11(1):77–101.  doi:10.1146/annurev-conmatphys-031218-013516.
  * 5. Marchetti MC, Joanny JF, Ramaswamy S, Liverpool TB, Prost J, Rao M, et al.  Hydrodynamics of soft active matter.  Reviews of Modern Physics. 2013;85(3):1143–1189.  doi:10.1103/RevModPhys.85.1143.
  * 6. Doostmohammadi A, Ignés-Mullol J, Yeomans JM, Sagués F.  Active nematics.  Nature Communications. 2018;9(1):3246.  doi:10.1038/s41467-018-05666-8.
  * 7. Balasubramaniam L, Mège RM, Ladoux B.  Active nematics across scales from cytoskeleton organization to tissue morphogenesis.  Current Opinion in Genetics & Development. 2022;73:101897.  doi:https://doi.org/10.1016/j.gde.2021.101897.
  * 8. Dombrowski C, Cisneros L, Chatkaew S, Goldstein RE, Kessler JO.  Self-concentration and large-scale coherence in bacterial dynamics.  Physical Review Letters. 2004;93(9):098103.  doi:10.1103/PhysRevLett.93.098103.
  * 9. Opathalage A, Norton MM, Juniper MPN, Langeslay B, Aghvami SA, Fraden S, et al.  Self-organized dynamics and the transition to turbulence of confined active nematics.  Proceedings of the National Academy of Sciences of the United States of America. 2019;116(11):4788–4797.  doi:10.1073/pnas.1816733116.
  * 10. Giomi L, Bowick MJ, Ma X, Marchetti MC.  Defect annihilation and proliferation in active Nematics.  Physical Review Letters. 2013;110(22):228101.  doi:10.1103/PhysRevLett.110.228101.
  * 11. Giomi L, Bowick MJ, Mishra P, Sknepnek R, Cristina Marchetti M.  Defect dynamics in active nematics.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences. 2014;372(2029):20130365.  doi:10.1098/rsta.2013.0365.
  * 12. Saw TB, Doostmohammadi A, Nier V, Kocgozlu L, Thampi S, Toyama Y, et al.  Topological defects in epithelia govern cell death and extrusion.  Nature. 2017;544:212–216.  doi:10.1038/nature21718.
  * 13. Kawaguchi K, Kageyama R, Sano M.  Topological defects control collective dynamics in neural progenitor cell cultures.  Nature. 2017;545(7654):327–331.  doi:10.1038/nature22321.
  * 14. Maroudas-Sacks Y, Garion L, Shani-Zerbib L, Livshits A, Braun E, Keren K.  Topological defects in the nematic order of actin fibres as organization centres of Hydra morphogenesis.  Nature Physics. 2021;17(2):251–259.  doi:10.1038/s41567-020-01083-1.
  * 15. Copenhagen K, Alert R, Wingreen NS, Shaevitz JW.  Topological defects promote layer formation in Myxococcus xanthus colonies.  Nature Physics. 2021;17(2):211–215.  doi:10.1038/s41567-020-01056-4.
  * 16. Meacock OJ, Doostmohammadi A, Foster KR, Yeomans JM, Durham WM.  Bacteria solve the problem of crowding by moving slowly.  Nature Physics. 2021;17(2):205–210.  doi:10.1038/s41567-020-01070-6.
  * 17. Doostmohammadi A, Thampi SP, Yeomans JM.  Defect-Mediated Morphologies in Growing Cell Colonies.  Physical Review Letters. 2016;117(4):048102.  doi:10.1103/PhysRevLett.117.048102.
  * 18. Killeen A, Bertrand T, Lee CF.  Polar Fluctuations Lead to Extensile Nematic Behavior in Confluent Tissues.  Phys Rev Lett. 2022;128:078001.  doi:10.1103/PhysRevLett.128.078001.
  * 19. de Gennes PG.  The physics of liquid crystals.  Oxford: Oxford University Press; 1974.
  * 20. DeCamp SJ, Redner GS, Baskaran A, Hagan MF, Dogic Z.  Orientational order of motile defects in active nematics.  Nature Materials. 2015;14(11):1110–1115.
  * 21. Ellis PW, Pearce DJG, Chang YW, Goldsztein G, Giomi L, Fernandez-Nieves A.  Curvature-induced defect unbinding and dynamics in active nematic toroids.  Nature Physics. 2018;14(1):85–90.  doi:10.1038/nphys4276.
  * 22. Duclos G, Erlenkämper C, Joanny JF, Silberzan P.  Topological defects in confined populations of spindle-shaped cells.  Nature Physics. 2017;13(1):58–62.  doi:10.1038/nphys3876.
  * 23. Balasubramaniam L, Doostmohammadi A, Saw TB, Narayana GHNS, Mueller R, Dang T, et al.  Investigating the nature of active forces in tissues reveals how contractile cells can form extensile monolayers.  Nature Materials. 2021;20:1156–1166.  doi:10.1038/s41563-021-00919-2.
  * 24. Armengol-Collado JM, Carenza LN, Eckert J, Krommydas D, Giomi L. Epithelia are multiscale active liquid crystals; 2022.  Available from: https://arxiv.org/abs/2202.00668.
  * 25. Hedlund E, Hedlund K, Green A, Chowdhury R, Park CS, Maclennan JE, et al.  Detection of islands and droplets on smectic films using machine learning.  Physics of Fluids. 2022;34(10):103608.  doi:10.1063/5.0117358.
  * 26. Colen J, Han M, Zhang R, Redford SA, Lemma LM, Morgan L.  Machine learning active-nematic hydrodynamics.  PNAS. 2021;118(10):e2016708118.  doi:10.1073/pnas.2016708118/-/DCSupplemental.y.
  * 27. Zhou Z, Joshi C, Liu R, Norton MM, Lemma L, Dogic Z, et al.  Machine learning forecasting of active nematics.  Soft Matter. 2021;17(3):738–747.  doi:10.1039/D0SM01316A.
  * 28. Beach MJS, Golubeva A, Melko RG.  Machine learning vortices at the Kosterlitz-Thouless transition.  Physical Review B. 2018;97(4):045207.  doi:10.1103/PhysRevB.97.045207.
  * 29. Walters M, Wei Q, Chen JZY.  Machine learning topological defects of confined liquid crystals in two dimensions.  Physical Review E. 2019;99(6):062701.  doi:10.1103/PhysRevE.99.062701.
  * 30. Minor EN, Howard SD, Green AAS, Glaser MA, Park CS, Clark NA.  End-to-end machine learning for experimental physics: Using simulated data to train a neural network for object detection in video microscopy.  Soft Matter. 2020;16(7):1751–1759.  doi:10.1039/c9sm01979k.
  * 31. Wenzel D, Nestler M, Reuther S, Simon M, Voigt A.  Defects in Active Nematics-Algorithms for Identification and Tracking.  Computational Methods in Applied Mathematics. 2021;21(3):683–692.  doi:10.1515/cmam-2020-0021.
  * 32. Schindelin J, Arganda-Carreras I, Frise E, Kaynig V, Longair M, Pietzsch T, et al.  Fiji: An open-source platform for biological-image analysis.  Nature Methods. 2012;9(7):676–682.  doi:10.1038/nmeth.2019.
  * 33. Sussman DM.  cellGPU: Massively parallel simulations of dynamic vertex models.  Computer Physics Communications. 2017;219:400–406.  doi:https://doi.org/10.1016/j.cpc.2017.06.001.
  * 34. Shaebani MR, Wysocki A, Winkler RG, Gompper G, Rieger H.  Computational models for active matter.  Nature Reviews Physics. 2020;2(4):181–199.  doi:10.1038/s42254-020-0152-1.
  * 35. Petrolli V, Le Goff M, Tadrous M, Martens K, Allier C, Mandula O, et al.  Confinement-Induced Transition between Wavelike Collective Cell Migration Modes.  Physical Review Letters. 2019;122(16):168101.  doi:10.1103/PhysRevLett.122.168101.
  * 36. Henkes S, Kostanjevec K, Collinson JM, Sknepnek R, Bertin E.  Dense active matter model of motion patterns in confluent cell monolayers.  Nature Communications. 2020;11(1):1405.  doi:10.1038/s41467-020-15164-5.
  * 37. Nielsen M.  Neural Networks and Deep Learning.  Determination Press; 2015.
  * 38. Glorot X, Bengio Y.  Understanding the difficulty of training deep feedforward neural networks.  In: Teh YW, Titterington M, editors. Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics. vol. 9 of Proceedings of Machine Learning Research. Chia Laguna Resort, Sardinia, Italy; 2010. p. 249–256.
  * 39. Huterer D, Vachaspati T.  Distribution of singularities in the cosmic microwave background polarization.  Physical Review D - Particles, Fields, Gravitation and Cosmology. 2005;72(4):043004.  doi:10.1103/PhysRevD.72.043004.
  * 40. Sanchez T, Chen DTN, Decamp SJ, Heymann M, Dogic Z.  Spontaneous motion in hierarchically assembled active matter.  Nature. 2012;491(7424):431–434.  doi:10.1038/nature11591.
  * 41. Turiv T, Krieger J, Babakhanova G, Yu H, Shiyanovskii SV, Wei QH, et al.  Topology control of human fibroblast cells monolayer by liquid crystal elastomer.  Science Advances. 2020;6(20):eaaz6485.  doi:10.1126/sciadv.aaz6485.
  * 42. Endresen KD, Kim MS, Pittman M, Chen Y, Serra F.  Topological defects of integer charge in cell monolayers.  Soft Matter. 2021;17(24):5878–5887.  doi:10.1039/D1SM00100K.
  * 43. Guillamat P, Blanch-Mercader C, Pernollet G, Kruse K, Roux A.  Integer topological defects organize stresses driving tissue morphogenesis.  Nature Materials. 2022;21(5):588–597.  doi:10.1038/s41563-022-01194-5.