t1t1footnotetext: This article is registered under preprint number: CERN-
TH-2016-047, DESY 16-042, UCI-HEP-TR-2016-05

# On Mono-$W$ Signatures in Spin-$1$ Simplified Models

Ulrich Haisch<EMAIL_ADDRESS>Felix Kahlhoefer
<EMAIL_ADDRESS>and Tim M.P. Tait<EMAIL_ADDRESS>Rudolf Peierls
Centre for Theoretical Physics, University of Oxford,
1 Keble Road, Oxford OX1 3NP, United Kingdom  CERN, Theory Division, CH-1211
Geneva 23, Switzerland  DESY, Notkestraße 85, D-22607 Hamburg, Germany
Department of Physics and Astronomy, University of California, Irvine,
California 92697, USA

###### Abstract

The potential sensitivity to isospin-breaking effects makes LHC searches for
mono-$W$ signatures promising probes of the coupling structure between the
Standard Model and dark matter. It has been shown, however, that the strong
sensitivity of the mono-$W$ channel to the relative magnitude and sign of the
up-type and down-type quark couplings to dark matter is an artefact of
unitarity violation. We provide three different solutions to this mono-$W$
problem in the context of spin-$1$ simplified models and briefly discuss the
impact that our findings have on the prospects of mono-$W$ searches at future
LHC runs.

###### keywords:

Mostly Weak Interactions: Beyond Standard Model

##  1 Introduction

At least two reasons make the process $pp\rightarrow W+E_{T,\rm miss}$
interesting in the context of dark matter (DM) searches at the LHC. First, for
leptonically decaying $W$ bosons, this process yields a distinct signature
that can be searched for in dedicated analyses, which suffer from
significantly less background than for instance mono-jet searches. Second,
hadronic $W$-boson decays will instead lead to a ${\rm jets}+E_{T,\rm miss}$
final state and therefore render an electroweak (EW) contribution to the mono-
jet channel.111Note that in spite of the name, mono-jet searches do not
usually veto events with two jets. Indeed, even higher jet multiplicities are
included in the latest LHC searches.

These so-called mono-$W$ searches have received significant interest because
they are potentially sensitive to the sign between the up-type and down-type
quark couplings to DM [1]. The reason is that two different diagrams
contribute to this process and therefore interference effects can be
important. Indeed, LHC Run I analyses [2, 3, 4] based on an effective field
theory (EFT) approach to parameterise the interactions of DM [5, 6, 7] found a
striking difference in the predicted fiducial cross sections between same-sign
(SS) and opposite-sign (OS) couplings. In the OS case, leading to constructive
interference between up-type and down-type quark contributions, the mono-$W$
results in fact turn out to set the strongest limits on the suppression scale
of the unknown mediating interaction, surpassing the EFT limits that arise
from mono-jet searches.222The ATLAS collaboration has very recently also
searched for a mono-$W$ signal in $13\,{\rm TeV}$ data [8], upgrading their
latest $8\,{\rm TeV}$ analysis [9]. In direct detection experiments, on the
other hand, OS couplings lead to destructive interference and correspondingly
smaller event rates, making the LHC a particularly promising probe for this
scenario.

Figure 1:  Prototypes of Feynman diagrams that lead to a mono-$W$ (left) and a
mono-jet signal (middle and right). While the graph in the middle exemplifies
the QCD contribution to this process, the diagram on the right-hand side
represents an EW correction.

Subsequently, it has however been pointed out in [10] that EFT interpretations
of mono-$W$ signals have issues with unitarity violation, if the effective
higher-dimensional operators break the $SU(2)_{L}$ invariance of the Standard
Model (SM). This is for example the case when considering pure vector DM-quark
operators with OS Wilson coefficients for up-type and down-type quarks, as
done in [2, 3, 4]. The enhancements of the mono-$W$ cross section found in
these analyses is thus spurious, because it is due to the emission of
longitudinal $W$ bosons. The EFT study [10] has recently been extended in [11]
to the case of simplified models with $t$-channel exchange of a new coloured
scalar. In particular, it has been shown that if these theories are formulated
before EW symmetry breaking (EWSB) in an $SU(2)_{L}$ gauge-invariant way,
effects associated with isospin breaking are generically small. This leads to
the conclusion that isospin-violating DM-quark couplings are not expected to
increase the sensitivity of mono-$W$ searches in $t$-channel simplified
models. It was also pointed out in [11] that a similar solution is expected to
be present in $s$-channel simplified models.

The goal of this note is to fully develop the reasoning presented in [10, 11]
for $s$-channel DM models with vector or axial-vector mediator exchange. In
Section 2 we demonstrate that a naive calculation of mono-$W$ signatures in
spin-$1$ simplified models leads to unphysical predictions, in particular the
violation of unitarity at large energies. We point out that this issue also
affects mono-jet searches once EW corrections are included in the signal
prediction. In Section 3, we propose a simple solution to this problem based
on imposing certain requirements on the coupling structure of the simplified
model. A different solution is derived in Section 4 by considering mono-$W$
signatures in the SM, where the potentially dangerous terms are cancelled via
interference with a diagram involving a triple gauge boson interaction. We
discuss how this solution can be extended to simplified DM models. In Section
5 it is then shown that this second solution is indeed the one naturally
incorporated in $Z^{\prime}$ models that obtain non-universal couplings to
quarks from mixing with the SM gauge bosons. In Section 6, we briefly discuss
how unitarity violation can be parameterised by non-renormalisable
interactions. The main findings of this note are summarised in Section 7,
while the impact of EW corrections in an existing mono-jet search are studied
in Appendix A.

##  2 The mono-$\boldsymbol{W}$ problem in simplified models

In order to illustrate the issues that can appear in the calculation of
mono-$W$ signals in simplified DM models, we consider a spin-1 mediator
$Z^{\prime}_{\mu}$ with mass $M_{Z^{\prime}}$ and a Dirac DM particle $\chi$
with mass $m_{\text{DM}}$. We write the relevant interaction terms in the
Lagrangian as

$\begin{split}\mathcal{L}=&-Z^{\prime}_{\mu}\,\bar{\chi}\left[g_{\text{DM}}^{V}\gamma^{\mu}+g_{\text{DM}}^{A}\gamma^{\mu}\gamma_{5}\right]\chi\\\\[2.84526pt]
&-\sum_{f=q,l,\nu}Z^{\prime}_{\mu}\,\bar{f}\left[g_{f}^{V}\gamma^{\mu}+g_{f}^{A}\gamma^{\mu}\gamma_{5}\right]f\,.\end{split}$
(1)

In the context of such a simplified model, the mono-$W$ signature arises at
the LHC from the process $pp\rightarrow W+Z^{\prime}$, followed by the decay
of the $Z^{\prime}$ into DM particles. A possible Feynman diagram leading to a
$l+E_{T,\rm miss}$ final state is shown on the left-hand side of Figure 1.
Notice that the interactions (1) also give rise to mono-jet signatures. There
are in fact both QCD and EW corrections to jets + $E_{\rm T,miss}$, which are
exemplified by the middle and right diagram in Figure 1, respectively.

The fundamental problem discussed in this paper is illustrated in Figure 2,
which shows the predicted $E_{T,\rm miss}$ spectrum for leptonically (left
panel) and hadronically (right panel) decaying $W$ produced in association
with an invisibly decaying $Z^{\prime}$. Our event generation has been
performed at leading order with MadGraph5_aMCNLO [12] starting from the
implementation of the spin-$1$ simplified model presented in [13] and utilises
NNPDF2.3 parton distribution functions (PDFs) [14]. We consider $pp$
collisions at a center-of-mass energy of $\sqrt{s}=8\,{\rm TeV}$, employ
$M_{Z^{\prime}}=1\,{\rm TeV}$, $\Gamma_{Z^{\prime}}=56.5\,{\rm GeV}$, $m_{\rm
DM}=10\,{\rm GeV}$, $g_{\rm DM}^{V}=1$ and set all axial-vector and leptonic
couplings to zero. The yellow curves in both panels correspond to SS couplings
$g^{V}_{u}=g^{V}_{d}=0.25$, while the red curves represent the OS coupling
choice $g^{V}_{u}=-g^{V}_{d}=-0.25$. It is evident from the left plot in
Figure 2 that the mono-$W$ predictions for the two coupling choices do not
simply differ by a rescaling factor, but that the predicted differential cross
sections are fundamentally different. Most notably the $E_{T,\rm miss}$
spectrum in the OS case is significantly harder than that for the SS choice.

Figure 2:  Example of a $E_{T,\rm miss}$ spectrum of a mono-$W$ (top panel)
and mono-jet (bottom panel) signal arising in the spin-$1$ simplified model.
See text for additional explanations.

The suspicious high-energy behaviour of the $pp\to W+Z^{\prime}$ amplitudes in
the OS case becomes particularly obvious in the right panel of Figure 2, where
we compare the $E_{T,\rm miss}$ spectrum for hadronically decaying $W$ to the
$E_{T,\rm miss}$ spectrum for a conventional mono-jet analysis, where the jets
arise from QCD interactions (see middle and right diagram in Figure 1). In the
SS case, one observes that the $E_{T,\rm miss}$ spectra for both processes are
similar, but that the overall cross section for the QCD process (blue curve)
is significantly larger than that of the EW process (yellow curve), the latter
amounting to a relative correction of a few percent only. In the OS case, on
the other hand, we find the EW contribution (red curve) to the differential
cross section of $pp\to{\rm jets}+E_{T,\rm miss}$ to be harder and to dominate
over the QCD process (green curve) at sufficiently large $E_{T,\rm miss}$
values. This finding turns out to be essentially unaffected by the parton
shower, hadronisation corrections and detector effects. For spin-$1$
simplified model realisations (1) with OS couplings, the EW channel $pp\to
W\hskip 0.7113pt(\to q\bar{q}^{\prime})+Z^{\prime}\hskip
0.7113pt(\to\chi\bar{\chi})$ and not the QCD process $pp\to{\rm
jets}+Z^{\prime}\hskip 0.7113pt(\to\chi\bar{\chi})$ is therefore expected to
give rise to the majority of events in the high-$E_{T,\rm miss}$ signal
regions (SRs) of LHC Run I mono-jet searches. That this expectation is indeed
correct is shown in Appendix A.

The paradoxical observation that an EW contribution appears to produce harder
mono-jet events than a QCD process casts doubt on the validity of the spin-$1$
simplified model introduced in (1). Indeed, we will show below that the harder
$E_{T,\rm miss}$ spectrum in the OS case is due to contributions that grow
with energy and therefore potentially violate unitarity. The aim of this note
is to accurately describe the nature of this problem and propose a number of
ways of solving it by appropriately modifying or restricting the spin-$1$
simplified model under consideration.333Note that, since the interactions (1)
of the $Z^{\prime}$ boson explicitly break the EW symmetry, the corresponding
Ward identities are no longer satisfied by $W$ bosons in the final state. As a
result, the Goldstone boson equivalence theorem does not hold, i.e. one does
not obtain the same result at high energies when replacing $W_{L}$ fields by
Goldstone bosons. In other words, since the gauge symmetry is broken, unitary
gauge and Feynman gauge are not equivalent and one cannot simply remove the
mono-$W$ problem by calculating cross sections in Feynman gauge.

