# Type-III see-saw: Search for triplet fermions in final states with multiple
leptons and fat-jets at 13 TeV LHC

Saiyad Ashanujjaman ID<EMAIL_ADDRESS>Kirtiman Ghosh<EMAIL_ADDRESS>Institute of Physics, Bhubaneswar, Sachivalaya Marg, Sainik School,
Bhubaneswar 751005, India  Homi Bhabha National Institute, Training School
Complex, Anushakti Nagar, Mumbai 400094, India

###### Abstract

The type-III see-saw model holding out a riveting rationale for the minuscule
neutrino masses caters for a well-to-do phenomenology at the Large Hadron
Collider (LHC). Several searches targetting the triplet fermions have been
performed at the LHC. Not only are the signals for the leptonic final states
considered in these searches suppressed by the branching fractions of the
Standard Model (SM) bosons, but they are also beset with considerably large SM
backgrounds. Thus, these searches are deemed not to be sensitive enough in
probing the triplet fermions much heavier than 1 TeV. To this end, we perform
a search for the triplet fermions in final states with multiple leptons and
fat-jets that are cleaner than the usual LHC searches and allow kinematic
reconstruction of the triplets. After performing a systematic and
comprehensive analysis with seven distinct final states, we project the
required luminosities for both $3\sigma$ and $5\sigma$ discoveries of the
triplet fermions as a function of their mass. The triplet fermions with mass
as large as 1265 (1380) and 1480 (1600) GeV could be discovered with $5\sigma$
($3\sigma$) significance at 500 and 3000 fb-1 luminosities, respectively.

Type-III See-saw, Triplet Fermions, Multi-lepton final states, Fat-jet
Signatures.

††preprint: IP/BBSR/2021-08

##  I Introduction

Despite being surpassingly successful, the Standard Model (SM) falls short of
proffering a natural and well-founded rationale for the observed sub-eV
neutrino masses. Diversely, the type-III see-saw model, one of the three UV
completions of the so-called Weinberg operator Weinberg:1979sa  at the tree-
level Ma:1998dn , based on the annexation of the SM by weak gauge triplet of
fermion field Foot:1988aq  seems to explain the minuteness of neutrino masses
readily. Though the see-saw anchors are naturally motivated to have very high
scale masses, if balanced with appropriate Yukawa couplings, nothing precludes
them from having mass at TeV scale, and thus a well-to-do phenomenology at the
Large Hadron Collider (LHC). Phenomenological outcome of this model has been
studied all-encompassingly in the literature Bajc:2006ia ;
FileviezPerez:2007bcw ; FileviezPerez:2007yji ; delAguila:2008hw ;
delAguila:2008cj ; Franceschini:2008pz ; Arhrib:2009mz ; Bandyopadhyay:2011aa
; Biggio:2011ja ; Goswami:2017jqs ; Cai:2017mow ; Das:2020gnt ; Jana:2020qzn ;
Das:2020uer ; Ashanujjaman:2021jhi ; Sen:2021fha . At the LHC, a number of
collider searches looking for the triplet fermions have been carried out
ATLAS:2015fma ; CMS:2017ybg ; ATLAS:2018ghc ; CMS:2019lwf ; ATLAS:2020wop ,
and the observations being consistent with the SM expectations, these searches
have derived exclusion limits at 95% confidence level (CL) on the mass of the
triplet fermions. The most stringent limit of 880 GeV has been put by the
recent CMS multi-lepton search CMS:2019lwf  assuming only one generation of
them in the flavour-democratic scenario. This limit extends to 1110 GeV when
three mass-degenerate triplets111Three generations of triplet fermions are
required to generate three non-vanishing neutrino masses. Note that one (two)
generation(s) of them would not suffice to explain more than one (two) non-
zero neutrino masses. are assumed Ashanujjaman:2021jhi .

Not only are the signals for the leptonic final states considered in the
abovementioned LHC searches ATLAS:2015fma ; CMS:2017ybg ; ATLAS:2018ghc ;
CMS:2019lwf ; ATLAS:2020wop  suppressed by the leptonic branching fractions of
the SM bosons ($W,Z,h$), but they are also beset with considerably large SM
backgrounds. Consequently, these searches are deemed not to be sensitive
enough in probing the triplet fermions much heavier than 1 TeV. Moreover,
these searches are not designed to reconstruct the triplet mass. This behoves
us to look for complementary final states, such as multiple leptons and fat-
jets resulting from the hadronically decaying SM bosons, which are not only
clean and thus sensitive enough in probing the TeV scale triplets but also
allow kinematic reconstruction of them. To this end, we perform a systematic
and comprehensive search for the triplet fermions in final states with
multiple leptons and fat-jets.

The rest of this work is structured as follows. In Section II, we briefly
discuss the productions and decays of the triplet fermions. We perform a
systematic and comprehensive collider analysis for the triplet fermions in
final states with multiple leptons and fat-jets at the 13 TeV LHC in Section
III followed by a summary in Section IV.

##  II  The triplet fermions

The type-III see-saw model employs, in addition to the SM field content, three
generations of right-handed $SU(2)_{L}$ triplet fermions with vanishing
hypercharge

$\Sigma=\left(\begin{array}[]{cc}\Sigma^{0}_{R}/\sqrt{2}&\Sigma^{+}_{R}\\\
\Sigma^{-}_{R}&-\Sigma^{0}_{R}/\sqrt{2}\end{array}\right);$

and the relevant triplet Yukawa and mass terms are given by

$-{\cal L}\supset\frac{1}{2}\left(M_{\Sigma}\right)_{ij}{\rm
Tr}\left(~{}\overline{\tilde{\Sigma}_{i}}\Sigma_{j}\right)+\sqrt{2}~{}\left(Y_{\Sigma}\right)_{i\alpha}~{}\tilde{H}^{\dagger}\overline{\Sigma_{i}}L_{\alpha}~{}+\rm{h.c.},$

where $i,j$ and $\alpha$ are the generation indices which run over 1,2 and 3,
$M_{\Sigma}$ is the Majorana mass matrix for $\Sigma$, $\tilde{\Sigma}$
denotes charge-conjugation of $\Sigma$, i.e.
$\tilde{\Sigma}=C\overline{\Sigma}^{T}$ with $C$ being the charge-conjugation
matrix, $Y_{\Sigma}$ is the Yukawa matrix,
$L=\left(\nu_{L},\ell_{L}\right)^{T}$ is the left-handed SM lepton doublet,
and $H$ is the SM Higgs doublet. For simplicity, we assume
$M_{\Sigma}=m_{\Sigma}\mathbb{1}_{3\times 3}$, i.e. all three generations of
triplets are mass-degenerate. The above Lagrangian leads to the so-called see-
saw formula for the light neutrinos:
$m_{\nu}\approx-\frac{v^{2}}{2}Y_{\Sigma}^{T}M_{\Sigma}^{-1}Y_{\Sigma}$. Fo a
given $M_{\Sigma}$, the well-known Casas-Ibarra parametrisation Casas:2001sr ;
Ibarra:2003up  enables one to find the most general texture of $Y_{\Sigma}$ to
reproduce the low-energy neutrino oscillation parameters Esteban:2018azc
which we avert to discuss in this work, see Jana:2020qzn ; Das:2020uer ;
Ashanujjaman:2021jhi .

The triplet fermions are pair produced copiously at the LHC by quark-antiquark
annihilation via $s$-channel $\gamma/Z$ and $W^{\textpm}$ exchanges

$q\bar{q^{\prime}}\to W^{*}\to\Sigma^{\pm}\Sigma^{0},\quad{\rm and}\quad
q\bar{q}\to\gamma^{*}/Z^{*}\to\Sigma^{+}\Sigma^{-}.$

They are also pair produced via photon-photon fusion processes:222They are
also pair produced via vector-boson fusion processes with one or two
associated forward jets at the LHC. However, such processes are rather sub-
dominant and can be neglected. $\gamma\gamma\to\Sigma^{+}\Sigma^{-}.$ Their
production at the LHC has been widely studied in the literature Bajc:2006ia ;
FileviezPerez:2007bcw ; FileviezPerez:2007yji ; delAguila:2008hw ;
delAguila:2008cj ; Franceschini:2008pz ; Arhrib:2009mz ; Bandyopadhyay:2011aa
; Biggio:2011ja ; Goswami:2017jqs ; Cai:2017mow ; Das:2020gnt ; Jana:2020qzn ;
Das:2020uer ; Ashanujjaman:2021jhi ; Sen:2021fha . In our analysis, we
evaluate the leading order (LO) production cross-sections using the UFO
modules generated from SARAH Staub:2013tta ; Staub:2015kfa  in MadGraph
Alwall:2011uj ; Alwall:2014hca  with the LUXqed17-plus-PDF4LHC15-nnlo-100
parton distribution function Manohar:2016nzj ; Manohar:2017eqh ;
Butterworth:2015oua . Following Ref. Ruiz:2015zca , we naively scale the LO
cross-section by an overall $K$-factor of 1.30.333Ref. Ruiz:2015zca  has
reported the leading QCD corrections to the production of triplet fermions at
hadron colliders. The resulting next-to-leading order (NLO) $K$-factors varies
from 1.24 to 1.37 over a triplet mass range of 1 TeV to 2 TeV. The total pair
production cross-sections varies from 585 fb to 0.02 fb for 400 GeV to 2 TeV
triplet mass.

After being produced, the triplet fermions undergo two types of decays — $(i)$
the heavy state transitions due to the radiative mass-splitting of $\Delta
m=m_{\Sigma^{+}}-m_{\Sigma^{0}}\sim 167$ MeV Cirelli:2005uq  , and $(ii)$ the
two-body decays into a SM boson ($W,Z,h$) and a lepton ($\ell,\nu$). For the
Yukawa couplings $\left(Y_{\Sigma}\right)_{i\alpha}$ not smaller than
$\mathcal{O}(10^{-8})$, the later decay modes subjugate the former ones
ensuring their prompt decays into two SM particles Jana:2020qzn ; Das:2020uer
; Ashanujjaman:2021jhi . The relevant partial decay widths are given by

$\Gamma(\Sigma^{0(\pm)}_{i}\to
XY^{(\prime)}_{\alpha})=\kappa\frac{m_{\Sigma}}{64\pi}|(Y_{\Sigma})_{i\alpha}|^{2}f\left(\frac{m_{X}}{m_{\Sigma}}\right),$

where $XY_{\alpha}\ni h\nu_{\alpha},Z\nu_{\alpha},W^{\pm}\ell^{\mp}_{\alpha}$,
$XY^{\prime}_{\alpha}\ni
h\ell^{\pm}_{\alpha},Z\ell^{\pm}_{\alpha},W^{\pm}\nu_{\alpha}$,
$f(r)=(1-r^{2})^{2}(1+\kappa^{\prime}r^{2}),\kappa(\kappa^{\prime})=1,1,2(0,2,2)$
for $X=h,Z,W^{\pm}$, respectively. The branching fractions for the $h,Z$ and
$W$ decay modes are, respectively, 25%, 25% and 50%. For simplicity, all three
generations of triplets are assumed to be mass-degenerate.444While, a priori,
nothing precludes the mass from varying across the generations, including such
a possibility only introduces additional parameters (and thus complications)
without adding anything qualitatively different to the discussion. Further, we
assume identical branching fractions across all lepton flavours (flavour-
democratic scenario).555We accentuate that the flavour democratic scenario
does not necessarily be an ad-hoc assumption; rather, it could be a
consequence of degenerate mass spectrum for both the light and heavy neutrinos
or could be a repercussion of the information lost in integrating the heavy
leptons out Ashanujjaman:2021jhi . Also, this assumption facilitates an easy
comparison with existing search paradigms adopted by the ATLAS/CMS
collaborations. However, relaxation of this assumption may significantly alter
the phenomenology, and thus the exclusion limit, see Ashanujjaman:2021jhi .
This is readily achieved by choosing the triplet Yukawa couplings that are
compatible with the existing constraints Abada:2008ea ; Goswami:2018jar ;
Biggio:2019eeo  as well as the neutrino oscillation data Esteban:2018azc , see
Jana:2020qzn ; Das:2020uer ; Ashanujjaman:2021jhi .

##  III Collider phenomenology

At the LHC, the triplet fermions are pair produced aplenty followed by their
prompt decays to a SM lepton ($\ell,\nu$) and a boson ($W,Z,h$) leading to
multifarious final state signatures. The characteristic signature of TeV scale
triplets666The recent CMS multi-lepton search CMS:2019lwf  has already
excluded triplet fermions with masses below 880 GeV assuming only one
generation of them in the flavour-democratic scenario. This limit extends to
1110 GeV when three mass-degenerate triplets are considered
Ashanujjaman:2021jhi . would be vastly boosted final sate lepton and boson
which are likely to be largely separated in the azimuth plane. The jets
stemming from the sufficiently highly boosted boson would be collimated. The
hadronically decaying bosons, thus, are more likely to manifest as a single
fat-jet rather than two resolved jets. In what follows, we perform a
systematic and comprehensive collider study for the triplet fermions in the
final states with multiple leptons and fat-jets. The processes contributing to
the final state signatures of our interest are summarised in the following:

$(i)$ One lepton and fat-jets: Both pair and associated production of charged
triplets contribute to this final state:

  1. (a)

$\Sigma^{\pm}\Sigma^{\mp}\to\ell^{\pm}Z/h\nu W^{\mp}\to\ell^{\pm}\nu
J_{Z/h}J_{W}$,

  2. (b)

$\Sigma^{+}\Sigma^{-}\to\nu W^{+}\nu W^{-}\to\ell^{\pm}\nu\nu\nu J_{W}$,

  3. (c)