##  3 Unitarity violation and coupling structures

It is well-known that the production of longitudinal gauge bosons (such as
$Z^{\prime}_{L}$) from fermions can potentially violate unitarity at large
energies (see for instance [15]). For example, it was shown in [16] that for
the spin-$1$ simplified model considered in (1) the process
$\chi+\bar{\chi}\rightarrow Z^{\prime}_{L}+Z^{\prime}_{L}$ violates unitarity
at large energies for non-zero axial-vector couplings, unless a second dark
Higgs is added to the theory. Since the amplitude is proportional to the
fermion mass, the corresponding process with light quarks only violates
unitarity at very large energies and can therefore be disregarded.

The process $u+\bar{d}\rightarrow W^{+}_{L}+Z^{\prime}_{L}$, however, can lead
to much lower scales of unitarity violation. Following [17, 16], we calculate
the helicity matrix elements for the case of a left-handed up-quark and a
right-handed anti-down quark (both of which are charged under $SU(2)_{L}$).
Since the initial state has a total spin of 1, the process cannot proceed at
$s$-wave ($\mathcal{M}^{0}=0$). For the $J=1$ partial wave (or $p$-wave), we
find in the limit $s\rightarrow\infty$

$\mathcal{M}^{1}=\frac{g\hskip 0.7113ptV_{ud}\hskip 0.7113pts}{96\pi\hskip
0.7113ptM_{W}\hskip
0.7113ptM_{Z^{\prime}}}\left(g^{A}_{u}-g^{A}_{d}-g^{V}_{u}+g^{V}_{d}\right)\,,$
(2)

with $g$ the $SU(2)_{L}$ gauge coupling and $M_{W}$ the mass of the $W$ boson.
The helicity matrix element is found to grow proportional to $s$ at large
energies. For

$s>\frac{48\pi\hskip 0.7113ptM_{W}\hskip
0.7113ptM_{Z^{\prime}}}{g\,|V_{ud}|\left|g^{A}_{u}-g^{A}_{d}-g^{V}_{u}+g^{V}_{d}\right|}\,,$
(3)

it will therefore violate the requirement for perturbative unitarity,
$|\mathcal{M}^{J}|<1/2$.444The matrix element corresponding to
$u+\bar{d}\rightarrow W^{+}_{L}+Z^{\prime}_{T}$ with $Z^{\prime}_{T}$ denoting
a transversal $Z^{\prime}$ boson scales as
$\sqrt{s}/M_{W}\left(g^{A}_{u}-g^{A}_{d}-g^{V}_{u}+g^{V}_{d}\right)$ in the
high-energy limit. In consequence, the resulting mono-$W$ cross section does
not grow with $s$, but becomes a constant in the $s\to\infty$ limit. Compared
to the purely longitudinal contribution (2), transversal modes thus provide
only a subleading source of unitarity violation.

In the SS case with vanishing axial-vector couplings, one finds for the
combination of couplings appearing in (3) that
$g^{A}_{u}-g^{A}_{d}-g^{V}_{u}+g^{V}_{d}=0$ and hence the potentially
unitarity-violating term vanishes identically. For OS couplings, on the other
hand, one has $g^{A}_{u}-g^{A}_{d}-g^{V}_{u}+g^{V}_{d}=-2g^{V}_{u}\neq 0$, if
the axial-vector couplings are zero. The resulting mono-$W$ and mono-jet cross
sections then grow proportional to $s$, which explains both the large
enhancement of the total fiducial cross section and the harder $E_{T,\rm
miss}$ spectrum observed in Section 2 for the OS case.

The most obvious way to avoid this problem is to simply impose the requirement

$g^{A}_{u}-g^{A}_{d}-g^{V}_{u}+g^{V}_{d}=0\,,$ (4)

on the coupling structure of the spin-$1$ simplified model (1). In order to
understand this requirement intuitively, let us rewrite the left-handed and
right-handed couplings in terms of vector and axial-vector couplings. For the
sign conventions adopted in (1), one obtains

$g_{q}^{L}=g_{q}^{V}-g_{q}^{A}\,,\qquad g_{q}^{R}=g_{q}^{V}+g_{q}^{A}\,.$ (5)

In terms of left-handed and right-handed couplings, the requirement (4) hence
simply becomes $g_{u}^{L}=g_{d}^{L}$, i.e. left-handed up-type and down-type
quarks should couple in the same way to the $Z^{\prime}$ mediator.

Notice that the latter requirement would be automatically fulfilled, if one
wrote the quark couplings to the spin-$1$ mediator in such a way so as to
preserve the EW symmetry,

$\displaystyle{\cal L}_{Z^{\prime}q\bar{q}}=$
$\displaystyle-\sum_{u,d}\,Z^{\prime}_{\mu}\left[g_{u}\hskip
1.42262pt\bar{u}_{R}\hskip 0.7113pt\gamma^{\mu}\hskip
0.28453ptu_{R}+g_{d}\hskip 1.42262pt\bar{d}_{R}\hskip
0.7113pt\gamma^{\mu}\hskip 0.28453ptd_{R}\right]$
$\displaystyle-\sum_{Q}\,Z^{\prime}_{\mu}\,g_{Q}\hskip
1.42262pt\bar{Q}_{L}\hskip 0.7113pt\gamma^{\mu}\hskip 0.28453ptQ_{L}\,,$ (6)

where $Q_{L}=(u_{L},d_{L})^{T}$ is the left-handed quark doublet and $u_{R}$
and $d_{R}$ are the right-handed quark singlets. Each quark flavour has
couplings described by the three parameters $g_{Q}$, $g_{u}$ and $g_{d}$, and
one can arrange for any desired relative sign between either the vector or
axial-vector couplings, but cannot choose them independently from one
another.555We emphasise that there are additional considerations that may
further restrict the possible choices of $g_{Q}$, $g_{u}$ and $g_{d}$. For
example, to be consistent with minimal flavour violation, the same couplings
should be chosen for all three quark families [18]. Moreover, models with non-
zero axial-vector couplings (i.e. $g_{u}\neq g_{Q}$ or $g_{d}\neq g_{Q}$)
require additional structure in the Higgs sector, which may lead to strong
constraints from Higgs measurements or EW precision observables [16].

It is however conceivable that the coupling structure of the $Z^{\prime}$
boson is modified by EWSB and that the spin-1 mediator ends up coupling
differently to left-handed up-type and down-type quarks. In this case, the
unitarity issue of amplitudes like $u+\bar{d}\to W^{+}+Z^{\prime}$ needs to be
solved in a different fashion.

##  4 Interference with additional diagrams

In the previous section, we have shown that unitarity violation in the process
$pp\rightarrow W+Z^{\prime}$ can be avoided if we require that the
$Z^{\prime}$ couples with equal strength to left-handed up-type and down-type
quarks. We can compare this requirement to the couplings of the SM $Z$ boson,
whose interactions with SM quarks can be written as

$\mathcal{L}_{Zq\bar{q}}=-\sum_{q}Z_{\mu}\,\bar{q}\left[g_{Z,q}^{V}\hskip
0.7113pt\gamma^{\mu}+g_{Z,q}^{A}\hskip
0.7113pt\gamma^{\mu}\gamma_{5}\right]q\,,$ (7)

with

$\begin{split}g_{Z,q}^{V}&=\frac{g}{2\hskip
0.7113ptc_{w}}\left(T_{3,q}-2Q_{q}\hskip
0.7113pts_{w}^{2}\right)\,,\\\\[2.84526pt] g_{Z,q}^{A}&=-\frac{g}{2\hskip
0.7113ptc_{w}}\hskip 1.42262ptT_{3,q}\,,\end{split}$ (8)

where $s_{w}$ and $c_{w}$ are the sine and cosine of the weak mixing angle,
$Q_{q}$ denotes the electromagnetic quark charge in units of the elementary
charge $e$ and the weak isospin is given by $T_{3,u}=1/2$ and $T_{3,d}=-1/2$.
It follows that the difference $\Delta g_{Z}$ of the left-handed couplings of
the $Z$ boson to up-type and down-type quarks reads

$\displaystyle\Delta g_{Z}$
$\displaystyle=g_{Z,u}^{L}-g_{Z,d}^{L}=g_{Z,u}^{V}-g_{Z,u}^{A}-g_{Z,d}^{V}+g_{Z,d}^{A}$
$\displaystyle=g\hskip 0.7113ptc_{w}\neq 0\,.$ (9)

In other words, the requirement (4) is not satisfied by the $Z$ boson within
the SM and as a result, one might expect that the process $pp\rightarrow W+Z$
also shows a unitarity-violating behaviour at high energies.

Indeed, when considering only the two $t$-channel diagrams shown on the left-
hand side and in the middle of Figure 3 the resulting amplitude grows with
energy. In the SM, however, a third diagram is present, which involves a
$s$-channel $W$ boson radiating off a $Z$ boson. The corresponding graph is
depicted on the right-hand side in the latter figure. It involves a triple
gauge boson vertex of the form

$\mathcal{L}_{WWZ}=ig_{WWZ}\hskip 0.7113pt{\cal T}_{WWZ}\,,$ (10)

with $g_{WWZ}=g\hskip 0.7113ptc_{w}=\Delta g_{Z}$ and

$\begin{split}{\cal
T}_{WWV}&=\Bigl{[}\left(\partial_{\mu}W_{\nu}^{+}-\partial_{\nu}W_{\mu}^{+}\right)W^{\mu-}V^{\nu}\\\\[2.84526pt]
&\phantom{xxx}-\left(\partial_{\mu}W_{\nu}^{-}-\partial_{\nu}W_{\mu}^{-}\right)W^{\mu+}V^{\nu}\\\\[2.84526pt]
&\phantom{xxx}+\frac{1}{2}\left(\partial_{\mu}V_{\nu}-\partial_{\nu}V_{\mu}\right)\\\\[2.84526pt]
&\phantom{xxxx}\times\left(W^{\mu+}W^{\nu-}-W^{\mu-}W^{\nu+}\right)\Bigr{]}\,.\end{split}$
(11)

Combining the three different contributions, one finds that in the high-energy
limit the square of the $u+\bar{d}\to W^{+}+Z$ matrix element is given by

$|{\cal
M}|^{2}=\frac{3g^{4}c_{w}^{4}|V_{ud}|^{2}}{32M_{W}^{2}}\left(d_{1}+d_{2}-2d_{3}\right)^{2}\,s^{2}\hskip
0.7113pt\sin^{2}\theta\,,$ (12)

where $V_{ud}$ denotes the relevant element of the Cabibbo-Kobayashi-Maskawa
matrix, $\theta$ is the scattering angle in the center-of-mass frame and
$d_{1}$, $d_{2}$ and $d_{3}$ label the contributions from the three different
graphs. From the result (12) it is readily seen that while the individual
diagrams all diverge, their sum remains finite at large energies, since
$d_{1}+d_{2}-2d_{3}=0$ for $d_{1}=d_{2}=d_{3}=1$. Notice that the cancellation
of unitarity-violating terms among the diagrams of Figure 3 is not at all
accidental, but a direct consequence of the local $SU(2)_{L}$ gauge invariance
of the SM.

Figure 3:  The three tree-level diagrams that contribute to $u+\bar{d}\to
W^{+}+Z$ in the SM. The left and middle graph involves only couplings between
quarks and gauge bosons, while the right diagram contains a triple gauge boson
vertex.

The above observation motivates a second solution to the mono-$W$
problem.666This solution was also briefly discussed in [11]. If the difference
$\Delta g$ between the left-handed couplings of the $Z^{\prime}$ boson to up-
type and down-type quarks takes the non-zero value

$\Delta g=g_{u}^{L}-g_{d}^{L}=g_{u}^{V}-g_{u}^{A}-g_{d}^{V}+g_{d}^{A}\,,$ (13)

it is still possible to avoid unitarity-violating processes by adding
additional gauge boson interactions of the form

$\displaystyle\Delta\mathcal{L}=ig_{WWZ^{\prime}}\hskip 1.42262pt{\cal
T}_{WWZ^{\prime}}\,,$ (14)

with $g_{WWZ^{\prime}}=\Delta g$ to the Lagrangian (1) of the spin-$1$
simplified model. Here ${\cal T}_{WWZ^{\prime}}$ denotes the $WWZ^{\prime}$
analogue of the $WWV$ interactions introduced in (11).

##  5 A ${\boldsymbol{Z}^{\prime}}$ boson with mixing

In this section we discuss a class of models where the additional interactions
(14) between the $W$ bosons and the $Z^{\prime}$ boson are automatically
generated in such a way that they cancel the potentially unitarity-violating
contributions to the mono-$W$ signal. The basic idea is that the mass
eigenstate $Z^{\prime}_{\mu}$ arises from an interaction eigenstate vector
field $X_{\mu}$ with the SM $U(1)_{Y}$ $B_{\mu}$ field and the neutral
component $W^{3}_{\mu}$ of the $SU(2)_{L}$ weak fields through mixing

$\displaystyle\left(\begin{array}[]{c}B_{\mu}\\\ W^{3}_{\mu}\\\
X_{\mu}\end{array}\right)=\left(\begin{array}[]{ccc}N_{11}&N_{12}&N_{13}\\\
N_{21}&N_{22}&N_{23}\\\
N_{31}&N_{32}&N_{33}\end{array}\right)\left(\begin{array}[]{c}A_{\mu}\\\
Z_{\mu}\\\ Z_{\mu}^{\prime}\end{array}\right)\,.$ (24)

Here $A_{\mu}$ and $Z_{\mu}$ denotes the physical photon and neutral massive
gauge boson fields of the SM. At this point, we leave the entries $N_{ij}$ of
the above mixing matrix unspecified. Examples for how such a mixing matrix can
be generated and a very brief discussion of the constraints on the $N_{ij}$
that arise from EW precision tests will be given below.

The couplings of the interaction eigenstate $X_{\mu}$ can be parameterised in
the same way as in the spin-$1$ simplified model introduced in Section 2.
Including only DM and quark couplings, we write

$\begin{split}\mathcal{L}=&-X_{\mu}\,\bar{\chi}\left[f_{\text{DM}}^{V}\gamma^{\mu}+f_{\text{DM}}^{A}\gamma^{\mu}\gamma_{5}\right]\chi\\\\[2.84526pt]
&-\sum_{q}X_{\mu}\,\bar{q}\left[f_{q}^{V}\gamma^{\mu}+f_{q}^{A}\gamma^{\mu}\gamma_{5}\right]q\,,\end{split}$
(25)

and assume that the direct couplings $f_{q}^{V}$ and $f_{q}^{A}$ satisfy the
requirement (4), i.e. they fulfill
$f_{u}^{L}-f_{d}^{L}=f^{A}_{u}-f^{A}_{d}-f^{V}_{u}+f^{V}_{d}=0$. In the
absence of mixing unitarity-violating contributions to the mono-$W$ process
will hence be absent.

However, as we have discussed in Section 4, the SM $Z$ boson couples
differently to left-handed up-type and down-type quarks and in the presence of
mixing the $Z^{\prime}$ can inherit this difference. The induced couplings
between the $Z^{\prime}$ and SM fermions are discussed for example in
[19].777Notice that the sign convention for the axial-vector couplings used in
[19] is opposite to the one adopted in the present work. For the left-handed
quark couplings, one obtains

$\begin{split}g_{u}^{L}&=f^{L}_{u}\hskip 0.7113ptN_{33}+\frac{g}{2}\hskip
0.7113ptN_{23}+\frac{e}{6c_{w}}\hskip 0.7113ptN_{13}\,,\\\\[2.84526pt]
g_{d}^{L}&=f^{L}_{d}\hskip 0.7113ptN_{33}-\frac{g}{2}\hskip
0.7113ptN_{23}+\frac{e}{6c_{w}}\hskip 0.7113ptN_{13}\,.\end{split}$ (26)

Since by assumption $f^{L}_{u}-f^{L}_{d}=0$ and due to the SS of the terms
proportional to $N_{13}$, the coupling difference (13) thus takes the simple
form

$\Delta g=g_{u}^{L}-g_{d}^{L}=g\hskip 0.7113ptN_{23}\,.$ (27)

In addition to the couplings (26), the mixing matrix $N$ entering (24) also
induces couplings of the $Z^{\prime}$ to other SM particles [20]. In
particular, the SM $WWZ$ vertex (10) will lead to $WWZ^{\prime}$ interactions.
Since the triple gauge boson vertex arises from the non-abelian $SU(2)_{L}$,
the induced coupling turns out to be proportional to

$g_{WWZ^{\prime}}=g\hskip 0.7113ptN_{23}\,,$ (28)

and of the form (11). The resulting $WWZ^{\prime}$ interaction is therefore
precisely that given in (14). As discussed in Section 4, the unitarity issue
of the $u+\bar{d}\to W^{+}+Z^{\prime}$ channel is then solved by the
interference of the two diagrams with $t$-channel quark exchange and the
single graph involving a $s$-channel $W$ boson.

It is also important to realise that in the class of models discussed in this
section, the mono-$W$ channel in addition receives a contribution from $pp\to
W+Z\hskip 0.7113pt(\to\chi\bar{\chi})$, since the SM $Z$-boson field $Z_{\mu}$
inherits the direct DM couplings $f_{\rm DM}^{V}$ and $f_{\rm DM}^{A}$ from
the interaction field $X_{\mu}$ through mixing. Like (12) this contribution
however remains unitary at large energies as a result of the $SU(2)_{L}$ gauge
invariance of the SM. This shows that the mono-$W$ problem is indeed absent in
$Z^{\prime}$ models with mixing.

Since we did not specify the origin of the gauge boson mixings $N_{ij}$, our
conclusions apply to a range of different scenarios. For example, a non-zero
entry $N_{23}$ can be generated via kinetic mixing [21]

$\mathcal{L}_{F^{\prime}B}=-\frac{1}{2}\hskip
0.7113pt\sin\epsilon\,F^{\prime}_{\mu\nu}B^{\mu\nu}\;,$ (29)

between the $U(1)^{\prime}$ field strength $F^{\prime}_{\mu\nu}$ and the SM
hypercharge field strength $B_{\mu\nu}$. In this case, one finds in the limit
of small mixing angles $\epsilon\ll 1$ the result

$N_{23}=-c_{w}\hskip
0.7113pts_{w}\,\epsilon\,\frac{M_{Z}^{2}}{M_{Z}^{2}-M_{Z^{\prime}}^{2}}\,,$
(30)

with $M_{Z}$ denoting the mass of the SM $Z$ boson. Another option is to
consider so-called mass mixing between the $Z$ and the $Z^{\prime}$ [22],
which is described by

$\mathcal{L}_{Z^{\prime}Z}=\delta m^{2}\,Z^{\mu}Z^{\prime}_{\mu}\,.$ (31)

The $N_{23}$ entry in the gauge boson mixing matrix takes for mass mixing the
form

$N_{23}=-c_{w}\,\frac{\delta m^{2}}{M_{Z}^{2}-M_{Z^{\prime}}^{2}}\,.$ (32)

We note for ${Z^{\prime}}$-boson masses in the TeV range both $\epsilon$ and
$\delta m^{2}/M_{Z}^{2}$ have to be very small (of the order of $10^{-2}$ or
below) in order not to violate bounds from EW precision measurements (see e.g.
[19]). Nevertheless, there are also other ways to generate gauge boson mixing,
for example from left-right symmetric models [23].

##  6 Parametrising unitarity violation with non-renormalisable interactions

In the previous section, we have seen that $Z^{\prime}$ models with mixing
remain unitary at high energies, if they are $SU(2)_{L}$ symmetric in the
absence of mixing. While this treatment covers most of the simplest models for
a spin-1 mediator, there may of course be more complicated models with
additional heavy particles at a new-physics scale $\Lambda$. In such a model,
the couplings of the mediator to the left-handed up-type and down-type quarks
may be realised via higher-dimensional operators such as

$\begin{split}{\cal
L}_{Z^{\prime}QH}=-\sum_{u,d}Z^{\prime}_{\mu}\Bigl{[}&\frac{1}{\Lambda_{u}^{2}}\,(\bar{Q}_{L}\tilde{H})\hskip
0.7113pt\gamma^{\mu}(\tilde{H}^{\dagger}Q_{L})\\\\[2.84526pt]
&+\frac{1}{\Lambda_{d}^{2}}\,(\bar{Q}_{L}H)\hskip
0.7113pt\gamma^{\mu}(H^{\dagger}Q_{L})\Bigr{]}\,,\end{split}$ (33)

which effectively decouple the left-handed up-type and down-type quark
couplings. Here $\Lambda_{u,d}$ are two different suppression scales and
$\tilde{H}=i\sigma_{2}H^{\ast}$ with $H$ denoting the $SU(2)_{L}$ Higgs
doublet. In such a case, the difference $\Delta g$ between the left-handed
quark couplings can be parameterised by [10]

$\Delta g=\frac{v^{2}}{\Lambda^{2}}\,,$ (34)

with $v\simeq 246\,{\rm GeV}$ the Higgs vacuum expectation value.888In the
specific case (33), the new-physics scale $\Lambda$ appearing in (34) is given
by $\Lambda=(\left|\Lambda_{u}^{-2}-\Lambda_{d}^{-2}\right|/2)^{-1/2}$. Since
the details of the ultraviolet completion are not specified, there is no
direct connection between the $Z^{\prime}$ coupling to SM quarks and the
$Z^{\prime}$ coupling to $W$ bosons.