$\Sigma^{\pm}\Sigma^{0}\to\ell^{\pm}Z/h\nu Z/h\to\ell^{\pm}\nu
J_{Z/h}J_{Z/h}$,

  4. (d)

$\Sigma^{\pm}\Sigma^{0}\to\nu W^{\pm}\ell^{\pm}W^{\mp}\to\ell^{\pm}\nu
J_{W}J_{W}$,

where $J_{W/Z/h}$ denotes a fat-jet emanating from $W/Z/h$. For process (b),
neutrinos result from both triplet fermions making kinematic reconstruction of
the latter challenging. For all other processes, one of the triplets decays to
$\ell^{\pm}$ and $W/Z/h$, where the latter manifests as $J_{W/Z/h}$. The
triplets being in the TeV regime, $\ell^{\pm}$ and $J_{W/Z/h}$ are deemed to
be largely separated in the azimuth plane. Therefore, $m_{\Sigma}$ can be
reconstructed from the invariant mass distribution of the lepton and the
farthest $J_{ZW//h}$ in the azimuth plane.

$(ii)$ Same-sign dilepton and fat-jets: Only the following process contributes
to this final state:

$\Sigma^{\pm}\Sigma^{0}\to\ell^{\pm}Z/h\ell^{\pm}W^{\mp}\to\ell^{\pm}\ell^{\pm}J_{Z/h}J_{W}.$

The invariant mass distribution of $J_{W/Z/h}$ and the farthest lepton in the
azimuth plane would peak at $m_{\Sigma}$.

$(iii)$ Opposite-sign dilepton and fat-jets: The processes contributing to
this final states are as follows:

  1. (a)

$\Sigma^{+}\Sigma^{-}\to\ell^{+}Z/h\ell^{-}Z/h\to\ell^{+}\ell^{-}J_{Z/h}J_{Z/h}$,

  2. (b)

$\Sigma^{\pm}\Sigma^{\mp}\to\ell^{\pm}Z/h\nu
W^{\mp}\to\ell^{\pm}\ell^{\mp}\nu\nu J_{Z/h}$,

  3. (c)

$\Sigma^{\pm}\Sigma^{0}\to\ell^{\pm}Z/h\ell^{\mp}W^{\pm}\to\ell^{\pm}\ell^{\mp}J_{Z/h}J_{W}$.

Similar to the above channel, the invariant mass distribution of $J_{W/Z/h}$
and the farthest lepton in the azimuth plane would yield $m_{\Sigma}$.

$(iv)$ Trilepton without fat-jet: Both pair and associated production of
charged triplets contribute to this final state:

  1. (a)

$\Sigma^{\pm}\Sigma^{0}\to\nu
W^{\pm}\ell^{\pm(\mp)}W^{\mp(\pm)}\to\ell^{\pm}\ell^{\mp}\ell^{\pm}\nu\nu$,

  2. (b)

$\Sigma^{\pm}\Sigma^{0}\to\nu W^{\pm}\nu
Z/h\to\ell^{\pm}\ell^{+}\ell^{-}\nu\nu$,

Here, both the triplets lead to a neutrino in the final state, making their
kinematic reconstruction challenging.

$(v)$ Trilepton and one fat-jet: Both pair and associated production of
charged triplets contribute to this final state:

  1. (a)

$\Sigma^{\pm}\Sigma^{\mp}\to\ell^{\pm}Z/h\nu
W^{\mp}\to\ell^{\pm}\ell^{+}\ell^{-}\nu J_{W}$,

  2. (b)

$\Sigma^{\pm}\Sigma^{0}\to\ell^{\pm}Z/h\nu Z/h\to\ell^{\pm}\ell^{+}\ell^{-}\nu
J_{Z/h}$,

  3. (c)

$\Sigma^{\pm}\Sigma^{0}\to\ell^{\pm}Z/h\ell^{\pm(\mp)}W^{\mp(\pm)}\to\ell^{\pm}\ell^{\pm}\ell^{\mp}\nu
J_{Z/h}$,

For processes (a) and (b), $\Sigma^{\pm}$ decays to $\ell^{\pm}$ and $Z/h$
which further decays to $\ell^{+}\ell^{-}$. For process (c), $\Sigma^{\pm}$
decays to $\ell^{\pm}$ and $J_{Z/h}$. Therefore, for the events with an on-$Z$
lepton pair,777An opposite-sign same flavour lepton pair with an invariant
mass within the $Z$-boson mass window, i.e. $m_{Z}\pm 15$ GeV is termed as
“on-$Z$ lepton pair”. the peak of the invariant mass distribution of the
trilepton system is expected to yield $m_{\Sigma}$, whereas, for those with no
on-$Z$ lepton pair, $m_{\Sigma}$ can be reconstructed from the invariant mass
distribution of $J_{Z/h}$ and the farthest lepton in the azimuth plane.

$(vi)$ Tetralepton with/without one fat-jet: Both pair and associated
production of charged triplets contribute to this final state:

  1. (a)

$\Sigma^{+}\Sigma^{-}\to\ell^{+}Z/h\ell^{-}Z/h\to\ell^{+}\ell^{+}\ell^{-}\ell^{-}J_{Z/h}$.

  2. (b)

$\Sigma^{\pm}\Sigma^{\mp}\to\ell^{\pm}Z/h\nu
W^{\mp}\to\ell^{\pm}\ell^{+}\ell^{-}\ell^{\mp}\nu\nu$.

  3. (c)

$\Sigma^{\pm}\Sigma^{0}\to\ell^{\pm}Z/h\ell^{\pm(\mp)}W^{\mp(\pm)}\to\ell^{\pm}\ell^{+}\ell^{-}\ell^{\pm(\mp)}J_{W}$.

All the events necessarily contain an on-$Z$ lepton pair. The invariant mass
distribution of the on-$Z$ lepton pair and the farthest lepton between the
ones which are not part of the on-$Z$ lepton pair would peak at $m_{\Sigma}$.

We next give a brief description of reconstruction and selection of various
physics objects, classification of events into several mutually exclusive
analysis channels and event selection.

###  III.1 Object reconstruction and selection

We use Delphes deFavereau:2013fsa  to reconstruct different objects, namely
photons, electrons, muons and jets. Constituents of the would-be fat-jets are
clustered using the anti-kT algorithm Cacciari:2008gp  with the winner-take-
all axis with the characteristic fat-jet radius $R=0.8$ as implemented in
FastJet Cacciari:2011ma . We use the jet pruning algorithm Ellis:2009su ;
Ellis:2009me , setting the parameters to their default values, viz.
$z_{cut}=0.1$ and $R_{cut}=0.5$ Ellis:2009su , to remove the softer and wide-
angle QCD emissions from the fat-jets. Further, an inclusive jet shape termed
as $N$-subjettiness, $\tau_{N}$, Thaler:2010tr ; Thaler:2011gf  is used to
unfold the multi-prong nature of the fat-jets. We choose the one-pass
$k_{T}$-axes for the minimization procedure, and take the thurst parameter
$\beta=1$. The reconstructed leptons (electrons and muons) and fat-jets are
required to have transverse momentum $p_{T}>10$ and 30 GeV, respectively, and
they must lie within the tracking system acceptance, pseudorapidity
$|\eta|<2.5$. The relative isolation888The relative isolation of a lepton is
defined as the scalar $p_{T}$ sum, normalized to the lepton $p_{T}$, of
photons and hadrons within a cone of $\Delta R=0.5$ around the lepton. is
required to be smaller than 12%(15%) for electrons (muons). Such stringent
lepton isolation requirements significantly suppress the reducible
backgrounds. The selected leptons and jets are ordered by decreasing $p_{T}$
($\ell_{0},\ell_{1},...$ and $J_{0},J_{1},...$). Often, some of the jets are
misidentified as leptons. Following the Ref. ATLAS:2016iqc , we take the
probability of 0.1–0.3% for a jet to be misidentified as a lepton.
Furthermore, bremsstrahlung interactions of electrons with the inner detector
material often lead to their charge misidentification. We adopt the charge
misidentification probability from Ref. ATLAS:2017xqs :
$P(p_{T},\eta)=\sigma(p_{T})\times f(\eta)$, where $\sigma(p_{T})$ and
$f(\eta)$ ranges from 0.02 to 0.1 and 0.03 to 1, respectively. Finally, the
missing transverse momentum vector $\vec{p}_{T}^{\rm~{}miss}$ (with magnitude
$p_{T}^{\rm miss}$) is estimated using all the reconstructed particle-flow
objects in an event.

###  III.2 SM Backgrounds

We next briefly mention the SM processes which could resemble the
abovementioned final states. The relevant processes can be classified into two
classes — reducible and irreducible backgrounds. The SM processes giving rise
to one or more prompt leptons, such as Drell-Yan processes, $VV$
($V=W,Z/\gamma^{*}$), $VVV$, $VVVV$, $t\bar{t}$, $t\bar{t}V$, single top
production ($tb,tW,tj$) and Higgsstrahlung ($Vh,VVh,t\bar{t}h$) processes,
constitute the irreducible backgrounds in this analysis. On the contrary, the
reducible backgrounds are from the SM processes like $Z/\gamma^{*}+$jets and
$t\bar{t}+$jets, where a jet is misidentified as lepton or additional leptons
originate from heavy quark decays. The reducible backgrounds are significantly
suppressed by applying stringent lepton isolation requirements. All the
background samples are generated in association with up to two jets using
MadGraph Alwall:2011uj ; Alwall:2014hca  at the LO using the 5 flavour scheme
followed by MLM matching in PYTHIA Sjostrand:2014zea . The background samples
are normalised at least to the NLO prediction Campbell:1999ah ;
Ciccolini:2003jy ; Brein:2003wg ; Catani:2007vq ; Campanario:2008yg ;
Balossini:2009sa ; Bredenstein:2009aj ; Catani:2009sm ; Kidonakis:2010ux ;
Campbell:2011bn ; Brein:2011vx ; Bevilacqua:2012em ; Garzelli:2012bn ;
Brein:2012ne ; Altenkamp:2012sx ; Nhung:2013jta ; Kidonakis:2013zqa ;
Denner:2014cla ; Harlander:2014wda ; Kidonakis:2015nna ; Muselli:2015kba ;
Shen:2015cwj ; Frederix:2017wme .

###  III.3  Event selection

After object selection, we make our identification of the $W/Z/h$-fat-jets,
denoted as $J_{W/Z/h}$, satisfying the conditions:999The N-subjettiness,
$\tau_{N}$, is a good measure of the number of subjets a jet is presumably
composed of. The ratio $\tau_{N}/\tau_{N-1}$ is, thus, an useful discriminant
between the $N$\- and $(N-1)$-prong jets. Being two-prong in nature, the jets
emanating from the sufficiently highly boosted $W/Z/h$-bosons tend to have
lower $\tau_{21}=\frac{\tau_{2}}{\tau_{1}}$ as compared to the overwhelming
one-prong QCD or top jets. The cut on $\tau_{21}$ significantly suppresses the
background contributions from the QCD or top-jets.

$m(J)\in[70,140]{\rm~{}GeV,~{}}p_{T}(J)>200{\rm~{}GeV~{}and~{}}\tau_{21}<0.7.$

Based on the number of leptons and $J_{W/Z/h}$ ($n_{\ell}$ and
$n_{J_{W/Z/h}}$) and the absolute value of the sum of charges of the leptons
$\left(|\sum Q_{\ell}|\right)$, the events are categorised into seven distinct
final states — 1L-1J, 1L-2J, SSD-1J, OSD-1J, 3L-0J, 3L-1J and 4L (see Table
1), while the events not falling under these channels are thrown away. As
mentioned earlier, the leptons and bosons coming from the decays of the TeV
sacle triplets would carry sufficiently high momenta that stronger cuts on the
leptons’ $p_{T}$ and/or $p_{T}^{\rm miss}$ would be utile in curtailing the SM
background while keeping the signal all but unharmed. The event selection
proceeds in two steps: the preselection and the channel-specific event
selection. The preselection requirements are primarily based on
$p_{T}(\ell_{0})$, $p_{T}(\ell_{1})$ and $p_{T}^{\rm miss}$, see Table 1.
Events containing an opposite-sign same-flavour (OSSF) lepton pair with
invariant mass below 12 GeV are thrown away in all the channels to reduce
background contributions from low-mass resonances. Furthermore, to suppress
the overwhelming background contributions from the Drell-Yan processes in the
SSD-1J and OSD-2J channels, events containing an opposite-sign same-flavour
lepton pair with an invariant mass within the nominal $Z$-boson mass window,
i.e. $m_{Z}\pm 15$ GeV are vetoed.101010Note that, in order to reduce the
contributions originating from electron charge misidentification, the Z-veto
is also applied to the SSD-1J events.

Selection criteria | 1L-1J | 1L-2J | SSD-1J | OSD-1J | 3L-0J | 3L-1J | 4L
---|---|---|---|---|---|---|---
$n_{\ell}$ | $=1$ | $=1$ | $=2$ | $=2$ | $=3$ | $=3$ | $=4$
$|\sum Q_{\ell}|$ | $=1$ | $=1$ | $=2$ | $=0$ | – | – | –
$n_{J_{W/Z/h}}$ | $=1$ | $=2$ | $\geq 1$ | $\geq 1$ | $=0$ | $=1$ | $\geq 0$
$p_{T}(\ell_{0})$ |  $>$ 200 |  $>$ 200 |  $>$ 400 |  $>$ 400 |  $>$ 300 |  $>$ 300 |  $>$ 400
$p_{T}(\ell_{1})$ | – | – |  $>$ 100 |  $>$ 100 |  $>$ 100 |  $>$ 100 |  $>$ 200
$p_{T}^{\rm miss}$ |  $>$ 100 |  $>$ 100 | – | – |  $>$ 100 | – | –
Low $m_{\ell\ell}^{\rm OSSF}$ veto | – | – | $>12$ | $>12$ | $>12$ | $>12$ | $>12$
$Z$-veto | – | – | Yes | Yes | – | – | –

Table 1:  The preselection criteria for the considered analysis channels. All
the selection cuts are in GeV. The symbol “–” means no requirement is made.

Figure 1: Normalised kinematic distributions for signal ($m_{\Sigma}=1150$
GeV) and background in the 1L-1J channel. Top panel: lepton $p_{T}$ (left) and
$p_{T}^{\rm miss}$ (right); middle panel: $H_{T}$ (left) and $m_{\rm eff}$
(right); bottom panel: $\Delta\phi(\ell_{0},p_{T}^{\rm miss})$ (left) and
$W/Z/h$-fat-jet $p_{T}$ (right).

We next briefly discuss the channel-specific selection criteria, which
appreciably reduce the remaining SM backgrounds without impinging much on the
signal strength. We use various kinematic distributions as guiding premises to
choose the appropriate selection cuts.

We display different normalised kinematic distributions for the 1L-1J events
in Figure 1 for a benchmark point: $m_{\Sigma}=1150$ GeV. In the top panel, we
show the $p_{T}(\ell_{0})$ (left) and $p_{T}^{\rm miss}$ (right)
distributions. The middle panel displays the $H_{T}$ (left) and $m_{\rm
eff}=L_{T}+H_{T}+p_{T}^{\rm miss}$ (right) distributions, where $L_{T}(H_{T})$
is the scalar sum of all the leptons’ (jets’) $p_{T}$. In the bottom panel,
the left plot shows the distribution of the azimuthal separation between
$\ell_{0}$ and $p_{T}^{\rm miss}$, $|\Delta\phi(\ell_{0},p_{T}^{\rm miss})|$,
whereas the right plot shows the $p_{T}(J_{0})$ distribution. As expected, the
kinematic distributions for the signal are much harder than those for the
background, see Figure 1. Consequently, harder cuts on these kinematic
variables are expected to suppress the remaining background without impinging
much on the signal strength. On the other hand,
$|\Delta\phi(\ell_{0},p_{T}^{\rm miss})|$ turns out to be efficient in
discriminating the background (in particular, the background from single
W-boson) and signal. We impose the following selection cuts:

SI-1: $\displaystyle p_{T}(\ell_{0})>400{\rm~{}GeV~{}and~{}}p_{T}^{\rm
miss}>400{\rm~{}GeV},$ SI-2: $\displaystyle H_{T}>500{\rm~{}GeV,~{}}m_{\rm
eff}>1800{\rm~{}GeV~{}and~{}}$ $\displaystyle p_{T}(J_{0})>300{\rm~{}GeV},$
SI-3: $\displaystyle|\Delta\phi(\ell_{0},p_{T}^{\rm miss})|>0.5~{}.$

Figure 2: Normalised distribution of $\tau_{21}(J_{0})$ after passing the
S1-1, S1-2 and S1-3 cuts for signal ($m_{\Sigma}=1150$ GeV) and background in
the 1L-1J channel.

Also displayed, in Figure 2, is the normalised distribution of
$\tau_{21}(J_{0})$ for the 1L-1J events passing the S1-1, S1-2 and S1-3 cuts.
As expected, an inordinately large fraction of the signal events have lower
$\tau_{21}$ as compared to the remaining background events. Therefore, to
further improve the signal-to-background ratio, we impose the following
selection cut:

SI-4: $\displaystyle\tau_{21}(J_{0})<0.5~{}.$

We display the effect of the selection cuts on the signal and background
strengths for the 1L-1J channel in Table 2. It shows the progression of the
number of expected background and signal ($m_{\Sigma}=1150$ GeV) events as the
subsequent selection cuts are imposed for 1000 fb-1 luminosity at the 13 TeV
LHC. As we expected, the preselection cuts and channel-specific selection cuts
effectively vanquish the SM backgrounds while keeping a large fractions of the
signal.

1L-1J channel
---
Event sample | Preselection |  SI-1 |  SI-2 |  SI-3 |  SI-4
Drell-Yan | $\sim 2.9\times 10^{5}$ |  5119 |  4057 |  221 |  98
$t\bar{t}$ | $\sim 1.7\times 10^{5}$ |  291 |  190 |  44 |  20
Other backgrounds | $\sim 1.7\times 10^{4}$ |  239 |  185 |  25 |  14
Total background | $\sim 4.8\times 10^{5}$ |  5649 |  4432 |  290 |  132
Signal | 330 |  126 |  93 |  90 |  71

Table 2:  The number of expected background and signal ($m_{\Sigma}=1150$ GeV)
events in the 1L-1J channel after passing various cuts for 1000 fb-1
luminosity at the 13 TeV LHC.

After discussing the cuts for the 1L-1J channel in detail, we next summarise
the channel specific selection cuts for all other channels in Table 3. Table 3
also shows the number of expected background and signal ($m_{\Sigma}=1150$
GeV) events after all the selection cuts are imposed for 1000 fb-1 luminosity
at the 13 TeV LHC.

Channel |  Selection cuts | Background | Signal
---|---|---|---
1L-2J |  $p_{T}(\ell_{0})>400$, $p_{T}^{\rm miss}>300$, | 74 | 63
|  $|\Delta\phi(\ell_{0},p_{T}^{\rm miss})|>0.5$ |  |
SSD-1J |  $L_{T}>600$, $H_{T}>500$, $m_{\rm eff}>1600$, | 81 | 32
|  $p_{T}(J_{0})>300$ |  |
OSD-1J |  $L_{T}>700$, $H_{T}>500$, $m_{\rm eff}>2000$, | 131 | 46
|  $p_{T}(J_{0})>400$, $\tau_{21}(J_{0})<0.6$ |  |
3L-0J |  $L_{T}>600$, $m_{\rm eff}>1500$, $L_{T}/H_{T}>1.0$, | 67 | 24
|  $p_{T}^{\rm miss}/H_{T}>0.4$ |  |
3L-1J |  $L_{T}>700$, $m_{\rm eff}>1600$, $p_{T}(J_{0})>300$, | 28 | 20
|  $L_{T}/H_{T}>1.0$ |  |
4L |  $L_{T}>1000$, $m_{\rm eff}>1600$ | 9 | 8

Table 3:  Summary of the channel-specific selection cuts along with the number
of expected background and signal ($m_{\Sigma}=1150$ GeV) events after all the
selection cuts are imposed for 1000 fb-1 luminosity at the 13 TeV LHC. All the
dimensionful cuts are in GeV.

In order to enhance the sensitivity of this search, the selected events in
each channel are divided into several independent search bins using a primary
kinematic discriminant between the signal and background. The discriminants
used in different channels are as follows:

  * •

1L-1J: the invariant mass of the lepton and the fat-jet: $m_{\ell
J}=m_{\ell_{0}J_{0}}$;

  * •

1L-2J: the invariant mass of the lepton and the farthest fat-jet in the
azimuth plane:

$\displaystyle m_{\ell J}$
$\displaystyle=m_{\ell_{0}J_{0}}{\rm~{}if~{}}\Delta\phi(\ell_{0},J_{0})>\Delta\phi(\ell_{0},J_{1})$
$\displaystyle=m_{\ell_{0}J_{1}}{\rm~{}otherwise};$

  * •

SSD-1J and OSD-1J: $(i)$ for events with exactly one $J_{W/Z/h}$: the
invariant mass of the fat-jet and the farthest lepton in the azimuth plane:

$\displaystyle m_{\ell J}$
$\displaystyle=m_{\ell_{0}J_{0}}{\rm~{}if~{}}\Delta\phi(\ell_{0},J_{0})>\Delta\phi(\ell_{1},J_{0})$
$\displaystyle=m_{\ell_{1}J_{0}}{\rm~{}otherwise};$

$(ii)$ for events with more than one $J_{W/Z/h}$:

$\displaystyle\quad\qquad m_{\ell J}$
$\displaystyle=m_{\ell_{0}J_{0}}{\rm~{}if~{}}|m_{\ell_{0}J_{0}}-m_{\ell_{1}J_{1}}|<|m_{\ell_{0}J_{1}}-m_{\ell_{1}J_{0}}|$
$\displaystyle=m_{\ell_{0}J_{1}}{\rm~{}otherwise};$

  * •

3L-0J: Missing neutrinos do not allow kinematic reconstruction of the triplets
in this channel, and, thus, we use $L_{T}+p_{T}^{\rm miss}$ as primary
kinematic discriminant;

  * •

3L-1J: $(i)$ for events with one on-$Z$ lepton pair, the invariant mass of the
trilepton system: $m_{\ell_{0}\ell_{1}\ell_{2}}$;
$(ii)$ for events with no on-$Z$ lepton pair: the invariant mass of the fat-
jet and the farthest lepton in the azimuth plane;

  * •

4L: the invariant mass of the on-$Z$ lepton pair and the farthest lepton in
the azimuth plane: $m_{(\ell\ell)_{Z}\ell}$.

For all the channels except the 3L-0J one, the distribution of the primary
discriminant peaks at the triplet mass and thus could be used to reconstruct
the same. We summarise the binning scheme in Table 4, which yields a total of
50 statistically independent search bins. The width of the bins is chosen to
ensure smooth behaviour for the expected background and signal distributions.
The underflow (overflow) events are contained in the first (last) bin for each
channel.

Channel | Kinematic discriminant | Range (GeV) | Number of bins
---|---|---|---
1L-1J | $m_{\ell_{0}J_{0}}$ | $[600,2000]$ | 7
1L-2J |  $m_{\ell_{0}J_{0}}$ or $m_{\ell_{0}J_{1}}$ | $[600,2000]$ | 7
SSD-1J |  $m_{\ell_{0}J_{0}}$ or $m_{\ell_{1}J_{0}}$ | $[600,2000]$ | 7
OSD-1J |  $m_{\ell_{0}J_{0}}$ or $m_{\ell_{1}J_{0}}$ | $[600,2000]$ | 7
3L-0J | $L_{T}+p_{T}^{\rm miss}$ | $[750,2750]$ | 8
3L-1J |  $m_{\ell_{0}J_{0}}$ or $m_{\ell_{1}J_{0}}$ or $m_{\ell_{0}\ell_{1}\ell_{2}}$ | $[600,2000]$ | 7
4L | $m_{(\ell\ell)_{Z}\ell}$ | $[600,2000]$ | 7

Table 4: Binning scheme of different channels.

Figure 3 displays the distributions of the primary kinematic discriminants for
the signal ($m_{\Sigma}=1150$ GeV) and background events for all the channels.
The events are weighted for 1000 fb-1 luminosity at the 13 TeV LHC. As
expected, these distributions (except the one for 3L-0J channel) peak in
1000–1200 GeV bin, and thereby reconstructing the triplet mass; a simple
rebinning of these distributions with smaller bin widths would make this mass
reconstruction apparent.

Figure 3: Distributions of the primary kinematic discriminants for for the
signal ($m_{\Sigma}=1150$ GeV) and background. The events are weighted for
1000 fb-1 luminosity at the 13 TeV LHC.

###  III.4 Discovery reach

We next estimate the discovery reach of the present search for the triplet
fermions. To this end, we use a hypothesis tester named Profile Likelihood
Number Counting Combination which uses the library RooFit Verkerke:2003ir  in
the ROOT Brun:1997pa  environment. This hypothesis tester considers all the
search bins as independent number counting channels with uncorrelated
systematic errors. The uncertainties are included via the Profile Likelihood
Ratio. For the sake of simplicity, we assume an overall 20% total uncertainty
on the estimated background. Further, to avoid spurious exclusions/discoveries
and to ensure robustness in statistical interpretations, we replace the less
than one per-bin expected background yield at 3000 fb-1 with one background
yield. This renders our estimated significances to be a bit conservative. In
Figure 4, we project the required luminosities for $3\sigma$ discovery of the
triplet fermions in different analysis channels as a function of their mass;
for brevity, we display the same for the most promising channels only. We see
that the 1L-2J channel is the most promising one with the discovery reach of
1480 GeV, while the rest of the channels have a similar discovery reach of
1330-1350 GeV.111111Note that the 1L-channels are not suppressed by the
leptonic branching fractions of the SM bosons, and thus, yield a sufficient
number of signal events even for larger triplet mass. Moreover, not only does
the requirement of two $J_{W/Z/h}$ in the final state vanquish the background,
but it also allows kinematic reconstruction of the triplets without any
ambiguity (see the top panel in Figure 3). On the contrary, final states with
higher lepton multiplicity (three or more) being suppressed by the leptonic
branching fractions of the SM bosons are less sensitive in probing triplets
with higher masses. Note that the recent CMS/ATLAS searches targetting the
triplet fermions in the type-III see-saw model rely on multi-lepton final
states only, and thus, are not sensitive enough in probing them with higher
masses.

Figure 4: Required Luminosity for $3\sigma$ discovery of the triplet fermions
in different analysis channels.

Since all the analysis channels considered in this search are mutually
exclusive, it is reasonable to combine them. For the combined channels, we
project the required luminosities for both $3\sigma$ and $5\sigma$ discoveries
of the triplet fermions as a function of their mass in Figure 5. The $5\sigma$
($3\sigma$) discovery reach of this search is about 1265 (1380) and 1480
(1600) GeV, respectively, at 500 and 3000 fb-1 luminosities.

Figure 5: Required Luminosity for $3\sigma$ and $5\sigma$ discovery of the
triplet fermions in the combined channels.

##  IV Summary

The CMS and ATLAS collaborations have performed several searches targetting
the triplet fermions in the type-III see-saw model. However, the final states
considered in these searches being beset with considerably large SM
backgrounds are deemed not sensitive enough to probe the triplet fermions much
heavier than 1 TeV. To this end, we perform a search for the triplet fermions
in final states with multiple leptons and fat-jets that are cleaner than the
usual LHC searches and allow reconstruction of the triplet mass. After
performing a systematic and comprehensive search with seven distinct analysis
channels, we project the required luminosities for both $3\sigma$ and
$5\sigma$ discoveries of the triplet fermions as a function of their mass. The
triplet fermions with mass as large as 1265 (1380) and 1480 (1600) GeV could
be discovered with $5\sigma$ ($3\sigma$) significance at 500 and 3000 fb-1
luminosities, respectively. 121212Thogh single production of the triplet
fermions, in principle, in association with a SM lepton is possible, such
processes are extremely suppressed both at $pp$ and $e^{-}e^{+}$ colliders on
account of their small mixing with the SM leptons. In view of pair production
of the triplet fermions, in the most optimistic scenario, a high energy
$e^{-}e^{+}$ collider could probe them very close to its kinematic reach
Goswami:2017jqs ; Das:2020gnt , i.e. 1500 GeV for a 3 TeV collider. Note that
the $3\sigma$ ($5\sigma$) reach of our search at 3000 fb-1 luminosity extends
beyond (very close to) the kinematic reach of a 3 TeV $e^{-}e^{+}$ collider.
In closing, we mention that a search similar to the one presented in this work
is prophesied to carry through for a large class of neutrino mass models
containing new weak gauge multiplets. We undertake a detailed analysis for the
same in a future work Ashanujjaman:toappear

###### Acknowledgements.

KG acknowledges support from the DST INSPIRE Research Grant
[DST/INSPIRE/04/2014/002158] and SERB Core Research Grant [CRG/2019/006831].
The simulations were supported in part by the SAMKHYA: High Performance
Computing Facility provided by Institute of Physics, Bhubaneswar.

## References

  * (1) S. Weinberg, Phys. Rev. Lett. 43, 1566-1570 (1979) doi:10.1103/PhysRevLett.43.1566
  * (2) E. Ma, Phys. Rev. Lett. 81, 1171-1174 (1998) doi:10.1103/PhysRevLett.81.1171 [arXiv:hep-ph/9805219 [hep-ph]].
  * (3) R. Foot, H. Lew, X. G. He and G. C. Joshi, Z. Phys. C 44, 441 (1989) doi:10.1007/BF01415558
  * (4) B. Bajc and G. Senjanovic, JHEP 08, 014 (2007) doi:10.1088/1126-6708/2007/08/014 [arXiv:hep-ph/0612029 [hep-ph]].
  * (5) P. Fileviez Perez, Phys. Lett. B 654, 189-193 (2007) doi:10.1016/j.physletb.2007.07.075 [arXiv:hep-ph/0702287 [hep-ph]].
  * (6) P. Fileviez Perez, Phys. Rev. D 76, 071701 (2007) doi:10.1103/PhysRevD.76.071701 [arXiv:0705.3589 [hep-ph]].
  * (7) F. del Aguila and J. A. Aguilar-Saavedra, Phys. Lett. B 672, 158-165 (2009) doi:10.1016/j.physletb.2009.01.010 [arXiv:0809.2096 [hep-ph]].
  * (8) F. del Aguila and J. A. Aguilar-Saavedra, Nucl. Phys. B 813, 22-90 (2009) doi:10.1016/j.nuclphysb.2008.12.029 [arXiv:0808.2468 [hep-ph]].
  * (9) R. Franceschini, T. Hambye and A. Strumia, Phys. Rev. D 78, 033002 (2008) doi:10.1103/PhysRevD.78.033002 [arXiv:0805.1613 [hep-ph]].
  * (10) A. Arhrib, B. Bajc, D. K. Ghosh, T. Han, G. Y. Huang, I. Puljak and G. Senjanovic, Phys. Rev. D 82, 053004 (2010) doi:10.1103/PhysRevD.82.053004 [arXiv:0904.2390 [hep-ph]].
  * (11) P. Bandyopadhyay, S. Choi, E. J. Chun and K. Min, Phys. Rev. D 85, 073013 (2012) doi:10.1103/PhysRevD.85.073013 [arXiv:1112.3080 [hep-ph]].
  * (12) C. Biggio and F. Bonnet, Eur. Phys. J. C 72, 1899 (2012) doi:10.1140/epjc/s10052-012-1899-z [arXiv:1107.3463 [hep-ph]].
  * (13) D. Goswami and P. Poulose, Eur. Phys. J. C 78, no.1, 42 (2018) doi:10.1140/epjc/s10052-017-5478-1 [arXiv:1702.07215 [hep-ph]].
  * (14) Y. Cai, T. Han, T. Li and R. Ruiz, Front. in Phys. 6, 40 (2018) doi:10.3389/fphy.2018.00040 [arXiv:1711.02180 [hep-ph]].
  * (15) A. Das, S. Mandal and T. Modak, Phys. Rev. D 102, no.3, 033001 (2020) doi:10.1103/PhysRevD.102.033001 [arXiv:2005.02267 [hep-ph]].
  * (16) S. Jana, N. Okada and D. Raut, [arXiv:1911.09037 [hep-ph]].
  * (17) A. Das and S. Mandal, Nucl. Phys. B 966, 115374 (2021) doi:10.1016/j.nuclphysb.2021.115374 [arXiv:2006.04123 [hep-ph]].
  * (18) S. Ashanujjaman and K. Ghosh, Phys. Lett. B 819, 136403 (2021) doi:10.1016/j.physletb.2021.136403 [arXiv:2102.09536 [hep-ph]].
  * (19) C. Sen, P. Bandyopadhyay, S. Dutta and A. KT, [arXiv:2107.12442 [hep-ph]].
  * (20) G. Aad et al. [ATLAS], Phys. Rev. D 92, no.3, 032001 (2015) doi:10.1103/PhysRevD.92.032001 [arXiv:1506.01839 [hep-ex]].
  * (21) A. M. Sirunyan et al. [CMS], Phys. Rev. Lett. 119, no.22, 221802 (2017) doi:10.1103/PhysRevLett.119.221802 [arXiv:1708.07962 [hep-ex]].
  * (22) [ATLAS], ATLAS-CONF-2018-020.
  * (23) A. M. Sirunyan et al. [CMS], JHEP 03, 051 (2020) doi:10.1007/JHEP03(2020)051 [arXiv:1911.04968 [hep-ex]].
  * (24) G. Aad et al. [ATLAS], Eur. Phys. J. C 81, no.3, 218 (2021) doi:10.1140/epjc/s10052-021-08929-9 [arXiv:2008.07949 [hep-ex]].
  * (25) J. A. Casas and A. Ibarra, Nucl. Phys. B 618, 171-204 (2001) doi:10.1016/S0550-3213(01)00475-8 [arXiv:hep-ph/0103065 [hep-ph]].
  * (26) A. Ibarra and G. G. Ross, Phys. Lett. B 591, 285-296 (2004) doi:10.1016/j.physletb.2004.04.037 [arXiv:hep-ph/0312138 [hep-ph]].
  * (27) I. Esteban, M. C. Gonzalez-Garcia, A. Hernandez-Cabezudo, M. Maltoni and T. Schwetz, JHEP 01, 106 (2019) doi:10.1007/JHEP01(2019)106 [arXiv:1811.05487 [hep-ph]].
  * (28) F. Staub, Comput. Phys. Commun. 185, 1773-1790 (2014) doi:10.1016/j.cpc.2014.02.018 [arXiv:1309.7223 [hep-ph]].
  * (29) F. Staub, Adv. High Energy Phys. 2015, 840780 (2015) doi:10.1155/2015/840780 [arXiv:1503.04200 [hep-ph]].
  * (30) J. Alwall, M. Herquet, F. Maltoni, O. Mattelaer and T. Stelzer, JHEP 06, 128 (2011) doi:10.1007/JHEP06(2011)128 [arXiv:1106.0522 [hep-ph]].
  * (31) J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli and M. Zaro, JHEP 07, 079 (2014) doi:10.1007/JHEP07(2014)079 [arXiv:1405.0301 [hep-ph]].
  * (32) A. Manohar, P. Nason, G. P. Salam and G. Zanderighi, Phys. Rev. Lett. 117, no.24, 242002 (2016) doi:10.1103/PhysRevLett.117.242002 [arXiv:1607.04266 [hep-ph]].
  * (33) A. V. Manohar, P. Nason, G. P. Salam and G. Zanderighi, JHEP 12, 046 (2017) doi:10.1007/JHEP12(2017)046 [arXiv:1708.01256 [hep-ph]].
  * (34) J. Butterworth, S. Carrazza, A. Cooper-Sarkar, A. De Roeck, J. Feltesse, S. Forte, J. Gao, S. Glazov, J. Huston and Z. Kassabov, et al. J. Phys. G 43, 023001 (2016) doi:10.1088/0954-3899/43/2/023001 [arXiv:1510.03865 [hep-ph]].
  * (35) R. Ruiz, JHEP 12, 165 (2015) doi:10.1007/JHEP12(2015)165 [arXiv:1509.05416 [hep-ph]].
  * (36) M. Cirelli, N. Fornengo and A. Strumia, Nucl. Phys. B 753, 178-194 (2006) doi:10.1016/j.nuclphysb.2006.07.012 [arXiv:hep-ph/0512090 [hep-ph]].
  * (37) A. Abada, C. Biggio, F. Bonnet, M. B. Gavela and T. Hambye, Phys. Rev. D 78, 033007 (2008) doi:10.1103/PhysRevD.78.033007 [arXiv:0803.0481 [hep-ph]].
  * (38) S. Goswami, K. N. Vishnudath and N. Khan, Phys. Rev. D 99, no.7, 075012 (2019) doi:10.1103/PhysRevD.99.075012 [arXiv:1810.11687 [hep-ph]].
  * (39) C. Biggio, E. Fernandez-Martinez, M. Filaci, J. Hernandez-Garcia and J. Lopez-Pavon, JHEP 05, 022 (2020) doi:10.1007/JHEP05(2020)022 [arXiv:1911.11790 [hep-ph]].
  * (40) J. de Favereau et al. [DELPHES 3], JHEP 02, 057 (2014) doi:10.1007/JHEP02(2014)057 [arXiv:1307.6346 [hep-ex]].
  * (41) M. Cacciari, G. P. Salam and G. Soyez, JHEP 04, 063 (2008) doi:10.1088/1126-6708/2008/04/063 [arXiv:0802.1189 [hep-ph]].
  * (42) M. Cacciari, G. P. Salam and G. Soyez, Eur. Phys. J. C 72, 1896 (2012) doi:10.1140/epjc/s10052-012-1896-2 [arXiv:1111.6097 [hep-ph]].
  * (43) S. D. Ellis, C. K. Vermilion and J. R. Walsh, Phys. Rev. D 80, 051501 (2009) doi:10.1103/PhysRevD.80.051501 [arXiv:0903.5081 [hep-ph]].
  * (44) S. D. Ellis, C. K. Vermilion and J. R. Walsh, Phys. Rev. D 81, 094023 (2010) doi:10.1103/PhysRevD.81.094023 [arXiv:0912.0033 [hep-ph]].
  * (45) J. Thaler and K. Van Tilburg, JHEP 03, 015 (2011) doi:10.1007/JHEP03(2011)015 [arXiv:1011.2268 [hep-ph]].
  * (46) J. Thaler and K. Van Tilburg, JHEP 02, 093 (2012) doi:10.1007/JHEP02(2012)093 [arXiv:1108.2701 [hep-ph]].
  * (47) [ATLAS], ATLAS-CONF-2016-024.
  * (48) M. Aaboud et al. [ATLAS], Eur. Phys. J. C 78, no.3, 199 (2018) doi:10.1140/epjc/s10052-018-5661-z [arXiv:1710.09748 [hep-ex]].
  * (49) T. Sjöstrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten, S. Mrenna, S. Prestel, C. O. Rasmussen and P. Z. Skands, Comput. Phys. Commun. 191, 159-177 (2015) doi:10.1016/j.cpc.2015.01.024 [arXiv:1410.3012 [hep-ph]].
  * (50) J. M. Campbell and R. K. Ellis, Phys. Rev. D 60, 113006 (1999) doi:10.1103/PhysRevD.60.113006 [arXiv:hep-ph/9905386 [hep-ph]].
  * (51) M. L. Ciccolini, S. Dittmaier and M. Kramer, Phys. Rev. D 68, 073003 (2003) doi:10.1103/PhysRevD.68.073003 [arXiv:hep-ph/0306234 [hep-ph]].
  * (52) O. Brein, A. Djouadi and R. Harlander, Phys. Lett. B 579, 149-156 (2004) doi:10.1016/j.physletb.2003.10.112 [arXiv:hep-ph/0307206 [hep-ph]].
  * (53) S. Catani and M. Grazzini, Phys. Rev. Lett. 98, 222002 (2007) doi:10.1103/PhysRevLett.98.222002 [arXiv:hep-ph/0703012 [hep-ph]].
  * (54) F. Campanario, V. Hankele, C. Oleari, S. Prestel and D. Zeppenfeld, Phys. Rev. D 78, 094012 (2008) doi:10.1103/PhysRevD.78.094012 [arXiv:0809.0790 [hep-ph]].
  * (55) G. Balossini, G. Montagna, C. M. Carloni Calame, M. Moretti, O. Nicrosini, F. Piccinini, M. Treccani and A. Vicini, JHEP 01, 013 (2010) doi:10.1007/JHEP01(2010)013 [arXiv:0907.0276 [hep-ph]].
  * (56) A. Bredenstein, A. Denner, S. Dittmaier and S. Pozzorini, Phys. Rev. Lett. 103, 012002 (2009) doi:10.1103/PhysRevLett.103.012002 [arXiv:0905.0110 [hep-ph]].
  * (57) S. Catani, L. Cieri, G. Ferrera, D. de Florian and M. Grazzini, Phys. Rev. Lett. 103, 082001 (2009) doi:10.1103/PhysRevLett.103.082001 [arXiv:0903.2120 [hep-ph]].
  * (58) N. Kidonakis, Phys. Rev. D 82, 054018 (2010) doi:10.1103/PhysRevD.82.054018 [arXiv:1005.4451 [hep-ph]].
  * (59) J. M. Campbell, R. K. Ellis and C. Williams, JHEP 07, 018 (2011) doi:10.1007/JHEP07(2011)018 [arXiv:1105.0020 [hep-ph]].
  * (60) O. Brein, R. Harlander, M. Wiesemann and T. Zirke, Eur. Phys. J. C 72, 1868 (2012) doi:10.1140/epjc/s10052-012-1868-6 [arXiv:1111.0761 [hep-ph]].
  * (61) G. Bevilacqua and M. Worek, JHEP 07, 111 (2012) doi:10.1007/JHEP07(2012)111 [arXiv:1206.3064 [hep-ph]].
  * (62) M. V. Garzelli, A. Kardos, C. G. Papadopoulos and Z. Trocsanyi, JHEP 11, 056 (2012) doi:10.1007/JHEP11(2012)056 [arXiv:1208.2665 [hep-ph]].
  * (63) O. Brein, R. V. Harlander and T. J. E. Zirke, Comput. Phys. Commun. 184, 998-1003 (2013) doi:10.1016/j.cpc.2012.11.002 [arXiv:1210.5347 [hep-ph]].
  * (64) L. Altenkamp, S. Dittmaier, R. V. Harlander, H. Rzehak and T. J. E. Zirke, JHEP 02, 078 (2013) doi:10.1007/JHEP02(2013)078 [arXiv:1211.5015 [hep-ph]].
  * (65) D. T. Nhung, L. Ninh and M. M. Weber, JHEP 12, 096 (2013) doi:10.1007/JHEP12(2013)096 [arXiv:1307.7403 [hep-ph]].
  * (66) N. Kidonakis, doi:10.3204/DESY-PROC-2013-03/Kidonakis [arXiv:1311.0283 [hep-ph]].
  * (67) A. Denner, S. Dittmaier, S. Kallweit and A. Mück, Comput. Phys. Commun. 195, 161-171 (2015) doi:10.1016/j.cpc.2015.04.021 [arXiv:1412.5390 [hep-ph]].
  * (68) R. V. Harlander, A. Kulesza, V. Theeuwes and T. Zirke, JHEP 11, 082 (2014) doi:10.1007/JHEP11(2014)082 [arXiv:1410.0217 [hep-ph]].
  * (69) N. Kidonakis, PoS DIS2015, 170 (2015) doi:10.22323/1.247.0170 [arXiv:1506.04072 [hep-ph]].
  * (70) C. Muselli, M. Bonvini, S. Forte, S. Marzani and G. Ridolfi, JHEP 08, 076 (2015) doi:10.1007/JHEP08(2015)076 [arXiv:1505.02006 [hep-ph]].
  * (71) Y. B. Shen, R. Y. Zhang, W. G. Ma, X. Z. Li, Y. Zhang and L. Guo, JHEP 10, 186 (2015) [erratum: JHEP 10, 156 (2016)] doi:10.1007/JHEP10(2015)186 [arXiv:1507.03693 [hep-ph]].
  * (72) R. Frederix, D. Pagani and M. Zaro, JHEP 02, 031 (2018) doi:10.1007/JHEP02(2018)031 [arXiv:1711.02116 [hep-ph]].
  * (73) W. Verkerke and D. P. Kirkby, eConf C0303241, MOLT007 (2003) [arXiv:physics/0306116 [physics]].
  * (74) R. Brun and F. Rademakers, Nucl. Instrum. Meth. A 389, 81-86 (1997) doi:10.1016/S0168-9002(97)00048-X
  * (75) S. Ashanujjaman, K. Ghosh and D. Choudhury, In preparation


# Rheology of Pseudomonas fluorescens biofilms: from experiments to DPD
mesoscopic modelling

José Martín-Roca  Dep. Est. de la Materia, Física Térmica y Electrónica,
Universidad Complutense de Madrid, 28040 Madrid, Spain  Valentino Bianco
Departamento de Quimica Fisica, Facultad de Ciencias Quimicas, Universidad
Complutense de Madrid, 28040 Madrid, Spain  Francisco Alarcón  Departamento de
Estructura de la Materia, Física Térmica y Electrónica, Facultad de Ciencias
Físicas, Universidad Complutense de Madrid, 28040, Madrid, Spain  Departamento
de Ingeniería Física, División de Ciencias e Ingenierías, Universidad de
Guanajuato, Loma del Bosque 103, 37150 León, Mexico.  Ajay K. Monnappa
Departamento de Quimica Fisica, Facultad de Ciencias Quimicas, Universidad
Complutense de Madrid, 28040 Madrid, Spain  Paolo Natale  Departamento de
Quimica Fisica, Facultad de Ciencias Quimicas, Universidad Complutense de
Madrid, 28040 Madrid, Spain  Francisco Monroy  Departamento de Quimica Fisica,
Facultad de Ciencias Quimicas, Universidad Complutense de Madrid, 28040
Madrid, Spain  Belen Orgaz  Sección Departamental de Farmacia Galénica y
Tecnología Alimentaria, Facultad de Veterinaria, Universidad Complutense de
Madrid, Madrid, Spain  Ivan López-Montero  Departamento de Quimica Fisica,
Facultad de Ciencias Quimicas, Universidad Complutense de Madrid, 28040
Madrid, Spain  Instituto de Investigación Biomédica Hospital 12 de Octubre
(imas12), Madrid, Spain  Chantal Valeriani  Dep. de Est. de la Materia, Física
Térmica y Electrónica, Universidad Complutense de Madrid, 28040 Madrid, Spain
GISC - Grupo Interdisciplinar de Sistemas Complejos 28040 Madrid, Spain

###### Abstract

The presence of bacterial biofilms in clinical and industrial settings is a
major issue worldwide. A biofilm is a viscoelastic matrix, composed of
bacteria producing a network of Extracellular Polymeric Substances (EPS) to
which bacteria crosslink. Modelling of complex biofilms is relevant to provide
accurate descriptions and predictions that include parameters as
hydrodynamics, dynamics of the bacterial population and solute mass transport.
However, up-to-date numerical modelling, even at a coarse-grained level, is
not satisfactorily. In this work, we present a numerical coarse-grain model of
a bacterial biofilm, consisting of bacteria immersed in an aqueous matrix of a
polymer network, that allows to study rheological properties of a biofilm. We
study its viscoelastic modulus, varying topology and composition (such as the
number of crosslinks between EPS polymers, the number of bacteria and the
amount of solvent), and compare the numerical results with experimental
rheological data of Pseudomonas fluorescens biofilms grown under static and
shaking conditions, as previously described by Jara et al, Frontiers in
Microbiology (2021).

##  I Introduction

Biofilms are synergistic colonies of one or more bacterial species growing on
a solid surface most often in contact with an aqeuous liquid. In contrast to
liquid-floating planktonic bacterial cells, bacteria embedded in solid-
supported biofilms do not have unlimited access to nutrients. The presence of
biofilms does not only provide the physical protection of the embedded cells
against chemical or mechanical challengesPeterson  _et al._ (2015); Flemming
and Wingender (2010), but it also facilitates the evolution of phenotypes that
eventually lead to the resistance to the treatment of antimicrobial agents
Donlan (2001a); Hoiby  _et al._ (2010). Antimicrobial resistance of bacteria
dwelled in biofilms compromises the process of bacterial removal giving rise
to a more persistent bacterial infection. Reason why biofilm formation is a
major health issue Khatoon  _et al._ (2018); Milchev and Binder (2013) and
industrial challenge worldwide Duan  _et al._ (2008); Mansour and Elshafei
(2016); Schultz and Swain (2000). Bacterial infections in plants and animals
lead to food contamination and poisoningZhao  _et al._ (2016) . A biofilm is a
complex and heterogeneous viscoelastic material J.N.Wilking  _et al._ (2011)
produced by bacterial cells. In more structural terms, a biofilm consists of a
network of extracellular polymeric substances (EPS) composed of
polysaccharides, proteins, extracellular DNA, and cells, everything surrounded
by an aqueous environment. EPS serves as a scaffold for matrix cross-linked
bacterial cells embedded within the biofilm. The amount of bacteria with
respect to EPS vary from one tenth to one fourth of the total biofilm
massGarrett  _et al._ (2008); Sveinbjörnsson  _et al._ (2012), depending on
the bacterial species. Mathematical modelling of bacterial biofilms started
more than 30 years ago Klapper and Dockery (2010); Wang and Thanou (2010).
Biofilms have been described via phase-field (continuum) models, in which one
phase represented the biofilm (EPS and bacteria) and the other phase the
surrounding matrix containing the nutrient substrates Zhang  _et al._ (2008).
The tension and compression of a biofilm have been computed using finite
elements simulations to obtain the Young’s modulus of an artificial biofilm
composed of bacteria embedded in an agarose hydrogel Kandemir  _et al._
(2018).

A different approach consisted in implementing an hybrid discrete-continuum
models to incorporate the flow of substrate, bacterial growth or biomass
spreading over the rough biofilm surfaces Picioreanu  _et al._ (1999, 2004).
More recently, mesoscale numerical models such as the immersed boundary-based
method (IBM) or the dissipative particle dynamics (DPD) have been proposed to
study biofilm maturation Horn and Lackner (2014); Eberl  _et al._ (2008); Liu
and Balazs (2018a); Xu  _et al._ (2011a), establishing a connection between
the mechanical stress and strain of the biofilm Bol  _et al._ (2012); Hammond
_et al._ (2014a); Stotsky  _et al._ (2016a). On the one side, to model a
complex system via the IBM model Peskin (1977, 2002) several assumptions need
to be taken into account, such as 1) not considering biofilm growth but
steadiness during the simulation time scale Hammond  _et al._ (2014b), 2)
mimicking biofilm elasticity via linear springs that connect individual
bacterial cells Stotsky  _et al._ (2016b) and 3) capturing the biofilm
viscosity via constitutive equations for stress. As a matter of
phenomenological evidence, macroscopic rheological experiments were considered
on Staphylococcus epidermidis biofilms grown on a plate rheometer Pavlovsky
_et al._ (2013). They were compared qualitatively to numerical results
obtained via the microscale spatial IBM model Stotsky  _et al._ (2018). On the
other side, as refers to dynamic methods Barai  _et al._ (2016), DPD conserve
mass and momentum under steady state conditions. Differently from grid-based
computational methods that do not consider any mechanical constraints, DPD can
be easily implemented when considering objects with complex geometries Xu  _et
al._ (2011b). Like Molecular Dynamics (MD) simulations Mukhi and Vishwanathan
(2021), DPD captures the time evolution of a many-body system with friction
explicitly considered. Unlike MD, DPD allows to model physical phenomena
occurring at longer time scales and length scale Raos and Casalegno (2011).
DPD was thus implemented to simulate biofilm formation on post-coated surfaces
Liu and Balazs (2018b), and the growth of two-dimensional biofilm in flow Xu
_et al._ (2011c). Recently, we have presented a combined numerical and
experimental study on how hydrodynamic stress affect structural features of
Pseudomonas fluorescens biofilms Jara  _et al._ (2021). By means of
rheological measurements on a biofilm grown in the presence and in the absence
of hydrodynamic stress, we demonstrated that bacteria are capable of adapting
to hostile deformations by tailoring the structure of the matrix and its
viscoelastic properties Jara  _et al._ (2021). The results were supported by
numerical results of the rheological properties of a mimetic biofilm simulated
as a first approximation via a DPD model. Thanks to computer simulations, we
learned that biofilms grown under stress were mechanically more stable due to
an increase of the number of crosslinks between EPS polymers.

In this work we probe the rheological features of a P. fluorescens biofilms
varying the frequency of the externally imposed shear flow. To numerically
study the biofilm viscoelastic modulus, we use a coarse-grain DPD model. To
better characterised the biofilm, we vary topology and composition parameters
such as the number of matrix crosslinks between EPS polymers, the number of
bacteria embedded and the amount of solvent. The numerical results are
compared to experimental rheological data of P. fluorescens biofilms grown
under static and shaking conditions.

##  II Numerical and experimental details

In what follows, we present the mesoscopic model used to simulate mature
biofilms, and the experimental details of the viscoelastic measurements of a
P. fluorescens biofilm.

###  II.1 Mesoscopic model for an in silico biofilm

The simulated biofilm consists of three different components: bacteria,
polymers and solvent. The biofilm is prepared in an initial configuration
consisting of a colony of monodisperse bacteria immersed in a periodic box
containing solvent and randomly distributed polymers.

####  II.1.1 Bacteria

Differently from the work described by Raos et al.Guido  _et al._ (2006),
where dispersed filler particles aggregated into a single bacterium, we use
realistic bacteria consisting of already clustered particles. We model the
bacterium as a bacillus (rod-shaped) composed by 440 particles (of type $b$),
each of diameter $\sigma_{b}$, held together via harmonic interactions.