As a result, the mono-$W$ problem may still be present in such models. Notice
that this is not a theoretical inconsistency, because unitarity violation is a
generic feature of EFTs once the partonic centre-of-mass energy $\sqrt{s}$ (or
more generically the momentum transfer relevant for the process under
consideration) approaches the new-physics scale $\Lambda$.999The effect of
unitarity violation in EFT models for DM has been discussed in [24, 25, 26].
The mono-$W$ problems thus does not preclude meaningful calculations within
the EFT parameterisation (34), provided the scale $\Lambda$ is sufficiently
large. In fact, one can translate the bound (3) into a limit on the new-
physics scale $\Lambda$ for a given partonic centre-of-mass energy:

$\Lambda>\left(\frac{s\hskip 1.42262ptv}{24\pi
M_{Z^{\prime}}}\right)^{1/2}\,.$ (35)

This means for instance that if one wants to study a model with
$M_{Z^{\prime}}=1\,{\rm TeV}$, one needs to impose $\Lambda>0.75\,{\rm TeV}$
in order to ensure that the model does not violate perturbative unitarity up
to $\sqrt{s}=13\,{\rm TeV}$. This in turn limits the magnitude of the
couplings in the OS case to $|g^{V}_{u}|=|g^{V}_{d}|<0.05$.

Of course, the $\sqrt{s}=13\,{\rm TeV}$ LHC does on average not probe partonic
center-of-mass energies of $13\,{\rm TeV}$, so a smaller new-physics scale
$\Lambda$ and therefore larger OS couplings may in practice be acceptable. A
proper treatment is likely to require a self-consistent truncation procedure,
which discards all events with a momentum transfer larger than $\Lambda$ (as
discussed for instance in [27, 28]). For example, to consistently study the OS
example considered in Section 2, one should only consider kinematic
configurations with $\sqrt{s}<6.1\,{\rm TeV}$ (for $M_{Z^{\prime}}=1\,{\rm
TeV}$). Such a truncation is expected to lead to a notable reduction of the
LHC sensitivity to mono-$W$ signatures.

##  7 Summary

In the context of simplified models with spin-1 mediators, we have discussed
the so-called mono-$W$ problem, i.e. the violation of perturbative unitarity
in the process $pp\rightarrow W+E_{T,\rm miss}$, and proposed three different
solutions. These are

  1. 1.

Restricted coupling structures: The mono-$W$ problem is absent if the coupling
structure of the simplified model is restricted in such a way that the spin-1
mediator couples equally to left-handed up-type and down-type quarks
$g_{u}^{L}=g_{d}^{L}$. This requirement is automatically satisfied if the
interaction between the $Z^{\prime}$ and the quarks are formulated in an
$SU(2)_{L}$ invariant way.

  2. 2.

Additional interactions: For $g_{u}^{L}\neq g_{d}^{L}$, the mono-$W$ problem
can still be solved if an additional interaction between the spin-1 mediator
and SM $W$ bosons is introduced, which is proportional to the difference
between the left-handed couplings. The interference between $t$-channel quark
exchange and $s$-channel $W$-boson diagrams then removes the potentially
dangerous contributions.

  3. 3.

Small couplings: If no additional interactions between the $Z^{\prime}$ and
$W$ bosons are present, the mono-$W$ process will violate perturbative
unitarity at large energies. The scale of unitarity violation may however be
larger than the energies probed in a given process, if the difference between
the couplings to up-type and down-type quark couplings is sufficiently small.

We have argued that the first solution is realised in the simplest
$Z^{\prime}$ models, whereas the second solution appears naturally in
$Z^{\prime}$-boson scenarios with mixing. If a more general model is to be
considered, the third solution allows for a general assessment of the validity
of a chosen coupling structure in a given experimental environment.

The arguments presented in our work strongly suggest that the sensitivity of
mono-$W$ searches to the parameter space of spin-$1$ simplified models cannot
exceed that of the mono-jet channel. This conclusion seems to hold
irrespectively of how the unitarity problem in the $pp\to W+E_{T,\rm miss}$
process is tamed. The same verdict has been previously reached in [10, 11] for
the EFT case and $t$-channel simplified DM models with coloured scalar
exchange. It has also been argued in [11] that the mono-$W$ problem is not
present in $s$-channel models if they are formulated in a gauge-invariant way,
and our note explicitly proves this conjecture.

## Acknowledgments

We thank the members of the LHC DM working group for discussions that
triggered this project and Nicole Bell for carefully reading a close-to-final
draft of the manuscript and for her useful comments. UH acknowledges the
hospitality and support of the CERN theory division. FK is supported by the
German Science Foundation (DFG) under the Collaborative Research Center (SFB)
676 Particles, Strings and the Early Universe. TMPT is supported in part by
NSF grant PHY-1316792 and by the University of California, Irvine through a
Chancellor’s Fellowship.

##  Appendix A EW contributions to mono-jet searches

In this appendix, we discuss the numerical impact of EW contributions of the
form $pp\to W\hskip 0.7113pt(\to q\bar{q}^{\prime})+Z^{\prime}\hskip
0.7113pt(\to\chi\bar{\chi})$ to mono-jet signatures at the LHC, comparing the
effects to the corrections associated to the QCD process $pp\to{\rm
jets}+Z^{\prime}\hskip 0.7113pt(\to\chi\bar{\chi})$. To make this exercise
concrete, we consider the recent ATLAS search [29], which is based on
$20.3\,{\rm fb}^{-1}$ of $\sqrt{s}=8\,{\rm TeV}$ data. We use MadGraph5_aMCNLO
with NNPDF2.3 PDFs to generate partonic events. The simulated parton-level
events were showered with PYTHIA 6 [30] and analysed with the publicly
available code CheckMATE [31], which relies on DELPHES 3 [32] as a fast
detector simulation. To cluster jets we used FastJet [33] employing the
anti-$k_{t}$ algorithm [34] with radius parameter $R=0.4$.

In the ATLAS analysis the following preselection criteria are imposed. Events
are required to have a reconstructed primary vertex, $E_{T,\rm miss}>150\,{\rm
GeV}$ and at least one jet with $p_{T}>30\,{\rm GeV}$ and $|\eta|<4.5$ in the
final state. Events that do not pass certain jet quality requirements or do
contain charged leptons or isolated tracks are rejected. Events having a
leading jet with $p_{T}>120\,{\rm GeV}$ and $|\eta|<2.0$ are selected, if the
leading-jet $p_{T}$ and the $E_{T,\rm miss}$ satisfy $p_{T}/E_{T,\rm
miss}>0.5$. Furthermore, the requirement $\Delta\phi({\rm jet},\vec{p}_{T,\rm
miss})>1.0$ on the azimuthal separation between the direction of the missing
transverse momentum and that of each of the selected jets is imposed. Nine
distinct SRs are considered with the following $E_{T,\rm miss}$ thresholds
$\\{150,200,250,300,350,400,500,600,700\\}\,{\rm GeV}$.

Figure 4:  Fiducial mono-jet cross sections for the SR1 to SR9 selections. The
top (bottom) panel shows the predictions for a spin-$1$ simplified model with
only vector (axial-vector) couplings. See text for further explanations.

Our results for the fiducial mono-jet cross section corresponding to SR1 to
SR9 are shown in the two panels of Figure 4. The plots are based on
$M_{Z^{\prime}}=1\,{\rm TeV}$, $m_{\rm DM}=10\,{\rm GeV}$ and all leptonic
couplings are taken to be zero. In the left (right) panel we present the case
of pure vector couplings $g_{\rm DM}^{V}=1$ and $|g^{V}_{u}|=|g^{V}_{d}|=0.25$
(pure axial-vector couplings $g_{\rm DM}^{A}=1$ and
$|g^{A}_{u}|=|g^{A}_{d}|=0.25$). The corresponding total decay width of the
$Z^{\prime}$ boson are $\Gamma_{Z^{\prime}}=56.5\,{\rm GeV}$ and
$\Gamma_{Z^{\prime}}=55.5\,{\rm GeV}$ assuming a minimal width. One first
observes that the obtained results are to good approximation independent of
whether vector or axial-vector exchange is considered. Second, while the
fiducial cross sections corresponding to the QCD contributions (blue and green
markers) are within statistical uncertainties identical for the SS and OS
choices, in the EW case the OS signal strengths (red markers) are, depending
on the considered SR, larger than the SS predicitons (yellow markers) by a
factor of around 25 to 400. One also sees that, as a result of the
constructive interference in the OS case, starting with SR4 the EW fiducial
cross sections $\sigma_{\rm fid}^{\rm EW,OS}$ surpass the corresponding QCD
predictions. For the SR7 selection corresponding to $E_{T,\rm miss}>500\,{\rm
GeV}$, we find for instance $\sigma_{\rm fid}^{\rm EW,OS}/\sigma_{\rm
fid}^{\rm QCD,OS}\simeq 330\%$ ($\sigma_{\rm fid}^{\rm EW,OS}/\sigma_{\rm
fid}^{\rm QCD,OS}\simeq 350\%$) in the vector (axial-vector) case. The ratio
of the SS fiducial cross sections in SR7 amount instead to $\sigma_{\rm
fid}^{\rm EW,SS}/\sigma_{\rm fid}^{\rm QCD,SS}\simeq 2\%$, independently of
the choice of mediator. This example nicely illustrates that the mono-$W$
problem is not only relevant for searches that look for a $W$ boson and large
$E_{T,\rm miss}$, but would in general also impact mono-jet searches, if it is
not resolved.

## References