$U_{b}^{\rm bond}=\frac{K_{b}}{2}(r_{i,j}-\sigma_{b})^{2}\qquad,$ (1)

being $r_{i,j}$ the distance between the bonded beads $i$ and $j$,
$\sigma_{b}$ the equilibrium distance, $K_{b}$ the harmonic coupling constant.
The central part of a bacterium is shaped as an empty cylinder of length
$l_{b}/\sigma_{b}=6$, while the extremes of the bacteria are shaped as
spherical caps of radius $R_{ext}/\sigma_{b}=4$. Both the central and extreme
parts have an external radius $R_{ext}/\sigma_{b}=4$ and internal radius
$R_{ext}/\sigma_{b}=2$, i.e. the bacterium membrane is formed by three layers
of particles.

####  II.1.2 Matrix: polymers and solvent

The bacteria are immersed in a bulk matrix made of polymers (EPS) and solvent.
The solvent is represented by $N_{\rm solv}$ particles of type $s$, whereas
the EPS is represented by linear chains of $l$ beads of type $p$ that are
bound via harmonic potential:

$U_{p}^{\rm bond}=\frac{K_{p}}{2}(r_{i,j}-\sigma_{p})^{2}\qquad,$ (2)

where $r_{i,j}$ is the distance between two consecutive beads $i$ and $j$,
$\sigma_{p}$ the equilibrium distance, and $K_{p}$ the coupling constant.

####  II.1.3 Cross-links and non-bonded interactions

The crosslink (CL) interactions, occurring between two beads of type $p$
(polymer-polymer crosslink), or between beads of type $b$ and $p$ (bacterium-
polymer crosslink), are described via harmonic potential:

$U^{\rm CL}=\frac{K}{2}(r_{i,j}-\sigma_{CL})^{2}\qquad.$ (3)

All non-bonded interactions, namely polymer-polymer (pp), polymer-bacterium
(pb), polymer-solvent (ps), bacterium-bacterium (bb), bacterium-solvent (bs),
and solvent-solvent (ss) interactions, are modelled according to the DPD force
field. In particular, the force between two particles $i$ and $j$,
respectively of type $\alpha$ and $\beta$, is expressed as the sum of a
conservative force $\vec{F}^{C}_{i,j}$, a dissipative force
$\vec{F}^{D}_{i,j}$ and a random force $\vec{F}^{R}_{i,j}$ given by

$\displaystyle\vec{F}^{C}_{i,j}=B_{\alpha,\beta}w(r)\hat{r}_{i,j}\qquad
r<r_{c}$ (4)

where $r_{c}$ is the cut-off distance, beyond which all the terms vanish;
$\hat{r}_{i,j}$ is the unit vector along the direction of the vector
$\vec{r}_{i}-\vec{r}_{j}$, indicating the positions of particles $i$ and $j$,
respectively; $B_{\alpha,\beta}$ is the amplitude of the conservative force
between particles of type $\alpha$ and $\beta$;