## References

  * [1] Y. Bai and T. M. P. Tait, Searches with Mono-Leptons,  Phys. Lett. B723 (2013) 384–387, [1208.4361].
  * [2] ATLAS, G. Aad et al., Search for dark matter in events with a hadronically decaying W or Z boson and missing transverse momentum in $pp$ collisions at $\sqrt{s}=$ 8 TeV with the ATLAS detector,  Phys. Rev. Lett. 112 (2014), no. 4 041802, [1309.4017].
  * [3] CMS, Search for dark matter in the mono-lepton channel with pp collision events at center-of-mass energy of 8 TeV,  CMS-PAS-EXO-13-004.
  * [4] CMS, V. Khachatryan et al., Search for physics beyond the standard model in final states with a lepton and missing transverse energy in proton-proton collisions at $\sqrt{s}$ = 8 TeV,  Phys. Rev. D91 (2015), no. 9 092005, [1408.2745].
  * [5] M. Beltran, D. Hooper, E. W. Kolb, Z. A. C. Krusberg, and T. M. P. Tait, Maverick dark matter at colliders,  JHEP 09 (2010) 037, [1002.4137].
  * [6] J. Goodman, M. Ibe, A. Rajaraman, W. Shepherd, T. M. P. Tait, et al., Constraints on Dark Matter from Colliders,  Phys. Rev. D82 (2010) 116010, [1008.1783].
  * [7] Y. Bai, P. J. Fox, and R. Harnik, The Tevatron at the Frontier of Dark Matter Direct Detection,  JHEP 12 (2010) 048, [1005.3797].
  * [8] ATLAS, Search for dark matter produced in association with a hadronically decaying vector boson in pp collisions at $\sqrt{s}=$ 13 TeV with the ATLAS detector at the LHC,  ATLAS-CONF-2015-080.
  * [9] ATLAS, G. Aad et al., Search for new particles in events with one lepton and missing transverse momentum in $pp$ collisions at $\sqrt{s}$ = 8 TeV with the ATLAS detector,  JHEP 09 (2014) 037, [1407.7494].
  * [10] N. F. Bell, Y. Cai, J. B. Dent, R. K. Leane, and T. J. Weiler, Dark matter at the LHC: Effective field theories and gauge invariance,  Phys. Rev. D92 (2015), no. 5 053008, [1503.07874].
  * [11] N. F. Bell, Y. Cai, and R. K. Leane, Mono-W Dark Matter Signals at the LHC: Simplified Model Analysis,  JCAP 1601 (2016), no. 01 051, [1512.00476].
  * [12] J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, et al., The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations,  JHEP 07 (2014) 079, [1405.0301].
  * [13] M. Backović, M. Krämer, F. Maltoni, A. Martini, K. Mawatari, et al., Higher-order QCD predictions for dark matter production at the LHC in simplified models with s-channel mediators,  Eur. Phys. J. C75 (2015), no. 10 482, [1508.05327].
  * [14] R. D. Ball et al., Parton distributions with LHC data,  Nucl. Phys. B867 (2013) 244–289, [1207.1303].
  * [15] K. S. Babu, J. Julio, and Y. Zhang, Perturbative unitarity constraints on general $W^{\prime}$ models and collider implications,  Nucl. Phys. B858 (2012) 468–487, [1111.5021].
  * [16] F. Kahlhoefer, K. Schmidt-Hoberg, T. Schwetz, and S. Vogl, Implications of unitarity and gauge invariance for simplified dark matter models,  JHEP 02 (2016) 016, [1510.02110].
  * [17] M. S. Chanowitz, M. A. Furman, and I. Hinchliffe, Weak Interactions of Ultraheavy Fermions. 2.,  Nucl. Phys. B153 (1979) 402.
  * [18] J. Abdallah et al., Simplified Models for Dark Matter Searches at the LHC,  Phys. Dark Univ. 9-10 (2015) 8–23, [1506.03116].
  * [19] M. T. Frandsen, F. Kahlhoefer, A. Preston, S. Sarkar, and K. Schmidt-Hoberg, LHC and Tevatron Bounds on the Dark Matter Direct Detection Cross-Section for Vector Mediators,  JHEP 07 (2012) 123, [1204.3839].
  * [20] E. J. Chun, J.-C. Park, and S. Scopel, Dark matter and a new gauge boson through kinetic mixing,  JHEP 02 (2011) 100, [1011.3300].
  * [21] B. Holdom, Two U(1)’s and Epsilon Charge Shifts,  Phys. Lett. B166 (1986) 196\.
  * [22] K. S. Babu, C. F. Kolda, and J. March-Russell, Implications of generalized $Z$–$Z^{\prime}$ mixing,  Phys. Rev. D57 (1998) 6788–6792, [hep-ph/9710441].
  * [23] J. Shu and J. Yepes, Diboson excess and $Z^{\prime}$-predictions via left-right non-linear Higgs,  (2016) 1601.06891.
  * [24] I. M. Shoemaker and L. Vecchi, Unitarity and Monojet Bounds on Models for DAMA, CoGeNT, and CRESST-II,  Phys. Rev. D86 (2012) 015023, [1112.5457].
  * [25] P. J. Fox, R. Harnik, R. Primulando, and C.-T. Yu, Taking a Razor to Dark Matter Parameter Space at the LHC,  Phys. Rev. D86 (2012) 015010, [1203.1662].
  * [26] M. Endo and Y. Yamamoto, Unitarity Bounds on Dark Matter Effective Interactions at LHC,  JHEP 06 (2014) 126, [1403.6610].
  * [27] G. Busoni, A. De Simone, E. Morgante, and A. Riotto, On the Validity of the Effective Field Theory for Dark Matter Searches at the LHC,  Phys. Lett. B728 (2014) 412–421, [1307.2253].
  * [28] G. Busoni, A. De Simone, J. Gramling, E. Morgante, and A. Riotto, On the Validity of the Effective Field Theory for Dark Matter Searches at the LHC, Part II: Complete Analysis for the $s$-channel,  JCAP 1406 (2014) 060, [1402.1275].
  * [29] ATLAS, G. Aad et al., Search for new phenomena in final states with an energetic jet and large missing transverse momentum in pp collisions at $\sqrt{s}=$ 8 TeV with the ATLAS detector,  Eur. Phys. J. C75 (2015), no. 7 299, [1502.01518]. [Erratum: Eur. Phys. J. C75 (2015), no. 9 408].
  * [30] T. Sjostrand, S. Mrenna, and P. Z. Skands, PYTHIA 6.4 Physics and Manual,  JHEP 05 (2006) 026, [hep-ph/0603175].
  * [31] M. Drees, H. Dreiner, D. Schmeier, J. Tattersall, and J. S. Kim, CheckMATE: Confronting your Favourite New Physics Model with LHC Data,  Comput. Phys. Commun. 187 (2014) 227–265, [1312.2591].
  * [32] DELPHES 3, J. de Favereau et al., DELPHES 3, A modular framework for fast simulation of a generic collider experiment,  JHEP 02 (2014) 057, [1307.6346].
  * [33] M. Cacciari, G. P. Salam, and G. Soyez, FastJet User Manual,  Eur. Phys. J. C72 (2012) 1896, [1111.6097].
  * [34] M. Cacciari, G. P. Salam, and G. Soyez, The anti-$k_{t}$ jet clustering algorithm,  JHEP 04 (2008) 063, [0802.1189].


Physics-based block preconditioning for   interaction

imcs]Max<EMAIL_ADDRESS>imcs]Ivo<EMAIL_ADDRESS>imcs]Alexander<EMAIL_ADDRESS>imcs,dsc]Matthias<EMAIL_ADDRESS>[imcs]Institute for Mathematics and Computer-Based Simulation, Universität der Bundeswehr München,
Werner-Heisenberg-Weg 39, D-85577 Neubiberg, Germany
[dsc]Data Science & Computing Lab, Universität der Bundeswehr München,
Werner-Heisenberg-Weg 39, D-85577 Neubiberg, Germany

[cor1]corresponding author

This paper presents a scalable physics-based block preconditioner for  models in  interaction and their application in engineering.
In particular, it studies the linear systems arising from a regularized mortar-type approach
for embedding geometrically exact beams into solid continua.

Due to the lack of block diagonal dominance of the arising $2\times 2$ block system,
an approximate block factorization preconditioner is used.
It exploits the sparsity structure of the beam  to construct a sparse approximate inverse,
which is then not only used to explicitly form an approximation of the Schur complement,
but also acts as a smoother within the prediction step of the arising SIMPLE-type preconditioner.
The correction step utilizes an algebraic multigrid method.
Although, for now, the beam  is tackled by a one-level method only,
the  nature of the computationally demanding correction step delivers a scalable preconditioner in practice.
In numerical test cases,
the influence of different algorithmic parameters on the quality of the sparse approximate inverse is studied
and the weak scaling behavior of the proposed preconditioner on up to 1000 MPI ranks is demonstrated,
before the proposed preconditioner is finally applied for the analysis of steel-reinforced concrete structures in civil engineering.

Physics-based block preconditioning algebraic multigrid sparse approximate inverse  modeling  interaction

§ PHYSICS-BASED BLOCK PRECONDITIONING FOR  INTERACTION

The construction of a preconditioner, that captures the coupling interactions properly and is tailored to the specific matrix properties,
is crucial for an efficient and scalable solution process. Preconditioners based on approximate block factorizations have shown to be
suited for similar problem types such as contact problems [57],
incompressible flow [37, 15],
ac:FSI [26, 23],
or magneto-hydro dynamics [8, 38, 39].

§.§ Schur complement based block factorization

Due to the lack of block diagonal dominance as discussed in sec:LinearSystems and the ill-conditioning of the matrix,
we resort to a block factorization preconditioner based on a Schur complement.
Specifically, we perform an approximate block factorization based on a SIMPLE-type approach [37].
To this end, the preconditioning matrix $\precMatrix$ can be written as

\begin{equation*}
\precMatrix
\left(
\begin{matrix}
\matBeam & 0 \\
\matSolidBeam & \approxMat{\schurMatrix} \\
\end{matrix}
\right)
\left(
\begin{matrix}
\idmat & \inv{\approxMat{\matBeam}} \trans{\matBeamSolid}\\
0 & \frac{1}{\simpleRelaxFac}\idmat \\
\end{matrix}
\right)
\left(
\begin{matrix}
\matBeam & \matBeam \inv{\approxMat{\matBeam}} \trans{\matBeamSolid}\\
\matSolidBeam & \frac{1}{\simpleRelaxFac} \matSolid + \left(1-\frac{1}{\simpleRelaxFac}\right)\matSolidBeam \inv{\approxMat{\matBeam}} \trans{\matBeamSolid}\\
\end{matrix}
\right)
\approx
\left(
\begin{matrix}
\matBeam & \trans{\matBeamSolid}\\
\matSolidBeam & \matSolid\\
\end{matrix}
\right)
\end{equation*}

with $\approxMat{\matBeam}$ and $\approxMat{\schurMatrix}$ denoting easy-to-invert approximations of the  $\matBeam$ from (<ref>)
and the Schur complement $\schurMatrix := \matSolid + \matSolidBeam \inv{\matBeam} \trans{\matBeamSolid}$, respectively,
and $\simpleRelaxFac > 0$ being a relaxation factor.
Expressing the application of the preconditioner as a  iteration over the index $\indSweep$,
one application of the preconditioner yields

\begin{equation*}
\left(
\begin{matrix}
\indexedBeamSweep{\mao{\sol}}{\indSweep+1} \\
\indexedSolidSweep{\mao{\sol}}{\indSweep+1} \\
\end{matrix}
\right)
\left(
\begin{matrix}
\indexedBeamSweep{\mao{\sol}}{\indSweep} \\
\indexedSolidSweep{\mao{\sol}}{\indSweep} \\
\end{matrix}
\right)
\left(
\begin{matrix}
\matBeam & \matBeam\inv{\approxMat{\matBeam}}\trans{\matBeamSolid}\\
\matSolidBeam & \frac{1}{\simpleRelaxFac} \matSolid + \left(1-\frac{1}{\simpleRelaxFac}\right)\matSolidBeam \inv{\approxMat{\matBeam}} \trans{\matBeamSolid}\\
\end{matrix}
\right)^{-1}
\left[
\left(
\begin{matrix}
\indexedBeam{\mao{\rhs}}\\
\indexedSolid{\mao{\rhs}}\\
\end{matrix}
\right)
\left(
\begin{matrix}
\matBeam & \trans{\matBeamSolid} \\
\matSolidBeam & \matSolid \\
\end{matrix}
\right)
\left(
\begin{matrix}
\indexedBeamSweep{\mao{\sol}}{\indSweep}\\
\indexedSolidSweep{\mao{\sol}}{\indSweep}\\
\end{matrix}
\right)
\right],
\end{equation*}