$\displaystyle\vec{F}^{D}_{i,j}=-\gamma
w^{2}(r)\left(\hat{r}_{i,j}\cdotp\hat{v}_{i,j}\right)\hat{r}_{i,j}$ (5)

where $\vec{v}_{i,j}\equiv\vec{v}_{i}-\vec{v}_{j}$ is the vector difference
between the velocity $\vec{v}_{i}$ of particles $i$, and the velocity
$\vec{v}_{j}$ of particle $j$; $\gamma$ the friction coefficient;

$\displaystyle\vec{F}^{R}_{i,j}=w(r)\left(\dfrac{2k_{B}T\gamma}{dt}\right)^{1/2}\Theta\,\,\,\hat{r}_{i,j}$
(6) $\displaystyle w(r)=1-r/r_{c}$ (7)

Here, $w(r)$ is a weighting factor varying from 0 to 1; $\Theta$ is a Gaussian
random number with zero mean and unit variance; $dt$ the integration time;
$k_{B}$ the Boltzman constant; $T$ the absolute temperature. The sum of the
three DPD forces allows to capturing the long-range correlations induced by
the hydrodynamic interactions, also accounting for thermal fluctuations.

###  II.2 Numerical details

We simulated the biofilm via the LAMMPS open source codePlimpton (1995). The
system evolves according to a velocity Verlet algorithm, keeping fixed the
volume of the simulation box, set to $32\times 32\times 32$ with periodic
boundary conditions. We simulate the biofilm changing the number of bacteria
$N_{b}$ from 90 to 184, the number of polymers $N_{p}$ (from 80 to 240), the
polymer lengths $l$ (from 50 to 200) or the number of solvent particles
$N_{s}$ (from 35000 to 75000). All parameters $N_{p}$, $N_{s}$ and $N_{b}$ are
chosen to guarantee the DPD condition for the densityHoogerbrugge and Koelman
(1992); Groot and Warren (1997): $(N_{p}+N_{b}+N_{s})/(32\times 32\times
32)>3$. Using internal units, we set $r_{c}=1$ for all DPD interactions,
$\sigma_{b}=\sigma_{p}=0.5$, $\gamma=4.5$, $K_{b}=K_{p}=30$. The time step is
set to $dt=0.05$, although we tested our results also against $dt=0.005$.
Moreover, to guarantee that bacteria and polymers are properly hydrated, we
choose the amplitudes of the solvent interactions in eq.4 smaller than the
remaining ones. In particular, we choose $B_{s,s}=B_{s,b}=B_{s,p}=25$ and
$B_{b,b}=B_{p,p}=B_{p,b}=30$. To convert those dimensionless units in real
world units, we consider $k_{B}T=4.11\times 10^{-21}$ J as the characteristic
energy scale in our simulation. Also, we set the system length scale to the
longest dimension of one P. fluorescens bacterium, which can be considered as
approximately $\sim 1.5\mu$m. In our simulation the bacterium is 15 beads
long, which results in a unit length corresponding to $l_{u}\sim 200$ nm
(hence $\sigma_{b}=\sigma_{p}\sim 100$ nm). Moreover, since the mass of a
single bacterium is $\sim 10^{-15}$ Kg, the mass of a single particle is set
to $m_{u}\sim 2.3\times 10^{-18}$ Kg. The time scale is derived accordingly as
$\tau_{\rm intrinsic}=\sqrt{m_{u}l_{u}^{2}/k_{B}T}\sim 5\times 10^{-6}$s. At a
time step of dt=0.01, a period of T=200 corresponds to 2000 time steps that is
within the total lifespan of the simulation runs.

###  II.3 Setting up the initial configuration

We first choose the parameter space to explore. By fixing $N_{b}$, $N_{p}$,
$l$ and $N_{s}$, we equilibrate the system in the absence of polymer cross-
links (CL) for $\sim 8\times 10^{6}$ time steps and perform simulations
changing the parameters within the following ranges: $N_{b}=$ 110, 120, 130,
140, 150, 160, 175, 184, 190, 195; $N_{p}=$ 50, 70, 80, 90, 110, 140, 170;
$N_{s}=$ 20000, 30000, 35000, 40000, 50000, number of CL=100, 300, 500, 700,
900, 1100, 1300, 1500, 1700, 2000, 2300, 2400, 2800, 3000, 3200, 3300, 3600,
3800, 4000, 4300, 4400, 4800, 5300, 6000, 7000.

Having prepared the initial configuration, we need to locate the crosslinks
within the polymer network. To create the network of cross-linked polymers and
bacteria, we randomly place $N_{\rm CL}$ new harmonic bonds 1) between $b$ and
$p$ types of particles, i.e. a bacterium-polymer CL and 2) between two $p$
particles (belonging or not to the same polymer chain), i.e. a polymer-polymer
CL. To avoid crosslinks agglomerating between nearest neighbours, leading to a
clump of $p$ particles of the same polymer chain, we forbid CL between the ten
closest neighbors along the same chain. Therefore, we fix the intermolecular
distance of the polymer loops within one polymer (intra) or between different
polymers (inter) to be at least 10 $\sigma$ apart. To prevent the formation of
too many cross-links per particle which would result in a globule-like
clusters of beads, we assume that any particle of type $b$ can form at most
one CL, while any $p$ particle can form at most two CL. For any choice of
$N_{b}$, $N_{p}$, we create an initial configuration and form a number of CL
ranging from 2300 up to 5400. After the formation of the CLs, a second
equilibration run of $\sim 10^{6}$ time steps is performed to relax the
polymer-bacteria network. A characteristic snapshot of the biofilm is shown in
Fig. 1.

Figure 1: Snapshots of the simulation box containing 184 bacteria (in green),
80 polymers (in yellow) and 35000 solvent particles (in red). The size of a
single particle has been chosen for the sake of visualization and does not
reflect the interaction radius.

###  II.4 Simulated rheology: numerical details

To study the rheological properties of the biofilm, we compute the shear
modulus by monitoring the shear stress response $\sigma(t)$ under a sinusoidal
deformation. We apply an external oscillatory shear deformation with frequency
$\nu$ along th $X$–$Y$ plane by changing the box size $L_{z}$Raos and
Casalegno (2011); Jara  _et al._ (2021) according to

$L_{z}(t)=L_{0}+A\sin(2\pi\nu t)\qquad,$ (8)

being $L_{0}=32$ the initial box size and $A$ the amplitude of the
oscillation. Under shearing conditions, the particle velocities are remapped
every time they cross periodic boundaries. Within our simulations, we have
explored the following values for the amplitude $A$ :A= 2.5, 3.2, 4.8, 8, 12,
16, 20, 24, 28, 32, 36; and the period $T$: T$\equiv\frac{2\pi}{\nu}$= 20, 50,
100, 200, 500, 1200, 2000, 2500, 5000, 10000, 20000, 100000, 200000.

According to Raos et al. Guido  _et al._ (2006), the resulting stress has been
calculated by fitting the $xy$ component of the shear stress with:

$\sigma_{xy}(t)=\sigma^{\prime}\cdot\sin(2\pi\nu
t)+\sigma^{\prime\prime}\cdot\cos(2\pi\nu t)$ (9)

and the amplitude values of $\sigma^{\prime}$ and $\sigma^{\prime\prime}$ were
used to compute the the in-phase and out-of-phase components of the complex
shear modulus

$G^{\prime}=\frac{\sigma^{\prime}}{(A/L_{0})}\qquad
G^{\prime\prime}=\frac{\sigma^{\prime\prime}}{(A/L_{0})}$ (10)

being $G^{\prime}$ the storage modulus and $G^{\prime\prime}$ the loss
modulus, respectively. Rheological properties of solid biomaterials are
characterized by a finite storage modulus $(G^{{}^{\prime}}>0$). Whereas fluid
viscous materials are characterized by near zero storage modulus
$(G^{{}^{\prime}}\sim 0$) and a large loss modulus
$(G^{{}^{\prime\prime}}>0$)Donlan (2001b); Alonci  _et al._ (2018).

From a phenomenological point of view, the computed stress-strain curves are
mainly characterized by two regimes a) The stress-strain curves are linear at
low strain, where the stress is proportional to the the strain. b) Beyond a
yield strain, where the stress-strain curves deviate from the linear regime.
Both regimes crossover at the yield point, which fixes the yield strain in
which the system softens becoming effectively plastic.

For the considered frequency ranging from a value of $3\times 10^{-5}$ to a
value of $0.3$ (in internal units), we run the simulations completing at least
3 box-shearing oscillations. We consider simulations amplitudes $A$ (ranging
from 4 to 36), corresponding to a deformation of up to the 112% of the box
size. However, we will mainly focus on data obtained for $A=24$ (75%), since:
the stress-strain curve is linear until around a deformation of 50% and starts
being non linear for a larger deformation; when varying all parameters
($N_{b}$, $N_{p}$, …) while the behaviour of $G^{{}^{\prime}}$ and
$G^{{}^{\prime\prime}}$ is not strongly affected by the choice of A, their
response is stronger and the noise lower for $A=75\%$. Thus, we have
considered a period $T$ ranging from $2\times 10$ to $2\times 10^{5}$.

###  II.5 Experimental details: Microorganisms and growth conditions

The strain P. fluorescens B52, originally isolated from raw milk Richardson
and Te Whaiti (1978), was used as model microorganism. Overnight precultures
and cultures were incubated at 20∘C under continuous orbital shaking (80 rpm)
in tubes containing 10 ml Trypticase Soy Broth (TSB, Oxoid). Cells were
recovered by centrifugation at 4,000 x  g for 10 min (Rotor SA-600; Sorvall
RC-5B-Refrigerated Superspeed Centrifuge, DuPont Instruments) and washed twice
with sterile medium. Cellular suspensions at OD600 were first adjusted to 0.12
(equivalent to 108 cfu/ml) and then diluted to start the experiments at 104
cfu/ml.

###  II.6 Experimental details: Biofilm development

Biofilms were developed using borosilicate glass surfaces ($20\times 20$ cm)
as adhesion substrate. Five glass plates were held vertically into the
sections of a tempered glass separating chamber provided with a lid. The whole
system was heat-sterilized as a unit before aseptically introducing 2 mL of
the inoculated culture medium. To check the effect of hydrodynamic stress on
biofilm mechanical properties, incubation was carried out for 96 h at 20∘C
both in an orbital shaker at 80 rpm (shaken sample) and statically (static
sample). For biofilm recovery, plates were aseptically withdrawn, rinsed with
sterile saline in order to eliminate weakly attached cells and then scraped to
remove the attached biomass (cells + matrix) from both faces of the plates.
For rheology measurements, biofilm material was casted every 24 h to be
directly poured onto the rheometer plate. Experiments were run in triplicate.