which is usually tackled via a a predictor-corrector scheme by first solving an equation related to the beam contribution and afterwards a second one related to the Schur complement.

Both quality and computational cost of the preconditioner are mainly governed by the choice $\approxMat{\matBeam}$ and $\approxMat{\schurMatrix}$ to approximate $\matBeam$ and $\schurMatrix$, respectively.
In traditional Schur complement based block preconditioners,
the inverse $\inv{\approxMat{\matBeam}} \approx \inv{\matBeam}$ appearing due to the block factorization and also in the Schur complement calculation itself
is often approximated by some diagonal matrix,
since the inversion of a diagonal matrix comes at a negligible computational cost.
The most simple approach is to base the inverse approximation on the diagonal part of the $\nRows\times\nRows$ matrix $\matBeam$
resulting in $\approxMat{\matBeam} = \diag{\matBeamEntry{\indRow}{\indRow}}, \indRow = 1,\hdots,\nRows$.
Another well-known method referred to as SIMPLEC takes the row sums of $\matBeam$ [55], reading

\begin{equation*}
\approxMat{\matBeam} = \diag{\sum_{\indCol=1}^\nRows \abs{\matBeamEntry{\indRow}{\indCol}}}, \quad \indRow = 1,\hdots,\nRows.
\end{equation*}

However, such simple diagonal approximations cannot be used in the present scenario
since $\matBeam$ lacks diagonal dominance,  sec:BlockDiagonalDominance.
Since $\matBeam$ resembles the  related to the beam equations, which themselves yield a block diagonal structure of $\matBeam$,
a more sophisticated explicit approximation scheme of the inverse can be applied taking this particular sparsity structure into account (see sec:SPAI).

To this end,
the choice of $\approxMat{\schurMatrix}$ is governed by $\approxMat{\matBeam}$ and the relaxation parameter $\simpleRelaxFac$.

§.§ Explicit sparse inverse approximation

Although in general the inverse of the sparse matrix $\matBeam$ cannot be expected to be sparse as well,
explicit sparse approximate inverses aim to create an explicit matrix representation $\spaiOf{\matBeam}$ of the approximation of the exact inverse $\inv{\matBeam}$,
that itself is still a sparse matrix.
Ideally, $\nnzOf{\spaiOf{\matBeam}}$, the number of  of $\spaiOf{\matBeam}$, does not exceed $\nnzOf{\matBeam}$,
since a matrix-vector product with $\spaiOf{\matBeam}$ must be performed at each Krylov iteration [18].
Although incomplete LU factorizations fall into that category,
they require considerable effort to parallelize [10].

To take advantage of the block diagonal structure of $\matBeam$,
we pursue a fully parallelizable approach to construct an explicit sparse approximate inverse $\spaiOf{\matBeam}$ of the matrix $\matBeam$
based on the minimization of the Frobenius norm of the residual matrix $\matBeam\spaiOf{\matBeam} - \idmat$,
see [18, 48].
By choosing an appropriate sparsity pattern $\graphOf{\spaiOf{\matBeam}}$ for the sparse approximate inverse $\spaiOf{\matBeam}$ from the set $\graph$ describing all known patterns,
the following least-squares problem needs to be solved:

\begin{equation}
\min_{\graphOf{\spaiOf{\matBeam}} \; \in \; \graph} \normFrobenius{\matBeam\spaiOf{\matBeam} - \idmat}.
\label{eq:SparseApproximateInverseMinimizationProblem}
\end{equation}

We represent the minimization procedure to compute $\spaiOf{\matBeam}$ by the operation $\spaiOf{\matBeam} \leftarrow \spaiMinimizeOp{\matBeam}{\graphOf{\spaiOf{\matBeam}}}$.

To this end, (<ref>) requires to select an appropriate sparsity pattern $\graphOf{\spaiOf{\matBeam}}$
and a practical approach to the minimization of the Frobenius norm in a distributed memory environment,
which we will address in sec:SPAIGraphSelectionsec:SPAIFrobeniusNorm, respectively.

§.§.§ Selection of a sparsity pattern for the sparse approximate inverse calculation

The main challenge in (<ref>) is the selection of a sparsity pattern $\graphOf{\spaiOf{\matBeam}}$
to be used as input into the minimization procedure.
An appropriate pattern needs to contain enough information of the inverse by retaining high values,
but should also act as a filter to remove small entries in order to reduce fill-in.

Straightforward approaches are based on a static sparsity pattern selection and include the choices $\graphOf{\spaiOf{\matBeam}} := \graphOf{\matBeam}$
and $\graphOf{\spaiOf{\matBeam}} := \graphOf{\trans{\matBeam}}$,
which are easy to obtain, but do not guarantee a good approximation quality.
Especially in cases with a partially known sparsity structure of the inverse,  for block-diagonal matrices such as $\matBeam$ in (<ref>),
more advanced static selections are able to deliver a satisfying approximation,
yet not requiring dynamic pattern selection approaches as proposed in [18].

In this work, we follow the static pattern selection proposed in [9]
and use powers of a sparsified version $\threshMat{\graphOf{\matBeam}}$ of the graph $\graphOf{\matBeam}$ of the original matrix $\matBeam$
to obtain an enriched sparsity pattern to be used for the minimization in (<ref>).
First, $\threshMat{\graphOf{\matBeam}}$ is obtained through a thresholding of $\graphOf{\matBeam}$ based on the entries in $\matBeam$
and using a drop-off tolerance $\spaiDropTol$.
We represent this thresholding by the filter operation

\begin{align*}
\threshMat{\graphOf{\matBeam}} \leftarrow \thresholdingOf{{\matBeam}}{\spaiDropTol}{\graph}
\end{align*}

delivering individual entries $\threshMat{\graphElementOf{\matBeam}{\indRow}{\indCol}}$
of the filtered graph via

\begin{equation}
\threshMat{\graphElementOf{\matBeam}{\indRow}{\indCol}} :=
\begin{cases}
1 \quad \text{if} \; \indRow=\indCol \; \text{or} \; \abs{\diagEntry{\indRow}^{-\frac{1}{2}} \matBeamEntry{\indRow}{\indCol} \diagEntry{\indRow}^{-\frac{1}{2}}} > \spaiDropTol, \\
0 \quad \text{otherwise},
\end{cases}
\quad \text{where} \quad
\diagEntry{\indRow} :=
\begin{cases}
\abs{\matBeamEntry{\indRow}{\indRow}} \quad \text{if} \; \abs{\matBeamEntry{\indRow}{\indRow}} > 0, \\
1 \quad \text{otherwise}.
\end{cases}
\label{eq:SparseApproximateInverseSparsification}
\end{equation}

The additional Jacobi scaling in (<ref>) simplifies the thresholding and choosing of $\spaiDropTol$
if $\matBeam$ is poorly scaled.
In a second step, a refined sparsity pattern $\graphOf{\matBeam^\spaiRefinementLevel}$
is calculated by taking powers $\spaiRefinementLevel$ of the sparsified graph $\threshMat{\graphOf{\matBeam}}$,

\begin{align*}
\graphOf{\spaiOf{\matBeam}} \leftarrow \spaiRefinementOp{\threshMat{\graphOf{\matBeam}}}{\spaiRefinementLevel}.
\end{align*}

The matrix $\matBeam^\spaiRefinementLevel$ is never calculated explicitly, but the powers are directly computed on its sparsified graph $\threshMat{\graphOf{\matBeam}}$.

§.§.§ Evaluation of the Frobenius norm on a parallel computer

Due to the 2-norm compatibility of the Frobenius norm, the problem can be decoupled into a sum of Euclidian norms, reading

\begin{equation*}
\min_{\graphOf{\spaiOf{\matBeam}} \; \in \; \graph} \normFrobenius{\matBeam\spaiOf{\matBeam} - \idmat}
= \sum_{\spaiIndRow=1}^\nRows \min_{\graphOf{\spaiOf{\matBeam}} \; \in \; \graph}
\normTwo{\left(\matBeam\spaiOf{\matBeam} - \idmat\right) \indexedRow{\mao{e}}{\spaiIndRow}}^2
= \sum_{\spaiIndRow=1}^\nRows \min_{\graphOf{\indexedRow{\spaiOf{\matBeam}}{\spaiIndRow}} \; \in \; \indexedRow{\graph}{\spaiIndRow}}
\normTwo{\matBeam\indexedRow{\spaiOf{\matBeam}}{\spaiIndRow} - \indexedRow{\mao{e}}{\spaiIndRow}}^2,
\end{equation*}

which can be solved independently for each row $\spaiIndRow=1,\hdots,\nRows$ of $\spaiOf{\matBeam}$.
Since the matrix $\matBeam$ is usually stored in a row-wise distribution,
where each parallel process stores a subset of all rows of $\matBeam$,
this decomposition renders the method, besides an initial communication step, inherently parallel.
To this end, calculating a row of the approximate inverse means solving a small least-squares problem by applying a dense QR factorization.

§.§.§ Practical algorithm

In practice, one usually does not solve (<ref>) directly,
but rather combines it with some pre- and post-operations,
foremost the graph computation from sec:SPAIGraphSelection and the handling of the Frobenius norm outlined in sec:SPAIFrobeniusNorm.
In addition, it appears beneficial to perform a post-filtering of $\spaiOf{\matBeam}$ by dropping all entries
with $\abs{\indexedRowCol{\spaiOf{\matBeamElement}}{\indRow}{\indCol}} < \spaiDropTol \quad \forall \indRow,\indCol=1,\hdots,\nRows$
to further reduce fill-in of $\inv{\approxMat{\matBeam}}$ and, thus, the cost of applying the preconditioner [9].
The post-filtering will be denoted by the operator $\thresholdingOf{\spaiOf{\matBeam}}{\spaiDropTol}{\matBeam}$.

We summarize all the steps to compute the sparse approximate inverse $\spaiOf{\matBeam}$ in fig:SpaiAlgorithm.

Optional (gray) and mandatory (orange) steps of the ac:SPAI computation
and its flow of information with computed data (in boxes) and user parameters (in circles)

Therein, optional steps are marked in gray, while the mandatory computation of $\spaiOf{\matBeam}$ is highlighted in orange.
User-given data to configure the individual steps is depicted in circles,
while computed input and output data is put into rectangular boxes.
The arrows indicate the flow of data between the individual steps.

Even tough the method can be steered quite effectively by the refinement level $\spaiRefinementLevel$ and threshold tolerance $\spaiDropTol$,
there are still several problems to avoid.
The method can still produce rather dense matrices with poor approximation quality.
Depending on the threshold parameter $\spaiDropTol$,
an aggressive dropping of values might result in a loss of information,
which again results in a poor result and the approximation of the inverse might even be singular.

§.§ One-level approach for the predictor step