###  II.7 Experimental rheology

As previously described Jara  _et al._ (2021), the viscoelastic response was
determined in a hybrid rheometer under oscillatory shear stress-control
(Discovery HR-2, TA instruments). We used a cone-plate geometry (40 mm
diameter) and a Peltier element to control the temperature through an external
an external water thermostat. Measurements were performed at a 1 mm gap
between the Peltier element and the plate tool (TA instruments), where a
sinusoidal strain $\gamma$ of amplitude $\gamma_{0}$ was performed at a
frequency $\omega$, as $\gamma(t)=\gamma_{0}sin(\omega t)$. The shear stress
exerted on the biofilm was monitored as $\sigma(t)=G^{*}\gamma(t)$, where
$G^{*}$ is the viscoelastic modulus
$G^{*}=G^{{}^{\prime}}+iG^{{}^{\prime\prime}}$, where $G^{{}^{\prime}}$ is the
storage modulus accounting for shear rigidity, $G^{{}^{\prime\prime}}$ is the
loss modulus accounting for viscous friction. These properties were considered
under the same definitions as in Eq. (10).

##  III Results

###  III.1 Matching the parameters of the biofilm model

To start with, we study the stress response of the biofilm containing cross-
linked polymer loops regardless of their intermolecular distances along the
polymer chain. We analyze how parameters affect the in-phase
($G^{{}^{\prime}}$) stress response function. In Fig. 2 we represent the
storage modulus $G^{{}^{\prime}}$ as a function of the shear amplitude,
keeping fixed $N_{b}=150$, $N_{p}=80$, $T=2\times 10^{5}$, $l=100$ and the
total number of $CL=1600$ cross-links ($N_{CL,pp}+N_{CL,pb}$).

Figure 2: Storage Modulus as a function of the amplitude for different cross-
link topologies in the presence (filled squares) and absence (empty circles)
of polymer loops. The color code indicates different numbers of polymer-
polymer (pp), polymer-bacteria (pb) cross-links: $N_{CL,pp}$ =200/
$N_{CL,pb}$=1400 (blue); $N_{CL,pp}$ =400/ $N_{CL,pb}$=1200 (red) and
$N_{CL,pp}$ =600/ $N_{CL,pb}$=1000 (green). All data refer to $N_{b}=150$,
$N_{p}=80$, $N_{s}=45000$ and $T=2\times 10^{5}$.

Figure 3: Stress components $G^{\prime}$ and $G^{\prime\prime}$ for $A=75\%$,
$CL=2300$ and $T=200$ as function of (a,d) the number of bacteria $N_{b}$;
(b,e) the number of polymers $N_{p}$ and solvent (c,f).

We consider two systems with different cross-link topologies: one with polymer
cross-links (filled squares) and one without polymer cross-links (empty
circles). The presence of polymer loops between $p$ particles belonging to the
same polymer chain affects the storage modulus. Since $G^{\prime}$ of the
cross-linked biofilm (filled squares, in Fig.2) is lower than the one of the
cross-linked biofilm that does not include polymer loops (Fig.2, empty
symbols).

This finding indicates the importance to include polymer loops with a biofilm
model. As expected, independently of the presence (or absence) of polymer
loops, the storage modulus increases with the strain-amplitude. In the
presence of polymer loops (for the chosen fixed parameters), we observe that
the shear response is not affected by the ratio between polymer-polymer ($pp$)
versus polymer-bacteria ($pb$) cross-links. On the other hand, when polymer
loops are not present within the same polymer chain, the system response
increases (empty symbols in Fig. 2). We think that this increase in system
response is due to the presence of a larger number of interactions
(connections) between bacteria and matrix polymers. This leads to a a denser
and stiffer biofilm and as a consequence generates a larger response function.

###  III.2 The shear stress dependence on $N_{b}$, $N_{p}$ and $N_{s}$

To check whether the shear stress depends on the number of bacterial cells
($N_{b}$), in Fig. 3, panel a)we study the $G^{{}^{\prime}}$ (in-phase) and
$G^{{}^{\prime\prime}}$ (out-of-phase) response as a function of $N_{b}$
(fixing all other parameters). We use an intermediate number of cross-links
($CL=2300$), and intermediate values $N_{p}=80$, $l_{p}=100$, $N_{s}=35000$
within the explored range. Moreover, we apply a deformation amplitude of
$A=75$ $\%$ and a period of the deformation of $T=200$. Fixing the different
oscillation’s periods to $T=2\times 10^{2}$, we observe a relevant response of
$G^{{}^{\prime\prime}}$. Whereas at lower frequencies we have found an
amplitude ($A$) depending limiting value of $G^{{}^{\prime\prime}}$ (as we
will also observe later in Fig 4).

Our results show that the in-phase response (accounting for the storage, or
rigidity modulus $G^{{}^{\prime}}$) remains constant through of the $N_{b}$
range (Fig. 3, panel a). This result is in agreement with previously reported
data J.N.Wilking  _et al._ (2011). Considering the bacterial biofilm as a
polymer-colloid viscoelastic medium J.N.Wilking  _et al._ (2011),
$G^{{}^{\prime}}$ should increase with the bacterial density until reaching a
solid-like plateau. Beyond this point, the system behaves as a viscous
(liquid-like) fluid. Our findings point out that the probed $N_{b}$-range
falls into the elastic response regime, where $G^{{}^{\prime}}$ is not
affected by $N_{b}$.

The loss modulus $G^{{}^{\prime\prime}}$ increases linearly with $N_{b}$
(Fig.3, panel d), reflecting the high friction imposed on the system by
dragging the large bacteria while increasing their density. The viscous loss
almost reaches the values of the storage modulus for the largest number of
bacteria considered. However, since $G^{{}^{\prime\prime}}$ is always smaller
than $G^{{}^{\prime}}$, the biofilm system behaves as a pasty material, with a
relative high structural rigidity $G^{{}^{\prime}}>G^{{}^{\prime\prime}}$
instead of a viscous fluid where $G^{{}^{\prime}}=0$ and
$G^{{}^{\prime\prime}}>0$. Thus our model shows a rheological behavior of a
quite resilient material with a relatively high mechanical compliance.

To visualize the effect of an increasing number of entangled polymers on the
shear response, we consider $N_{b}=150$, $l_{p}=100$, $N_{s}=35000$, $CL=2300$
with an amplitude of $A=75$ $\%$ and period of the deformation of $T=200$ as
boundary conditions for our simulations (Fig.3, panel b and e). When we
increase the number of polymers and maintain a constant bacterial density, we
observe that it is easier for the biofilm to generate a solid-like network
with an increased rigidity (Fig.3b). This suggests that bacteria colonize a
preexisting polymer template and modify it according to their necessities. The
out-phase stress response increases with increasing number of polymers (Fig.3
e) or solvent (Fig.3 f) present in the biofilm. Since this the increase of
$G^{{}^{\prime\prime}}$ is quite small when compared to the response caused by
the presence of bacteria (Fig.3 d), we consider $G^{{}^{\prime\prime}}$ for
polymer and solvent to be constant.

When studying the mechanical response of the system for different hydration
levels ($CL=2300$, $N_{b}=184$, $N_{p}=80$, $l_{p}=100$ with an amplitude of
$A=75$ % and deformation period of $T=200$), we observe that with increasing
number of solvent particles, the biofilm is fully hydrated and the storage
modulus decreases (Fig. 3c). This is rather expected as under these
experimental conditions the average distance between bacteria and polymers
increases and the network becomes mechanically fragile (shear softening). As
clearly observed from the data, the biofilm model is able to absorb the
solvent. The more solvent is present, the softer the storage modulus (Fig 3c).
On the other side, the hydration level seems not to affect
$G^{{}^{\prime\prime}}$, which only slightly increases. This softening is not
reflected in the loss modulus that remains more or less constant over the
tested frequencies (Fig 3d). This implies that the present cross-linked matrix
polymers might be entangeled, which allows expansion of the biofilm but at the
same time provides rigidity.

###  III.3 Shear stress dependence on the shearing frequency: simulations
versus experiments

Finally, we investigate the stress response as a function of the frequency of
the sinusoidal shearing strain $\omega$. On the one side, we perform a
rheology experiment on a P. fluorescens biofilm. We grow two P. fluorescens
biofilms for 24 hours: one under static conditions, and the second one under
shaking conditions. Next, we locate the biofilm ex vivo in the rheometer and
measure $G^{{}^{\prime}}$ (filled symbols in Fig.4, top panel) and
$G^{{}^{\prime\prime}}$ (empty symbols in Fig.4, top panel) for the static
biofilm (Fig.4, top panel, red curves) and for the shaking biofilm (Fig.4, top
panel, black curves). On the other side, we perform a numerical rheology
experiment on the model biofilm described earlier. After having studied the
parameter space in previous sections, we choose the following parameters:
$CL=2000$, $N_{b}=184$, $N_{p}=80$, $l_{p}=100$, $N_{s}=35000$ and $A=75$
$\%$. We compute $G^{{}^{\prime}}$ (filled symbols in Fig.4, bottom panel) and
$G^{{}^{\prime\prime}}$ (empty symbols in Fig.4, bottom panel) as a function
of the frequency of the externally imposed shear, and compare the experimental
wet-lab data with the biofilm simulations. P. fluorescens biofilm grown after
24 hours under shaking conditions (in black) present higher shear moduli (both
$G^{{}^{\prime}}$ and $G^{{}^{\prime\prime}}$) than those measured for P.
fluorescens biofilm grown under static conditions, independently on the
frequency. Qualitatively, we observe this behavior for any P. fluorescens
biofilms grown longer than 24h (data not shown).

Figure 4: Top) Experimental measurements of the elastic modulus (filled
symbols) and loss modulus (empty symbols) for two different biofilms: one
grown under shaking conditions (black) and one under static conditions (red).
Bottom) $G^{{}^{\prime}}$ (filled symbols) and $G^{{}^{\prime\prime}}$ (empty
symbols) computed for the biofilm model. Data refer to $N_{b}=184$,
$N_{p}=80$, $l=100$, $N_{s}=35000$ and 2000 cross-links.

By increasing the externally imposed stress (or frequency), we observed
increasingly larger absolute values of the $G^{{}^{\prime\prime}}$, which
become comparable to $G^{{}^{\prime}}$ for large frequencies (at least under
static conditions, in red in Fig. 4), as one would expect in a viscoelastic
material. In the biofilm grown under shaking conditions
$G^{{}^{\prime\prime}}$ does not reach the same value of $G^{{}^{\prime}}$ for
large frequencies (in black in Fig. 4). Due to the technical limitations of
the used rheometer, we underline that we could not measure the stress moduli
of the P. fluorescens biofilm grown under shaking conditions for frequencies
beyond 100 Hz (Fig. 4, top panel, black symbols). However, we expect that
eventually the trend of the stress response of the P. fluorescens biofilms
developed under shaken conditions will behave similarly to the observed stress
response of the biofilm developed under static conditions (Fig. 4, top panel).

The obtained values for the storage moduli (filled symbols) are always larger
than the ones for the loss moduli (empty symbols), both in experiments and in
simulations. This result is consistent with our already published workJara
_et al._ (2021), where we showed that P. fluorescens biofilms grow better
under shaking growth conditions than under static ones, the former leading to
a higher bacterial density and a thicker biofilm matrix. Augmented aeration
and nutrient concentration within the biofilm provides an optimal growth
conditions for the bacteria that eventually leads to the production of a
larger amount of cross-links within the EPS polymer matrix and produces a
stiffer biofilm structure Jara  _et al._ (2021).

According to the established experimental conditions, we numerically compute
the shear moduli of the biofilm model by changing the stress frequency over 4
order of magnitudes (Fig. 4, bottom panel) . We assume that the main
difference between the biofilm grown under static and shaking conditions is
only in in the matrix composition, being the one of the shaken-grown biofilm
richer in cross-links. In general terms, the numerical results (expressed in
internal units, i.u.) exhibits the same overall trend (Fig. 4, top panel)
observed for the P. fluorescens biofilms (both static and shaken). Overall,
for all tested frequencies, $G^{{}^{\prime}}$ is larger than
$G^{{}^{\prime\prime}}$ that remains almost constant. However,
$G^{{}^{\prime\prime}}$ increases with increasing frequency and eventually
converges with $G^{{}^{\prime}}$ at the highest simulated frequency. A
qualitative difference is the drop in G” shown by simulations at low
frequencies. A dominance of the flowing solvent is the halmark for these
simulations at low frequency of shearing. This is not however observed in
experiments (see Fig. 4; top panel), which a change of the number of cross-
links $CL$, bacteria $N_{b}$, polymers $N_{p}$ or hydration level $N_{s}$
within the studied parameter space does not change the trend of our obtained
computational results. A time-modulated solvent-matrix interaction could be
introduced to correct this viscoelastic bias of simulations not found in
experiments.

##  IV Discussion

Biofilms are composed of bacteria that secrete a network of extracellular
polymeric substances (EPS) producing a viscoelastic matrix to which they
eventually crosslink. Within a clinical or industrial environment, biofilm
formation faces major issues. The mathematical modelling of complex biofilms
might be of support, providing accurate descriptions and predictions that
include physical parameters as hydrodynamics, solute mass transport or
dynamics of the bacterial population within the biofilm. However, up-to-date
modelling of biofilms, even at a coarse-grained level, has not been completely
satisfactory. The parametrization of the individual components of a biofilm
increases the complexity of the algorithm beyond computational capacity Calude
and Longo (2017). Therefore, a simplification is needed that often leads to
covering relevant experimental information.