The first linear equation to be solved inside the block preconditioner resembling the predictor step consists of the beam equation with matrix $\matBeam$:

\begin{equation}
\matBeam \indexedBeamSweep{\mao{\sol}}{\indSweep+\frac{1}{2}} = \indexedBeamSweep{\mao{\rhs}}{\indSweep} - \trans{\matBeamSolid} \indexedSolidSweep{\mao{\sol}}{\indSweep}.
\label{eq:PredictorStep}
\end{equation}

Following the assumptions made in sec:BlockDiagonalDominance, traditional smoothing approaches are not applicable
to the linear system without a major loss in convergence properties.

Based on the explicit sparse approximate inverse calculation for the approximation of the Schur complement,
a rather good approximation $\inv{\approxMat{\matBeam}}$ is already available.
Following [7], the sparse approximate inverse is reused as smoother by applying
the following  iteration over index $\indSpaiSweep$ to (<ref>):

\begin{equation}
\indexedBeamMultiSweep{\mao{\sol}}{\indSweep+\frac{1}{2}}{\indSpaiSweep + 1}
= \indexedBeamMultiSweep{\mao{\sol}}{\indSweep+\frac{1}{2}}{\indSpaiSweep} - \inv{\approxMat{\matBeam}}\left(\indexedBeamMultiSweep{\mao{\rhs}}{\indSweep}{\indSpaiSweep} - \trans{\matBeamSolid} \indexedSolidMultiSweep{\mao{\sol}}{\indSweep}{\indSpaiSweep}\right).
\label{eq:SparseApproximateInverseSmoother}
\end{equation}

A comparison to more traditional smoothers,  Jacobi and  methods, can also be found in the original publication [7].

Expanding the solution procedure of the predictor step to an ac:AMG approach to solve the beam related
equation is challenging, as there are still several open questions around the construction of ac:AMG hierarchies for beam models:
the beams are represented as ac:1D elements, for which a suitable coarsening scheme to form aggregates
for coarser ac:MG levels is still part of active research.
Even with appropriate coarsening, fibers discretized with only a few elements would quickly form single-node
aggregates being suboptimal for the restrictor and prolongator construction as well as effectively stalling the coarsening process.
The addition of rotational components into the beam equation introduces another difficulty,
as nodes with a different number of ac:DOF exist, hindering the application of conventional aggregation strategies.
To project and, thus, treat important error modes correctly on coarser levels,
a proper near nullspace has to be constructed, which highly depends on the underlying beam formulation.
For the use cases presented in this publication, which mostly feature a small number of beam elements used to discretize a fiber,
the ac:SPAI smoother from (<ref>) proved to be very competitive
in terms of approximation quality of the inverse, computational efficiency and weak scaling behavior,
see sec:Experiments.

§.§ Multilevel approach for the correction step

The computationally most demanding part of the preconditioning algorithm involves the solution procedure for the
Schur complement equation given as

\begin{equation*}
\approxMat{\schurMatrix} \indexedSolidSweep{\mao{\sol}}{\indSweep+\frac{1}{2}} = \indexedSolidSweep{\mao{\rhs}}{\indSweep} - \matSolid \indexedSolidSweep{\mao{\sol}}{\indSweep} - \matSolidBeam \indexedBeamSweep{\mao{\sol}}{\indSweep+\frac{1}{2}}.
\end{equation*}

Due to the explicit sparse approximation of the inverse of $\matBeam$
used to form the approximate Schur complement $\approxMat{\schurMatrix} := \matSolid + \matSolidBeam \inv{\approxMat{\matBeam}} \trans{\matBeamSolid}$,
the resulting matrix $\approxMat{\schurMatrix}$ is rather
dense compared to calculations using one of the diagonal approximation approaches for $\matBeam$
given in sec:SchurComplementBlockFactorization.
The inverse of the Schur complement is approximated by a standard aggregation-based ac:AMG method.
As level smoother, a one-level domain decomposition with overlap $\iluOverlap$ and an incomplete LU factorization with
fill-in $\iluLevelFill$ and thresholding $\iluDropTol$ of small entries  is applied, often abbreviated as ILUT [45].
To accelerate convergence, a smoothing of the tentative prolongation operator basis functions can be done,
also known as ac:SAAMG [56].
This however leads to a higher fill-in of the coarse system matrix representations increasing the computational cost
of the preconditioner setup. Especially the Galerkin product for the calculation of the coarse level operator and
the incomplete LU factorization smoother are negatively influenced by the additional fill-in. Furthermore,
the operator smoothing is classically based on a Jacobi method, which heavily relies on diagonal dominance.
Thus, in some cases it makes sense to skip the prolongator smoothing and take advantage of the robustness of
ac:PAAMG by applying aggregate-wise constant basis functions [54].

§.§ Physics-based block preconditioning algorithm

The components described in sec:Preconditioner are now put together to form the final preconditioner summarized in alg:Preconditioner.
The calculation of the explicit sparse inverse of the beam matrix and the approximate Schur complement mark the pre-computation steps.
Based on the sweep index $\indSweep$ of the preconditioner,
the main computation loop consists of the solution of two linear systems related to the beam and Schur complement equation.
The beam equation is solved in the predictor step by applying a SPAI smoother.
For the Schur complement equation, an algebraic multigrid method is used to obtain the solution.
In a final step, the solution vector is updated and later returned.

Preconditioner for   coupling

Preconditioner($\simpleRelaxFac, \indSweep_{max}$)

// Pre-compute SPAI of $\matBeam$: $\inv{\approxMat{\matBeam}} \leftarrow \thresholdingOf{\spaiMinimizeOp{\matBeam}{\spaiRefinementOp{\thresholdingOf{{\matBeam}}{\spaiDropTol}{\graph}}{\spaiRefinementLevel}}}{\spaiDropTol}{\matBeam}$

// Form explicit, approximate Schur complement: $\approxMat{\schurMatrix} = \matSolid + \matSolidBeam \inv{\approxMat{\matBeam}} \trans{\matBeamSolid}$

$\indSweep = 1 \dots \indSweep_{max}$

// Prediction step: solve for $\indexedBeamSweep{\mao{\sol}}{\indSweep+\frac{1}{2}}$ with SPAI smoother

$\indSpaiSweep = 1 \dots \indSpaiSweep_{max}$

$\indexedBeamMultiSweep{\mao{\sol}}{\indSweep+\frac{1}{2}}{\indSpaiSweep + 1}
= \indexedBeamMultiSweep{\mao{\sol}}{\indSweep+\frac{1}{2}}{\indSpaiSweep} - \inv{\approxMat{\matBeam}}\left(\indexedBeamMultiSweep{\mao{\rhs}}{\indSweep}{\indSpaiSweep} - \trans{\matBeamSolid} \indexedSolidMultiSweep{\mao{\sol}}{\indSweep}{\indSpaiSweep}\right)$

// Correction step: solve for $\indexedSolidSweep{\mao{\sol}}{\indSweep+\frac{1}{2}}$ with AMG$\approxMat{\schurMatrix} \indexedSolidSweep{\mao{\sol}}{\indSweep+\frac{1}{2}} = \indexedSolidSweep{\mao{\rhs}}{\indSweep} - \matSolid \indexedSolidSweep{\mao{\sol}}{\indSweep} - \matSolidBeam \indexedBeamSweep{\mao{\sol}}{\indSweep+\frac{1}{2}}$

// Update step: update solution variables $\indexedSolidSweep{\mao{\sol}}{\indSweep+1} \leftarrow \indexedSolidSweep{\mao{\sol}}{\indSweep} + \simpleRelaxFac \indexedSolidSweep{\mao{\sol}}{\indSweep+\frac{1}{2}}$ $\indexedBeamSweep{\mao{\sol}}{\indSweep+1} \leftarrow \indexedBeamSweep{\mao{\sol}}{\indSweep+\frac{1}{2}} - \simpleRelaxFac \inv{\approxMat{\matBeam}} \trans{\matBeamSolid} \indexedSolidSweep{\mao{\sol}}{\indSweep+\frac{1}{2}}$

// Return solution vector $\left(\indexedBeamSweep{\mao{\sol}}{\indSweep_{max}}, \indexedSolidSweep{\mao{\sol}}{\indSweep_{max}}\right)$

§ ACKNOWLEDGEMENTS

The work described in this contribution has been funded by dtec.bw - Digitalization and Technology Research Center of the Bundeswehr
under the project “hpc.bw - Competence Platform for High Performance Computing”.
dtec.bw is funded by the European Union – NextGenerationEU.
Partial support has been provided by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)
within the project
“Stable discretization methods and scalable solvers for embedded fibersolid coupling” (project number 528397555)
and by the Bayerisches Verbundforschungsprogramm (BayVFP) – Digitalisierung within the project “IFB-Individuelle Fließfertigung für Betonfertigteile” (project number DIK-2201-0016 / DIK0411/03).
The authors gratefully acknowledge the computing resources provided by the Data Science & Computing Lab at the University of the Bundeswehr Munich.

§ DATA AVAILABILITY STATEMENT

Building blocks of the AMG preconditioner developed and applied in this study are openly available in Trilinos at <https://github.com/trilinos/Trilinos> [2].
All other data that support the findings of this study are available from the corresponding author upon reasonable request.

§ DECLARATION OF COMPETING INTEREST

The authors declare that they have no known competing financial interests or personal relationships
that could have appeared to influence the work reported in this paper.

[1]
BACI: A Comprehensive Multi-Physics Simulation Framework.
Accessed: 2023-07-14.

[2]
The Trilinos Project.
Accessed: 2023-07-14.

[3]
B. D. Agarwal, L. J. Broutman, and K. Chandrashekhara.
Analysis and Performance of Fiber Composites.
Wiley, 2017.

[4]
J. Ao, M. Zhou, and B. Zhang.
A dual mortar embedded mesh method for internal interface problems
with strong discontinuities.
International Journal for Numerical Methods in Engineering,
123(22):5652–5681, 2022.

[5]
E. Béchet, N. Moës, and B. I. Wohlmuth.
A stable lagrange multiplier space for stiff interface conditions
within the extended finite element method.
International Journal for Numerical Methods in Engineering,
78(8):931–954, 2009.

[6]
L. Berger-Vergiat, C. A. Glusa, G. Harper, J. J. Hu, M. Mayr, A. Prokopenko,
C. M. Siefert, R. S. Tuminaro, and T. A. Wiesner.
MueLu User's Guide.
Technical Report SAND2023-12265, Sandia National Laboratories,
Albuquerque, NM (USA) 87185, 2023.

[7]
O. Bröker and M. J. Grote.
Sparse approximate inverse smoothers for geometric and algebraic
Applied Numerical Mathematics, 41(1):61–80, 2002.

[8]
L. Chacón.
An optimal, parallel, fully implicit Newton–Krylov solver for
three-dimensional viscoresistive magnetohydrodynamics.
Physics of Plasmas, 15:056103, 2008.

[9]
E. Chow.
Parallel implementation and practical use of sparse approximate
inverse preconditioners with a priori sparsity patterns.
The International Journal of High Performance Computing
Applications, 15(1):56–74, 2001.