In our work, we present a mesoscopic coarse grained dissipative particle
dynamics (DPD) approach to model a mature bacterial biofilm. DPD models allow
to simulate short range interactions of any particle and thus give a detailed
picture of its close surroundings allowing an increase in spatial resolution.
This is an important advantage over generally used long-range interaction
modelsGroot (2004). Moreover the "soft" nature of DPD interactions allows to
increase the integration time by few orders of magnitude, which permits to
explore time scales normally inaccessible to atomistic simulations. Most
importantly, our model explicitly includes hydrodynamics and thermal
fluctuations, which are known to play an important role in biological systems
such as biofilms Krsmanovic  _et al._ (2021).

We generated a series of biofilm conformations by changing the composition,
the number of bacteria, the number of polymers or varied the cross-linking and
the hydration level of the biofilm. Next, the biofilm underwent deformations
when exposed to an external oscillatory shear with different amplitudes and
frequencies spanning 4 order of magnitudes. We computed the shear moduli and
presented the model predictions for the different systems finding the
computational results in agreement with experiments. When varying the number
of bacteria, we observed that $G^{{}^{\prime}}$ remains constant throughout
the $N_{b}$ range, in agreement with considering a biofilm as a polymer-
colloid viscoelastic medium. The loss modulus $G^{{}^{\prime\prime}}$
increased linearly with $N_{b}$ reaching the values of the storage modulus for
the largest number of bacteria. However, since $G^{{}^{\prime\prime}}$ is
always smaller than $G^{{}^{\prime}}$, the biofilm system behaves as a pasty
material, with a relative high structural rigidity
$G^{{}^{\prime}}>G^{{}^{\prime\prime}}$. When increasing the number of
entangled polymers, we observed that it was easier for the biofilm to generate
a solid-like network with an increased rigidity. The out-phase stress response
increases with increasing number of polymers. However, since the increase of
$G^{{}^{\prime\prime}}$ is quite small when compared to the response caused by
the presence of bacteria, we considered $G^{{}^{\prime\prime}}$ for polymer
and solvent to be constant. Therefore, the present cross-linked matrix
polymers might be entangeled, which allows expansion of the biofilm but at the
same time provides rigidity. When studying the mechanical response of the
system for different hydration levels, we observed that increasing the number
of solvent particles, the biofilm was fully hydrated and the storage modulus
decreased. The more the solvent, the softer the storage modulus.

We have chosen a set of parameter to compare our numerical data with
rheological experiments performed on P. fluorescens biofilms grown under
static and shaken conditions. We choose to study biofilms formed by P.
fluorescens for two reasons. 1) It is the same model system we chose in our
previous studyJara  _et al._ (2021), and it will be easier for us to compare
these results to the already published ones. 2)Although the model
microorganism used in this work for experimental data is considered non-
pathogenic, other members of the genus P. fluorescens such as Pseudomonas
aeruginosa is an opportunistic pathogen that is frequently associated with
chronic biofilm infections. Overall, the composition of the EPS is similar in
both species, with a high proportion of acetylated polysaccharides (alginate-
like) and extracellular DNA (Kives  _et al._ , 2006; Jennings  _et al._ , )
and therefore suggested to expose a similar mechanical behaviour. Knowledge on
the timing of mechanical transitions within biofilms may guide future
strategies that allow the penetration of antimicrobial agents into “soften”
biofilm matrices.

By accurately comparing the experimental to the numerical results, we conclude
that our proposed coarse grained model catches qualitatively the moduli
behaviors over four order of magnitudes as observed in Fig 4. The measured
elastic modulus is always higher in experiments and simulations than the loss
modulus. Although $G^{{}^{\prime}}$ and $G^{{}^{\prime\prime}}$ tend to
converge for high frequency, the simulations contain more viscous relaxation
than reflected by experiments. Overall, the qualitative trends observed as
soft solids for both studied biofilms (static and shaken) are recovered by our
numerical findings, which could be eventually enhanced by including sticking
details on the solvent-matrix interactions. Moreover, this concordance
exhibits the potential of the use of a coarse graianed DPD approach when
modelling complex biofilms. Even though the model strongly relies on an a
priori detailed study of the parameter space, in future prospective, this will
allow us to evaluate in more detail the biofilm transitions from a
predominantly solid-like system to predominantly liquid-like system, through
the interplay of the elastic and loss moduli upon changing the biofilm
composition.

## Conflicts of interest

There are no conflicts to declare.

## Acknowledgments

The authors acknowledge funding from Grant PID2019-105606RB-I00,
FIS2016-78847-P, PID2019-108391RB-I00 (to FM), and FIS2015-70339 (to FM and
ILM) ID2019-105343GB-I00 of the MINECO and the UCM/Santander PR26/16-10B-2.
Francisco Alarcón acknowledges support from the “Juan de la Cierva” program
(FJCI-2017-33580). AKM is recipient of a Sara Borrell fellowship (CD18/00206)
financed by the Spanish Ministry of Health. V.B. acknowledges the support from
the European Commission through the Marie Skłodowska-Curie Fellowship No.
748170 ProFrost. The authors acknowledge the computer resources from the Red
Española de Supercomputacion (RES) FI-2020-1-0015 and FI-2020-2-0032, and from
the Vienna Scientific Cluster (VSC).

## References

  * Peterson  _et al._ (2015) B. W. Peterson, Y. He, Y. Ren, A. Zerdoum, M. R. Libera, P. K. Sharma, A.-J. van Winkelhoff, D. Neut, P. Stoodley, H. C. van der Mei, and H. J. Busscher, FEMS Microbiology Reviews 39, 234 (2015).
  * Flemming and Wingender (2010) H. Flemming and J. Wingender, J. Nat Rev Microbiol 8, 623 (2010).
  * Donlan (2001a) R. Donlan, J. Emerg. Infect. Dis. 7, 277 (2001a).
  * Hoiby  _et al._ (2010) N. Hoiby _et al._ , Int J Antimicrobial Agents 35, 322 (2010).
  * Khatoon  _et al._ (2018) Z. Khatoon, C. D. McTiernan, E. J. Suuronen, T.-F. Mah, and E. I. Alarcon, Heliyon 4, e01067 (2018).
  * Milchev and Binder (2013) A. Milchev and K. Binder, Macromolecules , 131031083721001 (2013).
  * Duan  _et al._ (2008) J. Duan, S. Wu, X. Zhang, G. Huang, M. Du, and B. Hou, Electrochimica Acta 54, 22 (2008).
  * Mansour and Elshafei (2016) R. Mansour and A. Elshafei, British Biotechnology Journal 14, 1 (2016).
  * Schultz and Swain (2000) M. P. Schultz and G. W. Swain, Biofouling 15, 129 (2000).
  * Zhao  _et al._ (2016) G. Zhao, Y. J. Yu, J. L. Yan, M. C. Ding, X. G. Zhao, and H. Y. Wang, Physical Review B 93, 140203 (2016).
  * J.N.Wilking  _et al._ (2011) J.N.Wilking, T.E.Angelini, A.Seminara, M.P.Brenner, and D.A.Weitz, MRS Bulletin 36, 385 (2011).
  * Garrett  _et al._ (2008) T. R. Garrett, M. Bhakoo, and Z. Zhang, Progress in natural science 18, 1049 (2008).
  * Sveinbjörnsson  _et al._ (2012) B. R. Sveinbjörnsson, R. A. Weitekamp, G. M. Miyake, and Y. Xia, Proceedings Of The National Academy Of Sciences Of The United States Of America 109, 14332 (2012).
  * Klapper and Dockery (2010) I. Klapper and J. Dockery, SIAM Rev. 52, 221 (2010).
  * Wang and Thanou (2010) M. Wang and M. Thanou, Pharmacological research : the official journal of the Italian Pharmacological Society 62, 90 (2010).
  * Zhang  _et al._ (2008) T. Zhang, N. G. Cogan, and Q. Wang, SIAM J. Appl. Math 69, 641 (2008).
  * Kandemir  _et al._ (2018) N. Kandemir, W. Vollmer, N. S. Jakubovics, and J. Chen, Scientific Reports 8, 10893 (2018).
  * Picioreanu  _et al._ (1999) C. Picioreanu, M. C. M. van Loosdrecht, and J. J. Heijnen, Water Sci. Technol. 39, 115 (1999).
  * Picioreanu  _et al._ (2004) C. Picioreanu, J. U. Kreft, and M. C. M. van Loosdrecht, Appl. Environ. Microbiol. 70, 3024 (2004).
  * Horn and Lackner (2014) H. Horn and S. Lackner, Adv Biochem Eng Biotechnol 146, 53 (2014).
  * Eberl  _et al._ (2008) H. J. Eberl _et al._ ,  _Computing Intensive Simulations in Biofilm Modeling_ (22nd International Symposium on High Performance Computing Systems and Applications, 2008).
  * Liu and Balazs (2018a) Y. Liu and A. C. Balazs, Langmuir 34, 1807 (2018a).
  * Xu  _et al._ (2011a) Z. Xu _et al._ , Phys Rev E 83, 066702 (2011a).
  * Bol  _et al._ (2012) M. Bol _et al._ , Critical Reviews in Biotechnology 33, 145 (2012).
  * Hammond  _et al._ (2014a) J. Hammond _et al._ , Computer Modeling in Engineering and Sciences 98, 295 (2014a).
  * Stotsky  _et al._ (2016a) J. A. Stotsky _et al._ , J. Comp. Phys. 317, 204 (2016a).
  * Peskin (1977) C. Peskin, J. Comput.Phys. 81, 372 (1977).
  * Peskin (2002) C. Peskin, Acta Numer. 11, 479 (2002).
  * Hammond  _et al._ (2014b) J. Hammond, E. Stewart, J. Younger, M. Solomon, and D. Bortz, Comput. Model. Eng. Sci. 98, 295 (2014b).
  * Stotsky  _et al._ (2016b) J. A. Stotsky, J. F. Hammond, L. Pavlovsky, E. J. Stewart, J. G. Younger, M. J. Solomon, and D. M. Bortz, J. Comput. Phys 317, 204 (2016b).
  * Pavlovsky  _et al._ (2013) L. Pavlovsky, J. G. Younger, and M. J. Solomon, Soft matter 9, 122 (2013).
  * Stotsky  _et al._ (2018) J. A. Stotsky, V. Dukic, and D. M. Bortz, Euro. Jnl of Applied Mathematics , 1 (2018).
  * Barai  _et al._ (2016) P. Barai, A. Kumar, and P. P. Mukherjee, PLoS One 11, e0165593 (2016).
  * Xu  _et al._ (2011b) Z. Xu, P. Meakin, A. Tartakovsky, and T. D. Scheibe, Physical Review E 83, 066702 (2011b).
  * Mukhi and Vishwanathan (2021) M. Mukhi and A. Vishwanathan, Molecular Diversity , 1 (2021).
  * Raos and Casalegno (2011) G. Raos and M. Casalegno, J. Chem. Phys. 134, 054902 (2011).
  * Liu and Balazs (2018b) Y. Liu and A. C. Balazs, Langmuir 34, 1807 (2018b).
  * Xu  _et al._ (2011c) Z. Xu, P. Meakin, A. Tartakovsky, and T. D. Scheibe, Physical Review E 83, 066702 (2011c).
  * Jara  _et al._ (2021) J. Jara, F. Alarcon, A. K. Monnappa, J. I. Santos, V. Bianco, P. Nie, M. P. Ciamarra, A. Canales, L. Dinis, I. Lopez-Montero, C. Valeriani, and B. Orgaz, Frontiers in Microbiology 11, 588884 (2021).
  * Guido  _et al._ (2006) R. Guido, M. Margherita, and S. Elli, Macromolecules 39, 6744 (2006).
  * Plimpton (1995) S. Plimpton, Journal of computational physics 117, 1 (1995).
  * Hoogerbrugge and Koelman (1992) P. Hoogerbrugge and J. M. V. A. Koelman, EPL 19, 155 (1992).
  * Groot and Warren (1997) R. D. Groot and P. B. Warren, The Journal of chemical physics 107, 4423 (1997).
  * Donlan (2001b) R. M. Donlan, Emerging infectious diseases 7, 277 (2001b).
  * Alonci  _et al._ (2018) G. Alonci, F. Fiorini, P. Riva, F. Monroy, I. López-Montero, S. Perretta, and L. De Cola, ACS Applied Bio Materials 1, 1301 (2018).
  * Richardson and Te Whaiti (1978) B. C. Richardson and I. E. Te Whaiti, New Zealand Journal of Dairy Science and Technology 13, 172 (1978).
  * Calude and Longo (2017) C. S. Calude and G. Longo, Foundations of science 22, 595 (2017).
  * Groot (2004) R. D. Groot, in  _Novel Methods in Soft Matter Simulations_ (Springer, 2004) pp. 5–38.
  * Krsmanovic  _et al._ (2021) M. Krsmanovic, D. Biswas, H. Ali, A. Kumar, R. Ghosh, and A. K. Dickerson, Advances in Colloid and Interface Science 288, 102336 (2021).
  * Kives  _et al._ (2006) J. Kives, B. Orgaz, and C. Sanjosé, Colloids Surf B Biointerfaces 52, 123 (2006).
  * (51) L. K. Jennings, K. M. Storek, H. E. Ledvina, C. Coulon, L. S. Marmont, I. Sadovskaya, P. R. Secor, B. S. Tseng, M. Scian, A. Filloux, D. J. Wozniak, P. L. Howell, and M. R. Parsek, Proceedings of the National Academy of Sciences 112, 11353.