[10]
E. Chow and A. Patel.
Fine-Grained Parallel Incomplete LU Factorization.
SIAM Journal on Scientific Computing, 37(2):C169–C193, 2015.

[11]
E. C. Cyr, J. N. Shadid, and R. S. Tuminaro.
Teko: A Block Preconditioning Capability with Concrete Example
Applications in Navier–Stokes and MHD.
SIAM Journal on Scientific Computing, 38(5):S307–S331, 2016.

[12]
N. Dimola, M. Kuchta, K.-A. Mardal, and P. Zunino.
Robust Preconditioning of mixed-dimensional PDEs on 3d-1d domains
coupled with Lagrange multipliers.
preprint on arXiv: <https://arxiv.org/abs/2307.08431>, 2023.

[13]
D. C. Drucker and W. Prager.
Soil mechanics and plastic analysis for limit design.
Quarterly of Applied Mathematics, 10(2):157–166, 1952.

[14]
D. Durville.
Finite element simulation of textile materials at mesoscopic scale.
In Finite element modelling of textiles ans textile composites,
Saint Petersburg, 2007.

[15]
H. C. Elman, V. E. Howle, J. N. Shadid, R. Shuttleworth, and R. S. Tuminaro.
A taxonomy and comparison of parallel block multi-level
preconditioners for the incompressible Navier–Stokes equations.
Journal of Computational Physics, 227(3):1790–1808, 2008.

[16]
D. Feingold and R. Varga.
Block diagonally dominant matrices and generalizations of the
Gershgorin Theorem.
Pacific Journal of Mathematics, 12(4):1241–1250, 1962.

[17]
B. Flemisch and B. I. Wohlmuth.
Stable Lagrange multipliers for quadrilateral meshes of curved
interfaces in 3D.
Computer Methods in Applied Mechanics and Engineering,
196(8):1589–1602, 2007.

[18]
M. J. Grote and T. Huckle.
Parallel preconditioning with sparse approximate inverses.
SIAM Journal on Scientific Computing, 18(3):838–853, 1997.

[19]
N. Hagmeyer, M. Mayr, and A. Popp.
A fully coupled regularized mortar-type finite element approach for
embedding one-dimensional fibers into three-dimensional fluid flow.
International Journal for Numerical Methods in Engineering,
published online ahead of print:e7435, 2024.

[20]
N. Hagmeyer, M. Mayr, I. Steinbrecher, and A. Popp.
One-way coupled fluid-beam interaction: Capturing the effect of
embedded slender bodies on global fluid flow and vice versa.
Advanced Modeling and Simulation in Engineering Sciences, 9:9,

[21]
M. Hautefeuille, C. Annavarapu, and J. E. Dolbow.
Robust imposition of Dirichlet boundary conditions on embedded
International Journal for Numerical Methods in Engineering,
90(1):40–64, 2012.

[22]
M. A. Heroux, R. A. Bartlett, V. E. Howle, R. J. Hoekstra, J. J. Hu, T. G.
Kolda, Lehoucq, K. R. Long, R. P. Pawlowski, E. T. Phipps, A. G. Salinger,
H. K. Thornquist, R. S. Tuminaro, J. M. Willenbring, A. Williams, and K. S.
An Overview of the Trilinos Project.
ACM Transactions on Mathematical Software, 31(3):397–423,

[23]
D. Jodlbauer, U. Langer, and T. Wick.
Parallel block-preconditioned monolithic solvers for fluid-structure
interaction problems.
International Journal for Numerical Methods in Engineering,
117(6):623–643, 2019.

[24]
S. Kakaletsis, E. Lejeune, and M. Rausch.
The mechanics of embedded fiber networks.
Journal of the Mechanics and Physics of Solids, 181:105456,

[25]
U. Khristenko, S. Schuß, M. Krüger, F. Schmidt, B. Wohlmuth, and
C. Hesch.
Multidimensional coupling: A variationally consistent approach to
fiber-reinforced materials.
Computer Methods in Applied Mechanics and Engineering,
382:113869, 2021.

[26]
U. Langer and H. Yang.
Robust and efficient monolithic fluid-structure-interaction solvers.
International Journal for Numerical Methods in Engineering,
108(4):303–325, 2016.

[27]
B. Lé, G. Legrain, and N. Moës.
Mixed dimensional modeling of reinforced structures.
Finite Elements in Analysis and Design, 128:1–18, 2017.

[28]
F. Lespagnol, C. Grandmont, P. Zunino, and M. A. Fernandez.
A mixed-dimensional formulation for the simulation of slender
structures immersed in an incompressible flow.
Technical Report 98/2023, MOX Politecnico di Milano, 2023.

[29]
X. S. Li and J. W. Demmel.
SuperLU_DIST: A scalable distributed-memory sparse direct solver
for unsymmetric linear systems.
ACM Transactions on Mathematical Software, 29(2):110–140,

[30]
K.-A. Mardal and R. Winther.
Preconditioning discretizations of systems of partial differential
Numerical Linear Algebra with Applications, 18(1):1–40, 2011.

[31]
M. Mayr, L. Berger-Vergiat, P. Ohm, and R. S. Tuminaro.
Non-invasive multigrid for semi-structured grids.
SIAM Journal on Scientific Computing, 44(4):A2734–A2764, 2022.

[32]
M. Mayr, T. Klöppel, W. A. Wall, and M. W. Gee.
A Temporal Consistent Monolithic Approach to Fluid–Structure
Interaction Enabling Single Field Predictors.
SIAM Journal on Scientific Computing, 37(1):B30–B59, 2015.

[33]
M. Mayr, M. Noll, and M. W. Gee.
A hybrid interface preconditioner for monolithic fluid-structure
interaction solvers.
Advanced Modeling and Simulation in Engineering Sciences, 7:15,

[34]
M. Mayr and A. Popp.
Scalable computational kernels for mortar finite element methods.
Engineering with Computers, 39:3691–3720, 2023.

[35]
C. Meier, A. Popp, and W. A. Wall.
A locking-free finite element formulation and reduced models for
geometrically exact Kirchhoff rods.
Computer Methods in Applied Mechanics and Engineering,
290:314–341, 2015.

[36]
C. A. Meier, A. Popp, and W. A. Wall.
Geometrically Exact Finite Element Formulations for Slender Beams:
Kirchhoff–Love Theory Versus  Theory.
Archives of Computational Methods in Engineering,
26(1):163–243, 2019.

[37]
S. Patankar and D. Spalding.
A calculation procedure for heat, mass and momentum transfer in
three-dimensional parabolic flows.
International Journal of Heat and Mass Transfer,
15(10):1787–1806, 1972.

[38]
E. G. Phillips, H. C. Elman, E. C. Cyr, J. N. Shadid, and R. P. Pawlowski.
A Block Preconditioner for an Exact Penalty Formulation for
Stationary MHD.
SIAM Journal on Scientific Computing, 36(6):B930–B951, 2014.

[39]
E. G. Phillips, J. N. Shadid, E. C. Cyr, H. C. Elman, and R. P. Pawlowski.
Block Preconditioners for Stable Mixed Nodal and Edge finite element
Representations of Incompressible Resistive MHD.
SIAM Journal on Scientific Computing, 38(6):B1009–B1031, 2016.

[40]
A. Popp, M. Gitterle, M. W. Gee, and W. A. Wall.
A dual mortar approach for 3D finite deformation contact with
consistent linearization.
International Journal for Numerical Methods in Engineering,
83(11):1428–1465, 2010.

[41]
M. A. Puso.
A 3D mortar method for solid mechanics.
International Journal for Numerical Methods in Engineering,
59(3):315–336, 2004.

[42]
M. A. Puso and T. A. Laursen.
A mortar segment-to-segment contact method for large deformation
solid mechanics.
Computer Methods in Applied Mechanics and Engineering,
193(6–8):601–629, 2004.

[43]
M. A. Puso and T. A. Laursen.
A mortar segment-to-segment frictional contact method for large
Computer Methods in Applied Mechanics and Engineering,
193(45–47):4891–4913, 2004.

[44]
P. Raghavan and S. Ghosh.
A continuum damage mechanics model for unidirectional composites
undergoing interfacial debonding.
Mechanics of Materials, 37(9):955–979, 2005.

[45]
Y. Saad.
Ilut: A dual threshold incomplete lu factorization.
Numerical Linear Algebra with Applications, 1(4):387–402,

[46]
Y. Saad.
Iterative Methods for Sparse Linear Systems.
SIAM, Philadelphia, PA, USA, 2003.

[47]
Y. Saad and M. H. Schultz.
GMRES: A Generalized Minimal Residual Algorithm for Solving
Nonsymmetric Linear Systems.
SIAM Journal on Scientific and Statistical Computing,
7(3):856–869, 1986.

[48]
M. Sedlacek.
Sparse Approximate Inverse for Preconditioning, Smoothing, and
PhD thesis, Technische Universität München, 2012.

[49]
I. Steinbrecher, N. Hagmeyer, C. Meier, and A. Popp.
A consistent mixed-dimensional coupling approach for 1D cosserat
beams and 2D solid surfaces.
preprint on arXiv: <https://arxiv.org/abs/2210.16010>, 2022.

[50]
I. Steinbrecher, M. Mayr, M. J. Grill, J. Kremheller, C. Meier, and A. Popp.
A mortar-type finite element approach for embedding 1D beams into
3D solid volumes.
Computational Mechanics, 66(6):1377–1398, 2020.

[51]
I. Steinbrecher and A. Popp.
MeshPy – A general purpose 3D beam finite element input
<https://imcs-compsim.github.io/meshpy>, 2021.

[52]
I. Steinbrecher, A. Popp, and C. Meier.
Consistent coupling of positions and rotations for embedding 1D
Cosserat beams into 3D solid volumes.
Computational Mechanics, 69:701–732, 2022.

[53]
I. S. Steinbrecher.
Mixed-dimensional finite element formulations for beam-to-solid
PhD thesis, Universität der Bundeswehr München, 2022.

[54]
S. J. Thomas, S. Ananthan, S. Yellapantula, J. J. Hu, M. Lawson, and M. A.
A comparison of classical and aggregation-based algebraic multigrid
preconditioners for high-fidelity simulation of wind turbine incompressible
SIAM Journal on Scientific Computing, 41(5):S196–S219, 2019.

[55]
J. P. Van Doormaal and G. D. Raithby.
Enhancements of the simple method for predicting incompressible fluid
Numerical Heat Transfer, 7(2):147–163, 1984.

[56]
P. Vaněk, J. Mandel, and M. Brezina.
Algebraic Multigrid By Smoothed Aggregation For Second And Fourth
Order Elliptic Problems.
Computing, 56:179–196, 1996.

[57]
T. A. Wiesner, M. Mayr, A. Popp, M. W. Gee, and W. A. Wall.
Algebraic multigrid methods for saddle point systems arising from
mortar contact formulations.
International Journal for Numerical Methods in Engineering,
122(15):3749–3779, 2021.