††institutetext:  ∗ Institute for Theoretical Particle Physics (TTP),
Karlsruhe Institute of Technology (KIT), D-76131 Karlsruhe, Germany
† Università degli Studi and INFN Roma Tre, Via della Vasca Navale 84,
I-00146, Rome, Italy
‡ INFN, Sezione di Firenze Via G. Sansone 1, 50019 Sesto Fiorentino, Italy,

# Fusing photons into diphoton resonances at Belle II and beyond

Francesca Acanfora∗† Roberto Franceschini† Alessio Mastroddi† Diego Redigolo‡

###### Abstract

We propose a new search for a diphoton resonance in the
$e^{+}e^{-}+\gamma\gamma$ final state at Belle II that improves the expected
reach compared to the $\gamma+\gamma\gamma$ channel in most of the available
mass range. For simplicity we show our results in the simple parameter space
of an ALP coupled solely to Standard Model photons. In addition, we show how
an extension of the forward coverage of Belle II, or another similar
experiment at the high intensity frontier, could improve the reach in our
channel. We show that such a forward extension can be advantageous even with a
loss of a factor 100 in luminosity compared to Belle II.

##  1 Introduction

Diphoton resonances have been an extremely important final state in the
history of particle physics, which culminated with the Standard Model Higgs
boson discovery at ATLAS and CMS CMS:2012qbp ; ATLAS:2012yve . Without any
doubt, finding a new peak in the diphoton invariant mass distribution will
give a decisive hint of new physics. For this reason, independently on the
theory motivations, it is of utmost importance to think about new ways to
enhance the sensitivity of diphoton resonance searches at existing colliders.
This program has been pursued in the recent years in different directions and
lead to a significant extension of the reach on diphoton resonances for
different production channels, with a special emphasis in the relatively
unexplored region of low invariant mass for the diphoton pair Knapen:2016moh ;
Knapen:2017ebd ; Mariotti:2017vtv ; CidVidal:2018blh ; Aloni:2019ruo ;
Wang:2021uyb ; Knapen:2021elo ; ATLAS:2022abz ; Alonso-Alvarez:2023wni ;
CMS:2018erd ; Balkin:2023gya .

Figure 1: Expected sensitivity of Belle II at 95% C.L. to the ALP coupling to
photons $g_{a\gamma\gamma}$ as defined in Eq. (3). In red we show the expected
reach of $e^{+}e^{-}+(\gamma\gamma)_{\text{res}}$ derived here and in orange
the one of $\gamma+(\gamma\gamma)_{\text{res}}$ derived in Ref. Dolan:2017osp
both assuming the present Belle II detector coverage (see Eq. (5)) and
$50\text{ ab}^{-1}$ of integrated luminosity. The details of our analysis are
given in Sec. 2. The green dashed line corresponds to the reach of a
hypothetical experiment with increased forward acceptance (see Eq. (6)) and
identical luminosity. The details are given in Sec. 3. The red error band and
the green error band on both expected reaches correspond to the different
statistical treatment of the empty bins in the four-dimensional binned
likelihood, as detailed in App. A and App. B. The dotted blue lines are
isolines for the sample averaged ALP decay length measured in the lab-frame.
When the ALP decays become long enough, the expected sensitivity of
$e^{+}e^{-}+\text{invisible}$ Acanfora:2023gzr  and of
$\gamma+\text{invisible}$ Dolan:2017osp  are also shown. The gray shaded
region shows existing constraints, assuming $g_{a\gamma\gamma}$ is generated
above the electroweak scale by an ALP coupling to the $U(1)_{Y}$ gauge boson
(see Ref. Dolan:2017osp ; Ferber:2022rsf  for the case where the ALP couples
to $SU(2)$ gauge bosons). These include limits from beam dump experiments
Riordan:1987aw ; Dobrich:2015jyk , LEP and LHC limits recast in Ref.
Knapen:2016moh , constraints from PrimEX data derived in Ref. Aloni:2019ruo ,
the recent limits from the Belle II search on
$\gamma+(\gamma\gamma)_{\text{res}}$ Belle-II:2020jti , from BES III data
BESIII:2022rzz  and from ATLAS Knapen:2017ebd  and CMS CMS:2018erd  searches
on ultra-peripheral heavy ions collisions.

In this paper, we propose a new search for a diphoton resonance at Belle II
taking advantage of the “photon fusion” channel

$e^{+}e^{-}\to e^{+}e^{-}+(\gamma\gamma)_{\text{res}}\,,$ (1)

where $(\gamma\gamma)_{\text{res}}$ indicates the diphoton resonant pair. In
concrete scenarios where a new scalar (or pseudoscalar) is coupled to the SM
photons, we show how this channel can improve the reach with respect to the
more commonly considered “ALP-strahlung” channel Dolan:2017osp

$e^{+}e^{-}\to\gamma+(\gamma\gamma)_{\text{res}}\,.$ (2)

Despite the smaller production cross-section with respect to
$e^{+}e^{-}\to\gamma+(\gamma\gamma)_{\text{res}}$, the $e^{+}e^{-}\to
e^{+}e^{-}+(\gamma\gamma)_{\text{res}}$ process has characteristic features
that make this signal highly distinguishable from the QED background.
Ultimately, these features stem from the higher dimensional observable phase-
space with respect to the “ALP-strahlung” process.

For instance, in the SM process the photons are preferably emitted at small
angle with respect to the $e^{\pm}$ which radiated them and tend to carry
little energy, while the signal allows for a large angular separation and
energetic photons. Moreover, at low invariant masses the photons in the SM
process tend to carry very little energy compared to the signal ones. In the
following, we will elaborate further on the discrimination that can be
attained depending on the invariant mass of the diphoton resonance.

For concreteness, we consider the simple model of an axion-like particle (ALP)
coupled to SM photons, but we expect that our findings can be applied to a
much broader class of diphoton resonances. The ALP Lagrangian that we consider
is

$\mathcal{L}=\frac{1}{2}(\partial_{\mu}a)^{2}-\frac{m_{a}^{2}}{2}a^{2}-\frac{g_{a\gamma\gamma}}{4}aF_{\mu\nu}\tilde{F}^{\mu\nu}\
\ ,$ (3)

and in Fig. 1 we summarize the expected reach on $g_{a\gamma\gamma}$ at Belle
II with $50\text{ ab}^{-1}$ of luminosity. The red line shows the projected
sensitivity of our proposed $e^{+}e^{-}+(\gamma\gamma)_{\text{res}}$ channel
in the ALP parameter space compared with the one of the
$e^{+}e^{-}\to\gamma+(\gamma\gamma)_{\text{res}}$ channel derived in Ref.
Dolan:2017osp . The improvement of the reach in $g_{a\gamma\gamma}$ is
substantial for an ALP lighter than 6 GeV promptly decaying into di-photons.
For higher masses, the ALP searches at Belle II quickly become less sensitive
than existing constraint from ultraperipheral heavy ion collisions at ATLAS
Knapen:2017ebd  and CMS CMS:2018erd . In Sec. 2 we present the details of the
analysis that lead to the sensitivity curve at Belle II of Fig. 1.

Given the ALP life-time

$c\tau_{a}\simeq 1.5\text{ mm}\left[\frac{1\text{
GeV}}{m_{a}}\right]^{3}\left[\frac{5\times 10^{-6}\text{
GeV}^{-1}}{g_{a\gamma\gamma}}\right]^{2}\,,$ (4)

for sufficiently light ALP masses the decay length becomes macroscopic,
ultimately exceeding the length of the Belle II decay volume. For this reason,
the expected sensitivity in Fig. 1 accounts also for the possibility that the
ALP decay happens so far out from the interaction point its production leaves
no trace in the detectors, resulting in a $e^{+}e^{-}+\text{invisible}$ final
state. We studied this signal in our earlier Ref. Acanfora:2023gzr . The
details of how this effect is taken into account in our simulation are given
in App. D. In the intermediate region, a dedicated search for displaced
diphoton resonance could possibly further enhance the sensitivity. Here we do
not distinguish prompt or displaced diphoton resonance and ignore any possible
quality cut on the reconstruction of the diphoton vertex so that any displaced
diphoton event is considered as prompt.111In absence of such dedicated
analysis, our sensitivity estimate suffers a small, but appreciable, drop in
the mass range $m_{a}\simeq 0.1-1\text{ GeV}$. For masses around these values
and for the coupling strength at which we can put a bound the decay length
becomes comparable to the $\mathcal{O}(m)$ size of the detector, resulting in
a loss of rate for both the search presented in this work and for that of our
earlier Ref. Acanfora:2023gzr .

As already noted in Ref. Acanfora:2023gzr , the Belle II geometry does not
allow to fully take advantage of the photon-fusion production because of the
limited acceptance to forward electron and positrons. The current angular
coverage of Belle II can be written in terms of the electron (positron)
pseudorapidity in the center of mass frame as

$\text{Belle~{}II}:\quad|\eta^{*}_{e^{\pm}}|=\left|\log(\tan({\theta^{*}_{e^{\pm}}\over
2}))\right|\leq 1.64\ ,$ (5)

and corresponds to electron and positron with a _minimal_ angle of $22\degree$
with respect to the beam axis.

Since the limited forward acceptance of Belle II cuts away a large fraction of
the fusion production cross-section, it is interesting to explore how the
reach of diphoton resonances could be improved by a hypothetical experiment
with a wider angular coverage at an $e^{+}e^{-}$ facility. The result for this
hypothetical detector, which we call Belle-fwd, are presented in green in Fig.
1 assuming the same luminosity of Belle II. In order to fix a benchmark, we
imagine a pseudo-rapidity coverage in the center of mass frame for both
positrons and electrons of

$\text{Belle-
fwd}:\quad|\eta^{*}_{e^{\pm}}|=\left|\log(\tan({\theta^{*}_{e^{\pm}}\over
2}))\right|\leq 5\ ,$ (6)

which corresponds to a minimal angle with respect to the beam axis of
$0.76\degree$. For reference, this roughly corresponds to the forward
acceptance of the ATLAS calorimeter.

The possibility of extending the forward coverage at Belle II is limited by
the presence of collimating magnets around the beam line which do not allow
this region to be instrumented with detectors. The magnets are required to
reach the expected integrated luminosity of the Belle II experiment Belle-
II:2010dht . Thanks to the impressive boost in the signal rate for an extended
forward acceptance, even with reduced luminosity it is still possible to get
an improvement in the reach compared to the Belle II experiment. This of
course will also reduce the challenge of squeezing the beam to reach the
luminosity targeted at Belle II. A more detailed illustration of how the reach
in this channel depends on the forward acceptance and the luminosity is shown
in Fig. 4 and discussed in Sec. 3.

We conclude in Sec. 4 leaving the technical details of our study for a series
of appendices. In App. A we discuss how we overcome the challenges of
generating a smooth SM background prediction in our binned 4D likelihood
analysis. In particular, we discuss the different possible statistical
treatments of the empty bins in the MC, which result in the uncertainty in the
expected reach in Fig. 1. In App. B we further detail the features of our
likelihood. In App. C we quantify the range of applicability of the well-known
effective photon approximation (EPA) Fermi:1925fq ; Williams:1934ad ;
Budnev:1975poe ; Frixione:1993yw  that is often advocated to carry out
simplified cross-section computations and to develop intuition on the signal
kinematics. App. D discusses the event by event reweighting we performed to
account for the ALP lifetime.

##  2 Expected sensitivity at Belle II

In our signal we deal with $e^{\pm}$ and $\gamma$ detected in Belle II.
According to Belle II detector performance Belle-II:2018jsg  these objects are
well reconstructed if they satisfy the following energy and angular
requirements:

$E^{*}>0.25\text{ GeV},\quad\,\theta^{*}\in[22,158]\degree\,,$ (7)

where the starred quantities are measured in the collision center of mass
frame. In addition, as we want to discuss a resonance in the spectrum of
resolved photons, we require

$\Delta\theta_{\gamma\gamma}>0.048\quad{\text{ or
}}\quad\Delta\phi_{\gamma\gamma}>0.048\ ,$ (8)

where $\Delta\theta_{\gamma\gamma}$ and $\Delta\phi_{\gamma\gamma}$ indicate
the polar and azimuthal angular differences between the two resonant photons.
This requirement ensures that the two photons in the signal final state are
reconstructed as separated photons at Belle II Belle-II:2018jsg .

The photons need also to be separated from the electron to give rise to a
clean final state, thus we require for all $e^{\pm}$-photon pairs

$\Delta\theta_{\gamma e^{\pm}}>0.048\quad{\text{ or }}\quad\Delta\phi_{\gamma
e^{\pm}}>0.048\ .$ (9)

The angular separation and energy requirements of Eqs. 7, 8 and 9 guarantee
the convergence of the perturbative expansion so that both our signal and
background samples are generated with Madgraph5_aMC@NLO v2.7.3 Alwall:2014hca
at leading order in perturbation theory.

###  2.1 Signal vs Background

Two possible ALP production mechanisms at a lepton collider are

$\displaystyle e^{+}e^{-}$ $\displaystyle\to$
$\displaystyle\gamma_{\text{vis}}a\,,$ (10) $\displaystyle e^{+}e^{-}$
$\displaystyle\to$ $\displaystyle e^{+}_{\text{vis}}e^{-}_{\text{vis}}a\,,$
(11)

with $a\to\gamma\gamma$. We have denoted by the subscript $vis$ the final
states that are in the acceptance of the Belle II detector Belle-II:2018jsg
given in Eq. (7). Depending on the lifetime of the ALP, the production
mechanisms above can give rise to either a visible resonant diphoton pair
$a\to(\gamma\gamma)_{\text{res}}$, the processes Eq. (10) and Eq. (11), or to
missing energy and momentum $a\to\text{invisible}$.

Figure 2: Feynman diagram for the $e^{+}e^{-}\to e^{+}e^{-}a$ production.

The signal in Eq. (10) is the well studied “ALP-strahlung” production (see
Ref. Dolan:2017osp ), while the process in Eq. (11) gives rise to two
different processes: $i$) the “ALP-Dalitz” process, given by the ALP-strahlung
Eq. (10) with a photon conversion into $e^{+}e^{-}$; $ii$) the “photon-fusion”
into ALP, given by a photon line exchanged in the $t$-channel that radiates
the ALP depicted in Fig. 2. Similarly to what we found in Acanfora:2023gzr
for the invisible ALP decay, the sensitivity for an ALP decaying into a
diphoton pair at Belle II is dominated by the photon-fusion production. This
dominance would be further enhanced in an experiment with extended forward
detector coverage.

For each putative $m_{a}$ we require the two final state photons to have an
invariant mass in a range

$m_{\gamma\gamma}-m_{a}\in[-3,1.5]\cdot\sigma_{m_{a}}\,,$ (12)

where $\sigma_{m_{a}}$ denotes the expected resolution on a di-photon
resonance Belle-II:2018jsg , which is dominated by the uncertainty on the
photon energy depositions in the ECAL. This requirement removes large
fractions of the backgrounds for all $m_{a}$ of interest and entirely kills
the Dalitz contribution to the signal, resulting in a very small loss of
signal rate.

The most relevant background to the photon fusion production in Eq. (11) is
given by radiative corrections to Bhabha scattering

$e^{+}e^{-}\to e^{+}e^{-}\gamma\gamma\,.$ (13)

A large fraction of the cross-section for this process can be understood
qualitatively as $e^{+}e^{-}\to e^{+}e^{-}$ dressed by two $e^{\pm}\to
e^{\pm}\gamma$ splittings. As a consequence of the soft and collinear
enhancement of the splittings, a distinguishing feature of this background
process is the tendency of each photon to fly collinear to the fermion that
has emitted it, as well as being preferred to be soft compared to the parent
$e^{\pm}$. As a further consequence of the collinear preference of the photons
emerging from Bhabha scattering, we observe that a large fraction of the
background events tends to have photons flying in opposite directions. This is
expected when the two photons are arising from different fermionic lines in
the hard scattering.

Figure 3: Distributions of signal events for the ALP production in photon
fusion in Eq. (11) (in blue) and background events for the radiative Bhabha in
Eq. (13) (in red). A fixed diphoton invariant mass of
$m_{\gamma\gamma}=1\text{ GeV}$ within the range of Eq. (12) is imposed, as
well as the acceptance and isolation cuts at Belle II given in Eqs. 7, 8 and
9. Left: distributions as a function of the minimal angle between the
positron/electron direction and the photon and the total energy of the photon
pair, Right: distributions as a function of the minimal angle between the
positron/electron direction and the photon and the angular separation between
the two photons.

In Fig. 3 the distributions of signal and background at fixed invariant mass
of 1 GeV are shown. As discussed above, the background collinear enhancement
will favor a small angle between the photon and the corresponding fermionic
line (either the electron or the positron). As a consequence, in both panels
in Fig. 3 we can see how the background event distribution favors a small
_minimal_ angle between the photon and electron (or positron) in the final
state.

Moreover, for a small enough diphoton invariant mass, the signal events
passing the Belle II acceptance cut typically feature a boosted ALP resonance
(as already shown in Ref. Acanfora:2023gzr  but see also Fig. 8 in App. C),
which results in an energetic photon pair with relatively small opening angle
between the photons. Conversely, the enhancement of the splitting function for
soft radiation prefers the photons energies to be small and the angle between
them could easily be large to make the required invariant mass if for instance
they are split in opposite directions. The separation between signal and
background becomes less sharp at higher invariant masses, when the boost of
the signal resonance diminishes. This feature qualitatively explains the
degradation of the sensitivity at Belle II apparent in Fig. 1.

Before leaving this section, we comment on further backgrounds beyond the one
in Eq. (13). One possibility is the rare mis-reconstruction of hadrons as
electrons and fake photons, e.g. from mis-identification of neutral mesons,
e.g. a $\pi^{0}$ detected as a photon. The efficiency-vs-purity performance of
the Belle II experiment on particle-ID is a subject of active development (see
e.g. Ref. Sebastian:1821 ). A realistic estimate of these backgrounds is
challenging, especially when multiple hadrons arise from a copious two-body
hard scatterings such as $e^{+}e^{-}\to q\bar{q}$. The final background rate
depends mainly on the detector operation, e.g. isolation requirements, and
should be performed by the experimental collaboration.

###  2.2 Sensitivity

For an optimal signal vs background separation we pursue a full four
dimensional likelihood, which can take full advantage of their different
kinematic structure. In order to efficiently tessellate the four dimensional
phase space we define, in analogy with the usual Mandelstam variables, the
Lorentz invariants $t_{\pm}$ and $s_{\pm}$ as

$\begin{split}t_{+}&=\left(p_{e_{in}^{+}}-p_{e_{out}^{+}}\right)^{2}\,,\\\
t_{-}&=\left(p_{e_{in}^{-}}-p_{e_{out}^{-}}\right)^{2}\,,\\\
s_{-}&=\left(p_{e_{in}^{+}}+p_{e_{in}^{-}}-p_{e_{out}^{+}}\right)^{2}\,,\\\
s_{+}&=\left(p_{e_{in}^{+}}+p_{e_{in}^{-}}-p_{e_{out}^{-}}\right)^{2}\,,\end{split}$
(14)

where the positron, electron, and photon momenta are defined in Fig. 2. A
conversion between the Lorentz invariant variables and the measured energy and
angles of the positron, electron, and photons can be easily derived.
Neglecting the electron mass, we can write

$t_{\pm}\simeq-\sqrt{s}E_{e_{out}^{\pm}}\left(1\pm\cos\theta_{e_{out}^{\pm}}\right)\quad,\quad
s_{\pm}\simeq\sqrt{s}\left(\sqrt{s}-2E_{e_{out}^{\pm}}\right)\ ,$ (15)

where $\theta_{e_{out}^{\pm}}$ is the polar angle of the positron/electron
with respect to the beam axis. The above expressions show that $t_{\pm}\to 0$
controls the collinear singularity of the photon fusion production at
$\theta_{e_{out}^{\pm}}\to 0$. Using energy conservation, the sum of the other
two Mandelstam variables can be written in terms of the energy of the
resonance $s_{+}+s_{-}\simeq 2E_{a}\sqrt{s}$, while their difference gives the
asymmetry between the positron and the electron energy.

Starting from the energy and polar angle resolution of Belle II Belle-
II:2018jsg  we derive approximate resolutions for the measurement of $t_{\pm}$
and $s_{\pm}$. Each of these quantities contains energies and directions of
photons and electrons. The energies are the least well measured, so we can
take

$\frac{\delta\mu_{i}}{\mu_{i}}\simeq 2\frac{\delta E}{E}\simeq
4\%\quad\text{for}\quad\mu_{i}=\\{s_{+},s_{-},t_{+},t_{-}\\}\ .$ (16)

Based on the predicted differential rates for the SM background and for the
signal, we construct a log-likelihood for 50 ab-1 integrated luminosity at
Belle II

$\Lambda=-2\sum_{i,j}\ln\frac{L(S_{i,j},B_{i,j})}{L(0,B_{i,j})}\,,$ (17)

where $i$ and $j$ run on the bins of the 4D space spanned by $s_{+}$, $s_{-}$,
$t_{+}$, $t_{-}$ and $S$ and $B$ indicate respectively the expected number of
signal or background events in each bin. In Eq. (17) $L(S_{i,j},B_{i,j})$ is
the Poisson factor computed in each bin defined as follows

$L(S,B)=\frac{\left(S+B\right)^{B}}{B!}e^{-(S+B)}\,.$ (18)

The sensitivity shown in Fig. 1 corresponds to 95% C.L. and it is obtained by
requiring $\Lambda<4$.

Given the huge hierarchy between the signal rate and the background rate, our
sensitivity benefits from regions of phase space where the background is
strongly suppressed. To best exploit these features, we perform a full 4D
likelihood analysis. This likelihood is challenging to construct from the
viewpoint of the MC generation, but the improvement in the sensitivity is
worth the effort.222We note in passing that the expected sensitivity for a 3D
likelihood obtained by collapsing our 4D one in any of the four directions is
reduced by roughly one order of magnitude. This loss of sensitivity confirms
that the full phase space information is needed to distinguish the photon
fusion process from the background. For computational reasons we then
tessellate the 4D phase space with different bin sizes depending on the ALP
mass. Specifically for $m_{a}<1\text{ GeV}$ we take
$\delta\mu_{i}/\mu_{i}=15\%$, for $1\text{ GeV}\leq m_{a}<3\text{ GeV}$ we
take $\delta\mu_{i}/\mu_{i}=8\%$ and for $m_{a}\geq 3\text{ GeV}$ we take the
realistic Belle II resolution estimated in Eq. (16). In this sense, the
expected reach of Fig. 1 at low masses should be considered conservative and
could be improved by reducing the bin size further.

Even with coarse grained bins, our expected reach has a residual dependence on
the statistical treatment of the phase space cells with zero expected
background within our MC uncertainty. Evaluating this uncertainty ultimately
results in the uncertainty band in the expected reach shown in Fig. 1. A
detailed description of our methods to generate a smooth background sample in
the full 4D phase space and assigning a conservative uncertainty to it is
given in Appendix A.

We also checked the robustness of our result against possible systematic
uncertainties. These are particularly dangerous as one approaches a high
statistics regime and small couplings can be potentially probed. To compute
the loss of sensitivity due to systematics we computed the log-likelihood
keeping only bins for which the expected $S/B$ was larger than a certain
threshold $\epsilon_{\rm{sys}}$. For $\epsilon_{\rm{sys}}<10\%$ the expected
sensitivity does not change appreciably, especially at low masses. This is
because the likelihood is dominated by bins with very suppressed background
yield, as discussed in detail in Appendix B.

The variable $t_{+},(t_{-})$ in Eq. (14) is related to the degree of
collinearity of the initial and final $e^{+},(e^{-})$, which in turn can be
seen as the collinearity of the virtual photons and the initial
$e^{+},(e^{-})$. These quantities are typical markers of $t$-channel
scattering, which dominates the photon-fusion production of the ALP. As it can
be appreciated from Fig. 5 in Appendix B, the Belle II acceptance cuts off
this singularity reducing the signal acceptance. In the next section we show
how a hypothetical experiment with extended forward acceptance would allow
instead to explore kinematic configurations where both background and signal
are closer to the collinear singularity and can be effectively distinguished
by exploiting the full four dimensional phase space information.

Figure 4: The green contours show the expected reach for
$\log_{10}(g_{a\gamma\gamma})$ in the $e^{+}e^{-}\to
e^{+}e^{-}(a\to\gamma\gamma)$ channel as a function of the luminosity and the
final $e^{\pm}$ acceptance of a hypothetical lepton collider at
$\sqrt{s}=10.58~{}\text{ GeV}$, given $m_{a}=3$ GeV. The red and green stars
are the benchmarks for the Belle II and Belle-fwd reach in Fig. (1). As a
reference we show as a orange dashed line the Belle II reach for the
$e^{+}e^{-}\to\gamma(\gamma\gamma)_{\text{res}}$ channel. The red dashed line
is the isoline of the value of the Belle II reach.

##  3 Expected sensitivity of a forward flavor factory

In this section we entertain the possibility that a future experiment with
extended forward coverage will be available at an $e^{+}e^{-}$ facility
similar to SuperKEKB. Given that the expected luminosity of such a machine is
likely to be correlated with its forward instrumentation, we want to study how
the expected reach of the search discussed in Sec. 2 depends on the forward
acceptance, $\eta^{*}_{e^{\pm}}$, and the luminosity, $\mathcal{L}$. We do not
make a concrete technical proposal here on how to construct and operate this
extended detector. We remark however that this is an active area of research
for future collider detectors studies Accettura:2023ked ; Ruhdorfer:2023uea .

Our results are shown in Fig. 4 for benchmark ALP mass of 3 GeV and
$\sqrt{s}=10.58\text{ GeV}$, where the red and the green stars correspond to
the Belle II reach and the Belle-fwd reach of Fig. 1. The reach in
$g_{a\gamma\gamma}$ defined in Eq. (3) decreases roughly like
$\mathcal{L}^{1/4}$ as expected for a search that is not free from background.
Conversely, increasing the forward acceptance leads to a steep rise of the
sensitivity which decreases only at high forward acceptances. We find that,
thanks to the impressive gain in reach showed in Fig. 1, even an experiment
collecting 100 times less luminosity than Belle II, profiting of an extended
forward acceptance down to $\eta_{e^{\pm}}^{*}<5$, could lead to a substantial
improvement in the reach compared to Belle II.

##  4 Conclusion

In this work we developed a new search for diphoton resonances at Belle II
based on the photon fusion production in Eq. (1). By exploiting a full
tessellation of the four dimensional phase space of this $2\to 3$ process we
performed a detailed study of how this can be separated from the background
given by the radiative Bhabha process in Eq. (13). The result of this analysis
is summarized in Fig. 1 which shows how this channel can lead to a substantial
improvement in the reach compared to the “ALP-strahlung” production studied in
Ref. Dolan:2017osp  at least for ALP masses below 6 GeV. Above 6 GeV the Belle
II reach quickly becomes less sensitive than existing constraints from ultra-
peripheral heavy ions collisions Knapen:2017ebd ; CMS:2018erd . Along the way,
we illustrate in Appendix A a method to associate a conservative uncertainty
to the signal reach in analyses plagued by challenges in the MC sampling like
ours.

In addition, we explored how the expected reach for a diphoton resonance in
the $e^{+}e^{-}+\gamma\gamma$ channel can be further improved at a
hypothetical experiment with extended forward coverage compared to Belle II.
The results of Fig. 1 together with the ones in Fig. 4 show that with an
extended forward coverage down to $|\eta_{e^{\pm}}^{*}|<5$ the expected Belle
II reach could be improved even with 100 time less luminosity.

## Acknowledgments

We thank Torben Ferber, Enrico Graziani, Antonio Passeri and Laura Zani for
useful discussions about Belle II. RF and DR are supported in part by the
European Union - Next Generation EU through the MUR PRIN2022 Grant
n.202289JEW4.

##  Appendix A Generation of the background sample for the 4D likelihood
analysis

The generation of 4-fold differential cross-section requires intense
computational resources. A relatively gross-grained sampling of each dimension
with about 100 bins leads to a rather large set of $10^{8}$ bins in a single
histogram. Obtaining the value of the cross-section with precision $1/p$ in
each bin requires $p^{2}$ samplings of the differential cross-section, which,
even for leading order matrix elements, takes some substantial time.

Significant speed-up of the generation of the information necessary to
characterize our background in a 4D phase-space can be obtained i)
parallelizing the generation of multiple MC samples, ii) using a clever re-
weighting of the events, iii) applying suitable biases to the samples, and iv)
combining the information from the several runs.

In our calculation we have applied a bias to the differential cross-section as
to overpopulate tails at large $s_{+}$ and large $s_{-}$, and to increase the
population of the area of phase-space around $t_{+}\simeq t_{-}\simeq 0$.
These regions are most difficult to describe precisely with Madgraph5_aMC@NLO
both because they are tails and because they are rather narrow region with
respect to the whole phase-space of the process. In addition, $t_{+}\simeq
t_{-}\simeq 0$ is a boundary of the physical phase-space, which poses extra
difficulties. The accurate description of these regions of phase-space is of
utmost importance for our results as they give a large discrimination of
signal over background within our likelihood in Eq. (17).

Particularly challenging is the generation of background events for $m_{a}\leq
1$ GeV at Belle II, because the kinematics of the signal is such that the
background suppression in specific regions of phase space is larger than for
higher masses. In what follows we describe our MC strategy in this light mass
range and mention the simplification at higher masses at the end.

For $m_{a}\leq 1$ GeV, we produce ${\cal O}(10)$ MC samples with
Madgraph5_aMC@NLO with bias and cuts aimed towards several different
combinations of small and large $s_{+}+s_{-}$ and similarly for $t_{+},t_{-}$.
From each sample we obtain a weight, that corresponds to the cross-section for
each given cell of the 4D phase-space. We denote this weight as $w_{c,s}$,
where $c$ denotes the cell of the phase-space and $s$ the sample we are
dealing with. Each such weight can be endowed of an estimated uncertainty,
which originates from the finiteness of the MC samples that we generated. We
estimate an uncertainty $\delta w_{c,s}=w_{c,s}/\sqrt{n_{c,s}}$, where
$n_{c,s}$ is the number of weighted MC events in the cell $c$ from the sample
$s$.

We then combine the weights of a given cell $c$ as if they were the results of
different experiments into a combined value

$w_{c,C}=\frac{\sum_{s}w_{c,s}\cdot\xi_{c,s}}{\sum_{s}\xi_{c,s}}\,,$ (19)

where $\xi_{c,s}=\delta w_{c,s}^{-2}$, as in the usual combination of two or
several measurements described by an underlying Gaussian probability density
function. Linear uncertainty propagation is applied to this formula, finding
the standard result for the uncertainty on the combined value

$\delta w_{c,C}=\left(\sum_{s}\delta w_{c,s}^{-2}\right)^{-1/2}\,.$ (20)

The samples $s$ included in the combinations are chosen so that only
measurements that are compatible with each other are used. In particular, we
combine only samples that have non-zero weight on the cell $c$ and discard the
samples that returned zero weight, either because $c$ was entirely out of the
region of phase-space of the MC run that produced sample $s$, or because the
bias in that run was put to highlight other tails of the phase-space. The runs
with zero weight are discarded on the ground that they correspond to
incompatible measurements, and therefore cannot be combined with the other
group. In fact, the information from the samples returning zero weight would
imply zero cross-section on that cell, which is clearly incompatible with the
information from a dedicated sample that has found non-zero cross-section in
that cell.

When a cell is not populated by any of our background samples simulations, we
conservatively assume that the cross-section in that cell is non-zero and that
each sample has given us an upper-bound on the cross-section in that cell.
Each sample contributes an upper-bound on the weight $w_{c,s}<\Delta w_{s}$,
where $\Delta w_{s}$ is the minimum weight in the sample $s$. The combined
upper-bound then is

$w_{c,C}<\Delta w_{c,C}\quad\text{ for }\quad\Delta
w_{c,C}=\left(\sum_{s}\Delta w_{s}^{-2}\right)^{-1/2}\,.$ (21)

From the above procedure we obtain a combined histogram in 4D that has either
a measurement or an upper bound on the cross-section in each cell of phase-
space. In the following we describe a method to assess the impact of the
choice in the treatment of these cells for which we have obtained an upper
bound.

In order to give an estimate of the error on the bounds on $g_{a\gamma\gamma}$
due to the finiteness of our MC samples, we produce replicas of our combined
histogram. In each replica we fluctuate the number of expected events in the
cell $c$ based on the underlying uncertainty on the combined weight $\delta
w_{c,C}$. For the bins with non-zero background events the number of the
expected events follows a truncated normal distribution where the gaussian
mean and standard deviation are taken from Eq. (19) and Eq. (20).

For the empty bins we only have the upper bound in Eq. (21), which leaves us
the freedom of picking the probability distribution in the bin. Variations of
the underlying probability distribution for the bins gives the uncertainty in
the expected reach in Fig. (1). We identify three scenarios:

  * •

The weakest constraints come from assuming that the distribution is a Dirac
$\delta$ function centered at the upper bound in Eq. (21). This gives the
maximal number of expected background events compatible with having observed
zero events in our Monte Carlo.

  * •

The strongest constraints are obtained assuming a discrete Poisson
distribution with mean set at the upper bound in Eq. (21). The discrete nature
of this distribution gives a sizable probability to have zero expected
background, thus boosting the likelihood in Eq. (17).

  * •

Intermediate bounds between these two options are obtained by picking a
continuous uniform probability distribution from zero to the upper bound in
Eq. (21). The central line in the bands in Fig. (1) corresponds to this
choice, which we consider as the most reasonable estimate of the attainable
constraints.

The background generation for $m_{a}\geq 1\text{ GeV}$ is analogous to the one
explained above, with the technical simplification that a single MC sample was
sufficient to obtain a smooth background sample. Running multiple samples
would have allowed us to reduce the bin size of the likelihood for $1\text{
GeV}\leq m_{a}<3\text{ GeV}$ and in general reduce the uncertainty in our
projected sensitivity for all the ALP masses. The projected sensitivity and
its uncertainty shown in Fig. 1 should then be taken as a first conservative
assessment of the discriminating power of the $e^{+}e^{-}\to
e^{+}e^{-}(\gamma\gamma)_{\text{res}}$ channel.

Figure 5: Distribution of log-likelihood defined in Eq. (17) in the
$(t_{-},t_{+})$ plane for $m_{a}=3\text{ GeV}$. The 4D likelihood is maximized
in the $(s_{-},s_{+})$ plane to be projected on a 2D plane. The signal
strength is fixed to the expected reach of Belle-fwd with 50 $\text{ab}^{-1}$
of data which is $g_{a\gamma\gamma}=5.5\cdot 10^{-6}$ GeV-1 as shown in Fig.
1. The white region is kinematically forbidden because of the detector
acceptance. Left Belle II likelihood with acceptance defined in Eq. (5) Right
Belle-fwd likelihood with acceptance defined in Eq. (6).

##  Appendix B Features of the 4D likelihood

Figure 6: Likelihood distributions for $m_{a}=3\text{ GeV}$ in the plane
$(S/\sqrt{B},S/B)$ where $S$ is the expected number of signal events and $B$
the expected number of background events. The three panels correspond to the
three different assumption on the distribution of the background events in the
empty bins as illustrated in Appendix A. Left: The events in the empty bins
are distributed as a Poisson distribution with mean set to the upper bound of
the cross-section derived in our MC. The likelihood correspond to the best
expected sensitivity in Fig. 1 which is $g_{a\gamma\gamma}=5.49\cdot
10^{-6}\text{ GeV}^{-1}$ Center: The events on the empty bins are distributed
as a uniform distribution with maximum given by the upper bound derived in our
MC. This treatment correspond to the red line in Fig. 1 which gives
$g_{a\gamma\gamma}=2.88\cdot 10^{-5}\text{ GeV}^{-1}$. Right: The events on
the empty bins are distributed as a delta-function centered at the upper bound
derived in our MC. This assumption gives the most conservative expected reach
in Fig. 1 which is $g_{a\gamma\gamma}=5.89\cdot 10^{-5}\text{ GeV}^{-1}$.

We now comment on the features of the 4D likelihood which result in the
sensitivity at Belle II and at Belle-fwd in Fig. 1. Firstly, we show in Fig. 5
the distribution of our likelihood as defined in Eq. (17) in the
$(t_{-},t_{+})$ plane. For every $t_{-},t_{+}$ pair, we pick the maximal
likelihood value in the $(s_{-},s_{+})$ plane. In the left plot we see the
distribution at Belle II while in the right one at a hypothetical experiment
with extended forward coverage at an $e^{+}e^{-}$ facility and the same
luminosity of Belle II.

Figure 7: Dependence of the expected reach on the systematic uncertainties on
the background rate. The gray shaded region shows existing constraints as in
Fig. 1. Left: The red band shows the expected reach of
$e^{+}e^{-}+(\gamma\gamma)_{\text{res}}$ as in Fig. 1, and in magenta the same
reach with the requirement $S/B>0.1$, as detailed in Appendix B. The arrows of
the same color show how the uncertainty of our expected reach gets larger with
larger systematic uncertainties. Right: The green band corresponds to the
reach of Belle-fwd as defined in Eq. (6) and detailed in Sec. 3; the orange
band is the same reach with the requirement $S/B>0.1$.

The shape of the likelihood distribution can be explained by noticing that the
absolute values of the variables $t_{\pm}$ as defined in Eq. (15) are bounded
from below by a quantity that depends on the detector acceptance which
determines the minimal polar angle of the electron/positron that one can
measure, $\theta_{e_{out}^{\pm}}^{\text{min}}$. In formulas we can write the
bound for tip of the wedge with $t_{+}~{}\sim t_{-}$ as

$|t_{\pm}|>\frac{s}{4}(1-\cos\theta_{e_{out}^{\pm}}^{\text{min}})\ ,$ (22)

where $\theta_{e_{out}^{\pm}}^{\text{min}}$ can be extracted from Eq, (5) for
Belle II and from Eq. (6) for Belle-fwd. The comparison between the two panels
shows clearly how the extended forward coverage opens up a new region of phase
space where the likelihood is considerably larger, hence boosting the
sensitivity to the ALP.

We now illustrate the dependence of the likelihood on our treatment of the
empty bins illustrated in Appendix A. Calling $S$ the expected signal events
and $B$ the expected background events, we show in Fig. 6 the distribution of
the likelihood in the plane spanned by $S/\sqrt{B}$, which roughly controls
the statistical weight, and $S/B$, which is a good proxy for the sensitivity
of our result on the systematic uncertainties in the background rate. The
three panel illustrate the three different treatments of the empty bins, and
the coupling strength of the ALP in each panel is fixed at its corresponding
expected sensitivity at Belle II.

In the left panel we show the likelihood assuming that the probability
distribution of the empty bins is Poissonian and has its mean set at the upper
bound of the background cross-section obtained from our MC exploration. The
discrete nature of the likelihood makes most of the likelihood concentrated on
the tail of the Poisson distribution with zero background events which is
shown as a red square dot in the plot. The center panel shows the likelihood
obtained under the assumption that the events in the empty bins are uniformly
distributed with maximal value given by the upper bound derived in our MC.
Here we see that the likelihood distribution becomes nicely continuous in the
$(S/\sqrt{B},S/B)$ plane. In the right panel we show the likelihood
distribution where we assume the background rate in the empty bins to be fixed
at the lower bound derived from the MC. This assumption gives the most
conservative expected sensitivity in Fig. 1. For both the center and right
panel, we observe that a large fraction of the likelihood lives at large value
of $S/B$, anticipating the low dependence of our expected reach on systematic
uncertainties on the background.

The dependence of the expected reach on the systematic uncertainties is shown
in Fig. 7. To quantify this in our likelihood analysis we restrict our
attention to the bins with $S/B$ larger than a given value. The lower bound on
$S/B$ can be taken as a proxy for the expected systematic uncertainty. For
$S/B>0.1$ we see that the uncertainty in our expected reach is sensibly
increasing at Belle II (left panel) while it is left almost unchanged at
Belle-fwd (right panel). The sensitivity of the likelihood to the cut in $S/B$
depends very much on the statystical treatment of the empty bins as can bee
inferred from Fig. 6. In particular, for the we assume the events in the empty
bins to be distributed as a Poissonian the likelihood is dominated by bins
with $B=0$ and therefore essentially insensitive to systematic uncertainties.
This is confirmed by looking at the left panel of Fig. 6. Therefore,
accounting for systematic uncertainties does not shift our best sensitivity.
Conversely, systematic uncertainties affect our most conservative sensitivity
which is obtained by taking the upper bound on the cross-section derived in
our MC as the number of expected background events. At Belle II we see how
having systematic uncertainties at 10% level shifts the reach by a roughly a
factor of five on the coupling. This is confirmed by looking at the right
panel of Fig. 6.

##  Appendix C Effective photon approximation at flavor factories

The photon fusion process that we propose as a new channel to search for ALP
production in Belle II is traditionally linked to the equivalent vector boson
approximation Fermi:1925fq ; Williams:1934ad ; Budnev:1975poe ;
Frixione:1993yw . In the context of LHC physics this approximation provides a
simple computational tool to estimate the total cross-section of fusion
processes, see e.g. Ref. Dawson:1984gx ; Cahn:1983ip  for early applications
in the context of Higgs boson production in hadron colliders. While the
computational aspect is no longer as pressing as it used to be, the equivalent
boson picture can in principle be used to guide the isolation of the fusion
signal from the backgrounds.

Figure 8: Mean and median $p_{T,a}/m_{a}$ for Eq. 1 including phase-space up
to a minimal angle $\theta_{0}$ (maximal pseudo-rapidity $|\eta_{0}|$).

The basic idea of the approximation is to consider as on-shell the gauge boson
field that enters in a scattering, as long as its virtuality can be neglected
in comparison with the momentum exchanged in the scattering in which the boson
participates. In our case the photon, in order to be approximately on-shell,
must have negligible virtuality with respect to the mass of the $a$ that
emerges from the photon-fusion. The physical foundation of this approximation
relies on the fact that the process of photon fusion happens in a time of
order $1/m_{a}$, which is too short to probe the virtuality of the photon,
$q_{i}^{2}$ as defined in Fig. 2 if $q_{i}^{2}\ll m_{a}^{2}$ Borel:2012by .
The virtuality of the photon can only be “detected” in experiments sensitive
to much longer time-scales.

As in our kinematics the photon virtualities are linked to the final state
electron and positron transverse momenta

$q_{1}^{2}\simeq p^{2}_{T,e_{\text{out}}^{-}},\quad q_{2}^{2}\simeq
p^{2}_{T,e_{\text{out}}^{+}}\quad\text{ and }\quad
p_{T,e_{\text{out}}^{+}}+p_{T,e_{\text{out}}^{-}}=p_{T,a}\ ,$ (23)

it is instructive to look at the typical transverse momentum of the ALP
divided by its mass. In Fig. 8 we show the median and the mean $p_{T,a}/m_{a}$
for $m_{a}=1\text{ GeV }$ $m_{a}=8\text{ GeV }$ including phase-space up to a
minimal angle $\theta_{0}$ (maximal rapidity $|\eta_{0}|$) from the beam axis.
We observe that for 1 GeV ALP mass, and for similar ALP masses compared to the
center of mass energy, $p_{T,a}/m_{a}\sim 1$, which indicates that such light
ALP is not accurately described by the equivalent photon approximation. This
is especially true for $\theta_{0}\simeq 1\degree$ and $22\degree$, the
benchmarks for the Belle-fwd hypothetical detector and the actual Belle II
detector respectively.

To corroborate this expectation, we show in the left panel of Fig. 9 the exact
calculation of the fusion process Eq. 1 and the EPA prediction as computed at
its lowest order following the definition of Ref. Budnev:1975poe . The EPA
prediction is further equipped with an error band in light (dark) yellow,
which corresponds to an uncertainty given by the mean (median) $p_{T,a}/m_{a}$
shown in Fig. 8. The agreement between the EPA and the exact computation
degrades as $\theta_{0}\gtrsim 0.01\degree$ ($|\eta_{0}|\lesssim 10$), with an
order of magnitude disagreement expected for Belle-fwd and even worse for
Belle II. In the right panel of Fig. 9 we show the same quantities for a heavy
ALP mass of 8 GeV. In this case a much better agreement between EPA and the
exact computation is observed, in full compliance with the smallness of the
ratio $p_{T,a}/m_{a}\lesssim 0.1$ over the whole range of $\theta_{0}$.

From these analyses we conclude that the EPA works well in the cases in which
the observed final states corresponds to small $t$-channel photon
virtualities, but it is definitively possible to reach situations in which EPA
does not provide a good estimate of the ALP total production cross-section. In
these situations, also the kinematic picture of a signal dominated by small
$p_{T,a}$ can be flawed. In fact, for ALPs mass 1 GeV in Fig. 3 we have
observed a strong preference for large $E_{\gamma_{1}}+E_{\gamma_{2}}\gg
m_{a}$, in spite of the EPA expectation of an ALP produced nearly at rest.

Due to the difficulty to build forward detector in machines capable of
reaching high instantaneous luminosity, such as flavor factories, the EPA
should be used with great care to compute cross-sections at these machines.
Similarly, the physics intuition that follows from the validity of the EPA can
be inaccurate for experiments on these machines, due to the absence of
coverage in the forward phase-space.

Figure 9: Comparison of the EPA prediction and the exact cross-section for the
process Eq. 1, in the case of $m_{a}=1$ GeV (left) and $m_{a}=8$ GeV (right).
The black line is the cross-section for the exact computation of the process
Eq. 1 as a function of the lepton acceptance (see Eq. 5), expressed either as
a polar angle measured in degrees $\theta_{0}$ or as a pseudo-rapidity
$\eta_{0}$. The light yellow and the dark yellow error bands are the EPA with
mean and median $p_{T,a}/m_{a}$ relative error, as detailed in the text. The
vertical red and green lines correspond to the angular acceptance at Belle II
and Belle-fwd respectively. For reference we also show the
$\gamma+(\gamma\gamma)_{\text{res}}$ cross-section (which is invariant with
respect to $\theta_{0}$) as an orange dashed line.

##  Appendix D Influence of ALP lifetime

For small enough ALP masses, the decay length becomes so large that most of
the time the ALP decay happens far away from the interaction point, possibly
even outside the detector. In order to account for this effect, we reweigh
every event by including the probability of the ALP to decay outside the
detector length. The decay probability before some distance $L$ can be written
as

$P_{\text{decay}}=1-e^{-\frac{L}{\gamma_{a}\tau_{a}}}=1-e^{-\frac{L\Gamma_{a}m_{a}}{E_{a}}}\
,$ (24)

where $m_{a}$ and $E_{a}$ are respectively the ALP mass and energy and

$\Gamma_{a}=\frac{g_{a\gamma\gamma}^{2}m_{a}^{3}}{64\pi}$ (25)

its decay width. In order to make the most conservative estimate, we will
choose $L$ to be the shortest distance between the interaction point and the
ECL barrel which is $L=1.25\text{ m}$ at Belle II Belle-II:2018jsg . Each
event will then have a weight which is a function of $g_{a\gamma\gamma}$ and
also depends on the ALP mass and energy. We can then include this effect in
the likelihood of Eq. (17). The effects of the ALP lifetime, as expected, are
visible only for $m_{a}\lesssim 0.5\text{ GeV}$; for those masses there is
effectively a reduction of the signal rate. As shown in Fig. 1, however, for
ALP masses below $1\text{ GeV}$, invisible searches become quickly more
competitive, filling the gap of the visible ones.

Interestingly, the loss of promptly decaying events is less rapid at Belle-fwd
(defined in Eq. (6)) with respect to Belle II (see Eq. (5)). This can be
explained by computing the average boost of the particle as a function of the
opening angle of the detector. For a more forward detector like Belle-fwd the
average boost of the ALP will be smaller, making it easier for the ALP to
decay inside the detector.

## References

  * (1) CMS collaboration, S. Chatrchyan et al., _Observation of a New Boson at a Mass of 125 GeV with the CMS Experiment at the LHC_ , _Phys. Lett. B_ 716 (2012) 30–61, [1207.7235].
  * (2) ATLAS collaboration, G. Aad et al., _Observation of a new particle in the search for the Standard Model Higgs boson with the ATLAS detector at the LHC_ , _Phys. Lett. B_ 716 (2012) 1–29, [1207.7214].
  * (3) S. Knapen, T. Lin, H. K. Lou and T. Melia, _Searching for Axionlike Particles with Ultraperipheral Heavy-Ion Collisions_ , _Phys. Rev. Lett._ 118 (2017) 171801, [1607.06083].
  * (4) S. Knapen, T. Lin, H. K. Lou and T. Melia, _LHC limits on axion-like particles from heavy-ion collisions_ , _CERN Proc._ 1 (2018) 65, [1709.07110].
  * (5) A. Mariotti, D. Redigolo, F. Sala and K. Tobioka, _New LHC bound on low-mass diphoton resonances_ , _Phys. Lett. B_ 783 (2018) 13–18, [1710.01743].
  * (6) X. Cid Vidal, A. Mariotti, D. Redigolo, F. Sala and K. Tobioka, _New Axion Searches at Flavor Factories_ , _JHEP_ 01 (2019) 113, [1810.09452].
  * (7) D. Aloni, C. Fanelli, Y. Soreq and M. Williams, _Photoproduction of Axionlike Particles_ , _Phys. Rev. Lett._ 123 (2019) 071801, [1903.03586].
  * (8) D. Wang, L. Wu, J. M. Yang and M. Zhang, _Photon-jet events as a probe of axionlike particles at the LHC_ , _Phys. Rev. D_ 104 (2021) 095016, [2102.01532].
  * (9) S. Knapen, S. Kumar and D. Redigolo, _Searching for axionlike particles with data scouting at ATLAS and CMS_ , _Phys. Rev. D_ 105 (2022) 115012, [2112.07720].
  * (10) ATLAS collaboration, G. Aad et al., _Search for boosted diphoton resonances in the 10 to 70 GeV mass range using 138/fb of 13 TeV $pp$ collisions with the ATLAS detector_, _JHEP_ 07 (2023) 155, [2211.04172].
  * (11) G. Alonso-Álvarez, J. Jaeckel and D. D. Lopes, _Tracking axion-like particles at the LHC_ , 2302.12262.
  * (12) CMS collaboration, A. M. Sirunyan et al., _Evidence for light-by-light scattering and searches for axion-like particles in ultraperipheral PbPb collisions at $\sqrt{s_{\mathrm{NN}}}=$ 5.02 TeV_, _Phys. Lett. B_ 797 (2019) 134826, [1810.04602].
  * (13) R. Balkin, O. Hen, W. Li, H. Liu, T. Ma, Y. Soreq et al., _Probing axion-like particles at the Electron-Ion Collider_ , 2310.08827.
  * (14) M. J. Dolan, T. Ferber, C. Hearty, F. Kahlhoefer and K. Schmidt-Hoberg, _Revised constraints and Belle II sensitivity for visible and invisible axion-like particles_ , _JHEP_ 12 (2017) 094, [1709.00009].
  * (15) F. Acanfora, R. Franceschini, A. Mastroddi and D. Redigolo, _Fusing photons into nothing, a new search for invisible ALPs and Dark Matter at Belle II_ , 2307.06369.
  * (16) T. Ferber, A. Filimonova, R. Schäfer and S. Westhoff, _Displaced or invisible? ALPs from B decays at Belle II_ , _JHEP_ 04 (2023) 131, [2201.06580].
  * (17) E. M. Riordan et al., _A Search for Short Lived Axions in an Electron Beam Dump Experiment_ , _Phys. Rev. Lett._ 59 (1987) 755.
  * (18) B. Döbrich, J. Jaeckel, F. Kahlhoefer, A. Ringwald and K. Schmidt-Hoberg, _ALPtraum: ALP production in proton beam dump experiments_ , _JHEP_ 02 (2016) 018, [1512.03069].
  * (19) Belle-II collaboration, F. Abudinén et al., _Search for Axion-Like Particles produced in $e^{+}e^{-}$ collisions at Belle II_, _Phys. Rev. Lett._ 125 (2020) 161806, [2007.13071].
  * (20) BESIII collaboration, M. Ablikim et al., _Search for an axion-like particle in radiative J/ $\psi$ decays_, _Phys. Lett. B_ 838 (2023) 137698, [2211.12699].
  * (21) Belle-II collaboration, T. Abe et al., _Belle II Technical Design Report_ , 1011.0352.
  * (22) E. Fermi, _On the theory of collisions between atoms and electrically charged particles_ , _Nuovo Cim._ 2 (1925) 143–158, [hep-th/0205086].
  * (23) E. J. Williams, _Nature of the high-energy particles of penetrating radiation and status of ionization and radiation formulae_ , _Phys. Rev._ 45 (1934) 729–730.
  * (24) V. M. Budnev, I. F. Ginzburg, G. V. Meledin and V. G. Serbo, _The Two photon particle production mechanism. Physical problems. Applications. Equivalent photon approximation_ , _Phys. Rept._ 15 (1975) 181–281.
  * (25) S. Frixione, M. L. Mangano, P. Nason and G. Ridolfi, _Improving the Weizsacker-Williams approximation in electron - proton collisions_ , _Phys. Lett. B_ 319 (1993) 339–345, [hep-ph/9310350].
  * (26) Belle-II collaboration, W. Altmannshofer et al., _The Belle II Physics Book_ , _PTEP_ 2019 (2019) 123C01, [1808.10567].
  * (27) J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer 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].
  * (28) S. Stengel, C. Sfienti and T. Ferber, _Optimization of the pi0 reconstruction selections for the Belle II experiment_.  PhD thesis, Mainz, Johannes Gutenberg University, Mainz, 2019.
  * (29) C. Accettura et al., _Towards a muon collider_ , _Eur. Phys. J. C_ 83 (2023) 864, [2303.08533].
  * (30) M. Ruhdorfer, E. Salvioni and A. Wulzer, _Invisible Higgs boson decay from forward muons at a muon collider_ , _Phys. Rev. D_ 107 (2023) 095038, [2303.14202].
  * (31) S. Dawson, _The Effective W Approximation_ , _Nucl. Phys. B_ 249 (1985) 42–60.
  * (32) R. N. Cahn and S. Dawson, _Production of Very Massive Higgs Bosons_ , _Phys. Lett. B_ 136 (1984) 196.
  * (33) P. Borel, R. Franceschini, R. Rattazzi and A. Wulzer, _Probing the Scattering of Equivalent Electroweak Bosons_ , _JHEP_ 06 (2012) 122, [1202.1904].


# Surface Tension and Energy Conservation in a Moving Fluid

Tomas Bohr ID<EMAIL_ADDRESS>Department of Physics, Technical
University of Denmark, DK-2800 Kgs. Lyngby, Denmark  Bernhard Scheichl ID
<EMAIL_ADDRESS>Institute of Fluid Mechanics and Heat Transfer,
Technische Universität Wien, Getreidemarkt 9, 1060 Vienna, Austria  AC2T
research GmbH, 2700 Wiener Neustadt, Austria

###### Abstract

The transport of energy in a moving fluid with a simply connected free surface
is analyzed, taking into account the contribution of surface tension. This is
done by following a “control volume” with arbitrary, specified velocity,
independent of the flow velocity, and determining the rates of energy passing
through the boundaries, as well as the energy dissipation in the bulk. In
particular, a simple conservation equation for the surface area is written
down, which clearly shows the contribution of the Laplace pressure at the free
surface and the tangential surface tension forces at its boundary. It emerges
as the mechanical conservation law for the surface energy in its general form.
For a static control volume, all contributions from surface tension disappear,
except that the pressure has to be modified by the Laplace contribution.

Capillary effects have been included in flows with free surfaces or interfaces
for around 150 years. Some famous pioneers are Rayleigh [1, 2], Landau and
Levich [3], Taylor [4], Culick [5], and Bretherton [6] and it is still an
active field of research with outstanding survey papers [7, 8] and the
textbook [9]. In contrast to flows having free surfaces with objects
protruding, where surface tension gives rise to tangential forces at the
boundaries, or Marangoni flows, where gradients of the surface tension
contribute to the surface-stress balance, free-surface flows include a uniform
surface tension only via the discontinuity the normal component of the fluid
stresses on the interface, commonly referred to as Laplace pressure and
described, e.g., carefully by Batchelor [11] (§ 1.9).

Although this is merely a consequence of the governing equations of motion and
conventional kinematic and dynamic boundary conditions, this procedure has
been challenged recently by Bhagat et al. [12] and Bhagat and Linden [13].
These authors claim that the tangential surface tension contributes to a power
term in the integral energy budget and should, therefore, be taken into
account explicitly. For a static control volume in a steady (stationary) flow,
this power is claimed to be

$\oint_{\mathcal{C}_{s}}{\bf f}_{\gamma}\cdot{\bf
u}\,\mathrm{d}\mathcal{C}_{s}.$ (1)

Here the closed curve $\mathcal{C}_{s}$ has the infinitesimal arc length
$\mathrm{d}\mathcal{C}_{s}$ and encloses a portion of the free surface that
forms the control area and confines the control volume, ${\bf u}$ is the flow
velocity on $\mathcal{C}_{s}$ and ${\bf f}_{\gamma}$ is the external surface
tension force (pr. length). Indeed, inclusion of such a power term, referring
to the external work of surface tension, seems intuitively correct and its
omission therefore a serious error in all the preceding works. However, in
this paper we demonstrate that within the continuum hypothesis, where the
interface is taken as infinitely thin and mass- and inertialess, this is not
the case. If the control surface moves with velocity ${\bf U}$, such a term
appears whith ${\bf u}$ is now replaced by ${\bf U}$, so for a static control
volume there is no such term. Even in the general situation of an unsteady
flow, the appearance of this term is inconsistent with the advection of
surface energy combined with the standard work-energy balance for the bulk
flow, derived from the momentum (Navier–Stokes) equations. This was already
noticed by Duchesne et al. [14] and Scheichl [15], but in the present paper we
show how this can be seen directly from from the conservation law for surface
area that we derive below.

We consider a free-surface flow of an incompressible fluid (liquid) of uniform
density $\rho$, as sketched in Fig. 1. The free surface forming the interface
with the (gaseous) environment is characterised by a uniform surface tension
$\gamma$. The environment is typically taken as nominally at rest and under a
constant ambient pressure $p_{0}$. We assume that the dynamic viscosity $\mu$
for the gas is much smaller than that of the liquid, so the deviatoric
stresses from the external gas flow can be neglected (for simplicity). Let $t$
denote the time and ${\bf
x}=x\,{\bf\hat{x}}+y\,{\bf\hat{y}}+z\,{\bf\hat{z}}=(x,y,z)$ the space vector
expressed in static Cartesian laboratory coordinates $x$, $y$, $z$ and the
associated unit vectors (indicated with hats). For partial derivatives we use
$\nabla=(\partial_{x},\partial_{y},\partial_{z})$ as well as the subscripts
$t$, $x$, $y$, $z$,. We assume gravity to be the only body force at play, with
the constant gravitational acceleration ${\bf g}=(0,0,-g)$. In the bulk of the
flow, the Eularian “material” velocity field ${\bf u}({\bf x},t)=(u,v,w)$, the
pressure $p({\bf x},t)$, and the (symmetric) tensor $\bf T$ of the deviatoric
(viscous) Cauchy stresses, subject to a constitutive relationship, then
satisfy the continuity and momentum equations

$\displaystyle\nabla\cdot{\bf u}$ $\displaystyle=0,$ (2)
$\displaystyle\rho({\bf u}_{t}+{\bf u}\cdot\nabla{\bf u})$
$\displaystyle=-\nabla p-\rho g\,{\bf\hat{z}}+\nabla\cdot{\bf T}.$ (3)

For a Newtonian liquid, ${\bf T}=\mu\bigl{[}\nabla{\bf u}+(\nabla{\bf
u})^{T}\bigr{]}$ with the dynamic viscosity $\mu$.

${\bf
U}_{C}$${\bf\hat{m}}_{C}$${\bf\hat{n}}_{s}$$C$$y$$x$$\mathcal{V}$$z$${\bf\hat{n}}_{s}$$\mathcal{C}_{s}$$B$${\bf
U}$$u_{n}\>(U_{n})$${\bf\hat{m}}_{s}$${\bf\hat{l}}_{s}$${\bf f}_{\gamma}$$\bf
g$${\bf u}\>(\bf U)$${\bf u}_{s}\>({\bf U}_{s})$$z=h(x,y,t)$$\mathcal{A}_{s}$,

Figure 1: Control volume $\mathcal{V}$ with top face $\mathcal{A}_{s}$ lying
in free surface, projection $B$ of $\mathcal{A}_{s}$ onto $(x,y)$-plane
(bottom); body and surface forces, velocities, unit vectors; dotted:
projection and tangent lines.

We look at the simply-connected control area $\mathcal{A}_{s}$ of the moving
fluid surface with the boundary $\mathcal{C}_{s}=\partial\mathcal{A}_{s}$ and
consider also their normal projections $B$ and $C=\partial B$ onto on the
$(x,y)$-plane (the “bottom”). The free surface is described simplistically as
a function $z=h(x,y,t)$. The value of ${\bf u}$ on the surface can be
decomposed into a component in the surface (i.e., in the tangent plane of the
instantaneous surface) ${\bf u}_{s}(x,y,t)$ and a normal component
$u_{n}(x,y,t)$ along the outward normal direction, $n_{s}$, and surface unit
normal

${\bf\hat{n}}_{s}(x,y,t)\\!=\\!\frac{(-h_{x},-h_{y},1)}{\sqrt{\smash[b]{1+|\nabla
h|^{2}}}}\\!=\\!\frac{{\bf\hat{z}}-{\nabla h}}{\sqrt{\smash[b]{1+|\nabla
h|^{2}}}}.$ (4)

We furthermore introduce the vectors

${\bf\hat{l}}_{s}=\mathrm{d}{\bf
x}/\mathrm{d}\mathcal{C}_{s},\quad{\bf\hat{m}}_{s}={\bf\hat{l}}_{s}\times{\bf\hat{n}}_{s}\quad{\bf
f}_{\gamma}=\gamma\,{\bf\hat{m}}_{s}$ (5)

where $\mathrm{d}\mathcal{C}_{s}$ is the infinitesimal arc length on
$\mathcal{C}_{s}$, ${\bf\hat{l}}_{s}$ a unit tangent vectior on
$\mathcal{C}_{s}$, ${\bf\hat{m}}_{s}$ is a unit vector perpendicular to
$\mathcal{C}_{s}$, and ${\bf f}_{\gamma}$ is the tangential capillary force
pr. length from the surface region outside $\mathcal{C}_{s}$. Similarly, in
the $(x,y)$-plane, with infinitesimal arc length $\mathrm{d}C$ of $C$, its
outer normal ${\bf\hat{m}}_{C}$, and the projection of ${\bf\hat{n}}_{s}$ are

${\bf\hat{m}}_{C}=\frac{(\mathrm{d}y,-\mathrm{d}x,0)}{\mathrm{d}C},\quad{\bf
n}_{b}=\frac{-\nabla h}{\sqrt{\smash[b]{1+|\nabla h|^{2}}}}.$ (6)

Further $\bf u$ satisfies the kinematic free-surface boundary condition

$h_{t}=w-uh_{x}-vh_{y}=u_{n}\sqrt{1+|\nabla h|^{2}}$ (7)

on $\mathcal{A}_{s}$, so $u_{n}$ is only nonzero when $h$ is explicitly time-
dependent. The dynamic boundary condition expresses continuity of the total
stress across $\mathcal{A}_{s}$ in the form (see, e.g., [10], § 61)

$(p-p_{0})\,{\bf{\hat{n}}_{s}}-{\bf T}\cdot{\bf{\hat{n}}_{s}}=\Delta
p_{L}\,{\bf{\hat{n}}_{s}},$ (8)

where $\Delta p_{L}$ is the Laplace pressure involving the double local mean
curvature $\kappa$ of $\mathcal{A}_{s}$ as

$\Delta
p_{L}=\gamma\kappa=\gamma\nabla\cdot{\bf\hat{n}}_{s}=\gamma\nabla\cdot{\bf
n}_{b}$ (9)

(where the sign of $\kappa$ is chose so that $\kappa>0$ and $\nabla^{2}h<0$
for a convex, droplet-shaped surface). The equivalence of the local stress
equilibrium (8) on $\mathcal{A}_{s}$ with its global counterpart

$\int_{\mathcal{A}_{s}}\bigl{[}(p-p_{0})\,{\bf{\hat{n}}_{s}}-{\bf
T}\cdot{\bf{\hat{n}}_{s}}\bigr{]}\,\mathrm{d}A_{s}=-\\!\oint_{\mathcal{C}_{s}}{\bf
f}_{\gamma}\,\mathrm{d}\mathcal{C}_{s}$ (10)

is obtained using Stokes’ theorem in the form

$\\!\oint_{\mathcal{C}_{s}}{\bf\hat{m}}_{s}\,\mathrm{d}\mathcal{C}_{s}=\\!\oint_{\mathcal{C}_{s}}{\bf\hat{l}}_{s}\times{\bf\hat{n}}_{s}\,\mathrm{d}\mathcal{C}_{s}=-\int_{\mathcal{A}_{s}}(\nabla\cdot{\bf\hat{n}}_{s}){\bf\hat{n}}_{s}\,\mathrm{d}A_{s}$
(11)

as shown, e.g., in [16], with the counterintuitive result that an integral
over tangential vectors points in the direction of the normals.

It is our goal to combine the (hitherto unappreciated) conservation law for
the surface energy contained in $\mathcal{A}_{s}$ with that for the mechanical
energy of the bulk flow in an arbitrarily moved, simply-connected control
volume $\mathcal{V}$ bounded by $\mathcal{A}_{s}$ at the free surface. We
denote the velocity field moving the closed boundary $\partial\mathcal{V}$ as
${\bf U}({\bf x},t)=(U,V,W)$ and denote its components tangential and normal
to $\mathcal{A}_{s}$ as $U_{s}$ and $U_{n}$, respectively. Analogously, $C$
moves with ${\bf U}_{C}$, the projection of $\bf U$ evaluated on
$\mathcal{C}_{s}$ on the $(x,y)$-plane. Importantly, $\bf U$ must be
kinematically compatible to the material flow field $\bf u$ in the sense that
any point on the free surface considered to move with $\bf U$ remains in the
surface for all times. That is, $\bf U$ satisfies

$h_{t}=W-Uh_{x}-Vh_{y}=U_{n}\sqrt{1+|\nabla h|^{2}}$ (12)

on $\mathcal{A}_{s}$, in analogy to (7). Thus a completely static control
volume $\mathcal{V}$ is of course only possible for a steady flow.

The surface energy stored in $\mathcal{A}_{s}$ is simply $E_{s}=\gamma A_{s}$
where

$A_{s}=\int_{\mathcal{A}_{s}}\,\mathrm{d}A_{s}=\int_{B}\sqrt{1+\left|\nabla
h\right|^{2}}\,\mathrm{d}x\mathrm{d}y$ (13)

is the time-dependent area of $\mathcal{A}_{s}$. The integral conservation law
for the surface energy is found by analysing
$\gamma\,\mathrm{d}A_{s}/\mathrm{d}t$ in terms of $\bf U$, which transports
$\mathcal{V}$ and thus $\mathcal{A}_{s}$. Using Leibnitz’ integral rule, the
time derivative of (13) is

$\displaystyle\frac{\mathrm{d}A_{s}}{\mathrm{d}t}=$
$\displaystyle\int_{B}\partial_{t}\sqrt{1+|\nabla
h|^{2}}\,\mathrm{d}x\mathrm{d}y$ $\displaystyle+$
$\displaystyle\oint_{C}\sqrt{1+|\nabla h|^{2}}\,{\bf
U}_{C}\cdot{\bf\hat{m}}_{C}\,\mathrm{d}C.$ (14)

The first term to the right can be rewritten as

$\displaystyle\int_{B}\frac{\nabla h_{t}\cdot\nabla
h}{\sqrt{\smash[b]{1+|\nabla h|^{2}}}}\,\mathrm{d}x\mathrm{d}y$
$\displaystyle=\int_{B}\left[-h_{t}\nabla\cdot\left(\frac{\nabla
h}{\sqrt{\smash[b]{1+|\nabla
h|^{2}}}}\right)+\nabla\cdot\left(\frac{h_{t}\nabla
h}{\sqrt{\smash[b]{1+|\nabla h|2}}}\right)\right]\mathrm{d}x\mathrm{d}y$
$\displaystyle=\int_{B}h_{t}\,\kappa\,\mathrm{d}x\mathrm{d}y+\oint_{C}\frac{h_{t}}{\sqrt{\smash[b]{1+|\nabla
h|^{2}}}}\nabla h\cdot{\bf\hat{m}}_{C}\,\mathrm{d}C.$ (15)

The last equality follows from (9) and the divergence theorem. All in all, we
have

$\displaystyle\frac{\mathrm{d}A_{s}}{\mathrm{d}t}$
$\displaystyle=\int_{B}h_{t}\,\kappa\,\mathrm{d}x\mathrm{d}y$
$\displaystyle+\oint_{C}\left(\sqrt{1+|\nabla h|^{2}}\,{\bf
U}+\frac{h_{t}}{\sqrt{\smash[b]{1+|\nabla h|^{2}}}}\nabla
h\right)\cdot{\bf\hat{m}}_{C}\,\mathrm{d}C.$ (16)

The expressions (5) combined with (4) yield the explicit representation

${\bf\hat{m}}_{s}=\frac{\left[(1+h_{y}^{2})\mathrm{d}y,-(1+h_{x}^{2})\mathrm{d}x,h_{x}\mathrm{d}y-h_{y}\mathrm{d}x)\right]}{\sqrt{\smash[b]{1+|\nabla
h|^{2}}}\,\mathrm{d}\mathcal{C}_{s}}.$ (17)

and together with the kinematic boundary condition (12) this leads to the
simple result

$\frac{\mathrm{d}A_{s}}{\mathrm{d}t}=\int_{\mathcal{A}_{s}}U_{n}\,\kappa\,\mathrm{d}A_{s}+\oint_{\mathcal{C}_{s}}{\bf
U}\cdot{\bf\hat{m}}_{s}\,\mathrm{d}\mathcal{C}_{s}.$ (18)

Actually this result could have been obtained more directly and intuitively by
considering independently the three infinitesimally small temporal (linear)
variations of the shape of the material surface element $\mathrm{d}A$: 1.
dilatation normal to $\mathcal{A}_{s}$ (inflation), 2. distortion in its
tangential directions (stretching and straining) and 3. tangential solid-body
rotation. Since the latter has no effect, this leaves us with the two
contributions on the right side of (18), where the first term originates in
the relation $(\delta\mathrm{d}A_{s}/\delta n_{s})/\mathrm{d}A_{s}$ (see,
e.g., [18] and [10], § 61) and the last term describes the tangential change
of $A$ as the change of area adjacent to $\mathcal{C}_{s}$.

It is instructive to write down the analog of (18) for a planar flow where the
free surface is a line $z=h(x,t)$ as shown Fig. 2. The control volume (area)
has endpoints $x=x_{1}(t)$ and $x=x_{2}(t)>x_{1}$ on the free surface move
with the imposed velocity $U(x,t)$. The analog of the surface are, i.e. the
length of the free surface segment is

$L_{s}=\int_{\mathcal{L}_{s}}\mathrm{d}L_{s}=\int_{x_{1}}^{x_{2}}\sqrt{1+h_{x}^{2}}\,\mathrm{d}x$
(19)

and the curvature $\kappa_{1}$ is

$\kappa_{1}=-\frac{\partial}{\partial
x}\left(\frac{h_{x}}{\sqrt{\smash[b]{1+h_{x}^{2}}}}\right)=-\frac{h_{xx}}{(1+h_{x}^{2})^{3/2}}.$
(20)

If $d$ denotes its depth in the $y$-direction, $E_{s}=\gamma L_{s}d$ is the
corresponding surface energy. In order to consider the advection of a point on
this segment by a velocity field ${\bf U}(x,y,t)=(U,0,W)$, we take its
$x$-coordinate as time-dependent, thus its $z$-coordinate as $h(x(t),t)$. On
$\mathcal{C}_{s}$, $U=\dot{x}$ and $W=h_{t}+\dot{x}h_{x}$ by (12). It is
evident that ${\bf U}\cdot{\bf\hat{m}}_{s}$ vanishes on the points of
$\mathcal{C}_{s}$ with $x$-coordinates located between $x_{1}$ and $x_{2}$.
Here
${\bf{\hat{m}}}_{s}={(\bf\hat{x}}+h_{x}\,{\bf\hat{z}})/\sqrt{\smash[b]{1+h_{x}^{2}}}$
and $\bf U$ are collinear if $h_{t}$ vanishes. These results help us to
rewrite (18) in this case as

$\displaystyle\frac{\mathrm{d}L_{s}}{\mathrm{d}t}$
$\displaystyle=\int_{x_{1}}^{x_{2}}h_{t}\,\kappa_{1}\,\mathrm{d}x+\left[\sqrt{1+h_{x}^{2}}\,\dot{x}+\frac{h_{t}\,h_{x}}{\sqrt{\smash[b]{1+h_{x}^{2}}}}\right]_{x=x_{1}}^{x=x_{2}}$
$\displaystyle=\int_{\mathcal{L}_{s}}U_{n}\,\kappa_{1}\,\mathrm{d}L_{s}+\left[{\bf
U}\cdot{\bf\hat{m}}_{s}\right]_{x=x_{1}}^{x=x_{2}}$ (21)

in direct analogy with (16) by noticing the projection ${\bf
U}_{C}=\dot{x}_{1,2}\,{\bf\hat{x}}$ of $\bf U$ onto the endpoints.
Alternatively, one obtains (21) from (18) in straightforward manner under the
assumptions of a planar flow and $\mathcal{C}$ choosen such that $C$ is a
rectangle with a side of length $x_{2}-x_{1}$ parallel to the $x$-axis.

$z=h(x,t)$${\bf\hat{m}}_{s}$$\mathcal{A}_{s}$$B$$h_{t}$$U_{n}$${\bf
U}$$z$$x_{1}$$x_{2}$$\dot{x}$${\bf\hat{m}}_{s}$${\bf f}_{\gamma}$${\bf
f}_{\gamma}$${\bf U}_{s}$$W$${\bf U}_{C}$$x$${\bf U}_{C}$

Figure 2: Planar flow: projection onto $(x,z)$-plane; for legend see Fig. 1.

We next remark some important consequences of (18). First we derive the
differential version using the divergence theorem. With the normal projection
$\nabla_{s}=\nabla-{\bf\hat{n}}_{s}\partial_{n}$ of $\nabla$ onto
$\mathcal{A}_{s}$ acting on the tangential flow field ${\bf U}_{s}={\bf
U}-U_{n}{\bf\hat{n}}_{s}$,

$\oint_{\mathcal{C}_{s}}{\bf
U}\cdot{\bf\hat{m}}_{s}\,\mathrm{d}\mathcal{C}_{s}=\int_{\mathcal{A}_{s}}\nabla_{s}\cdot{\bf
U}_{s}\,\mathrm{d}A_{s}.$ (22)

Using $\partial_{n}\equiv{\bf\hat{n}}_{s}\cdot\nabla$, we expand the
divergence in the surface integral in (22) as $\nabla_{s}\cdot{\bf
U}_{s}=\nabla\cdot{\bf
U}-{\bf\hat{n}}_{s}\cdot~{}({\bf\hat{n}}_{s}\cdot\nabla{\bf
U})-\nabla_{s}\cdot(U_{n}{\bf\hat{n}}_{s})$. Since $\nabla_{s}U_{n}$ lies in
the surface and
$\nabla_{s}\cdot{\bf\hat{n}}_{s}\equiv\nabla\cdot{\bf\hat{n}}_{s}=\kappa$, cf.
(9), (18) takes on the form

$\frac{\mathrm{d}A_{s}}{\mathrm{d}t}=\int_{\mathcal{A}_{s}}[U_{n}\nabla\cdot{\bf\hat{n}}_{s}+\nabla_{s}\cdot{\bf
U}_{s}]\,\mathrm{d}A_{s}=\int_{\mathcal{A}_{s}}\bigl{[}\nabla\cdot{\bf
U}-{\bf\hat{n}}_{s}\cdot\left(\nabla{\bf
U}\right)\cdot{\bf\hat{n}}_{s}\bigr{]}\,\mathrm{d}A_{s}.$ (23)

As $\mathcal{A}_{s}$ is chosen arbitrarily, the relative temporal rate of
change of an infinitesimal area element must be

$(\delta\mathrm{d}A_{s}/\delta
t)/\mathrm{d}A_{s}=U_{n}\nabla\cdot{\bf\hat{n}}_{s}+\nabla_{s}\cdot{\bf
U}_{s}=\nabla\cdot{\bf U}-{\bf\hat{n}}_{s}\cdot\left(\nabla{\bf
U}\right)\cdot{\bf\hat{n}}_{s}$ (24)

For $\bf U$ identical to the material flow field $\bf u$, this local form of
(23) was already derived by Batchelor [11], eqs. (3.1.5)–(3.1.8), without
taking notice of its integral formulation. For an incompressible material flow
field, it takes the form $(\delta\mathrm{d}A_{s}/\delta
t)/\mathrm{d}A_{s}=-{\bf\hat{n}}_{s}\cdot\left(\nabla{\bf
u}\right)\cdot{\bf\hat{n}}_{s}$ and these local forms makes explicit the
required invariance under the Galilean boost, i.e., for (24): ${\bf
U}\mapsto{\bf U}+{\bf C}$, ${\bf x}\mapsto{\bf x}+{\bf C}t$, which is evident
neither from (18) nor from (16).

Finally, for a stationary velocity field $\bf U$, the curvature term in (18)
is zero. Interestingly and counterintuitively, then
$\mathrm{d}A_{s}/\mathrm{d}t$ does not necessarily vanish even if ${\bf
U\equiv\bf u}$ represents an incompressible flow satisfying (2). This
suppresses the dilatation of the fluid particles but not their stretching and
straining in the surface governed by the surface divergence in (22).
Therefore, the streamlines lying in the surface then can diverge/converge,
which explains the in general non-vanishing $\mathrm{d}A_{s}/\mathrm{d}t$.

In the remainder of this letter, we discuss the dynamic implications of the
kinematic relationship (18). Multiplying (18) with $\gamma$ and using the
expressions for the Laplace pressure in (8) and the vectorial surface tension,
${\bf f}_{\gamma}=\gamma\,{\bf\hat{m}}_{C}$, results in the integral transport
equation for the surface energy convected by $\bf U$:

$\frac{\mathrm{d}E_{s}}{\mathrm{d}t}=\int_{\mathcal{A}_{s}}\Delta
p_{L}U_{n}\,\mathrm{d}A_{s}+\oint_{\mathcal{C}_{s}}{\bf f}_{\gamma}\cdot{\bf
U}\,\mathrm{d}\mathcal{C}_{s}.$ (25)

This has an appealing and convincing physical interpretation as a law of
mechanical energy conversion for the piece $\mathcal{A}_{s}$ of the
inertialess free fluid surface: the work done by the Laplace pressure on
$\mathcal{A}_{s}$ and by the tangential surface tension forces at its boundary
$\mathcal{C}_{s}$, on the surface moving with the imposed velocity field ${\bf
U}$, changes the potential internal energy, i.e. the surface energy $E_{s}$,
contained in $\mathcal{A}_{s}$. Therefore, (25) represents the precise
justification for the Laplace pressure and the existence of the tangential
surface tension forces. It is readily seen that (25) can also be derived from
the pointwise relationships (9), after multiplication with $U_{n}$, and (24).
Note that the inclusion of Marangoni effect (inhomogeneous $\gamma$) would
lead to an extension of (25) incorporate the transport of $\gamma$ as a non-
uniform passive scalar [19].

The conservation law (25) is the basis for the subsequent careful re-
examination of the work-energy balance for the entire volume $\mathcal{V}$
including its top face $\mathcal{A}_{s}$ once ${\bf U}$ is specified
appropriately. To do this, we first briefly recall the derivation of the
classical work-energy theorem as stated, e.g., in eq. (16.2) in [10]. More
precisely, it denotes the transport of the overall mechanical energy of the
bulk flow,

$E_{m}=\int_{\mathcal{V}}e_{m}\,\mathrm{d}V,\quad e_{m}=\frac{\rho|{\bf
u}|^{2}}{2}+\rho gz,$ (26)

where the energy density $e_{m}$ consists of the kinetic and the potential
energy densities. Taking the inner product of ${\bf u}$ and the momentum
equation (3) together with the incompressibility condition (2) and the
identity $w\equiv\nabla\cdot(z{\bf u})$ leads to the differential form of the
energy equation:

$\partial_{t}e_{m}=\nabla\cdot[{\bf T}\cdot{\bf u}-(p+e_{m}){\bf u}]-{\bf
T}\cdot\nabla{\bf u}.$ (27)

Integration of (27) over $\mathcal{V}$ and applying the divergence theorem
gives

$\int_{\mathcal{V}}\partial_{t}e_{m}\,\mathrm{d}V=\oint_{\partial\mathcal{V}}[{\bf
T}\cdot{\bf u}-(p+e_{m}){\bf
u}]\cdot{\bf\hat{n}}\,\mathrm{d}A-\\!\int_{\mathcal{V}}{\bf T}\cdot\nabla{\bf
u}\,\mathrm{d}V$ (28)

where the first term represents the advected energy, and the last term
represents the (non-negative) integral dissipation per time unit due to the
deviatoric stresses, commonly denoted with $\dot{\Phi}$. We have also
introduced the infinitesimal area element $\mathrm{d}A$ of
$\partial\mathcal{V}$ and its outward unit normal $\bf\hat{n}$. We note that,
according to the derivation of (28), the stresses acting on
$\partial\mathcal{V}$ are taken as their limiting values approached from
inside $\mathcal{V}$. For the bulk flow this makes no difference due to
continuity, but when we include the free surface it becomes important. The
transport equation then follows from the classical Leibnitz–Reynolds transport
theorem (cf. [20], eqs. (3.1.3)–(3.1.5) in [11]):

$\displaystyle\frac{\mathrm{d}E_{m}}{\mathrm{d}t}$
$\displaystyle=\int_{\mathcal{V}}\partial_{t}e_{m}\,\mathrm{d}V+\oint_{\partial\mathcal{V}}e_{m}{\bf
U}\cdot{\bf\hat{n}}\,\mathrm{d}A$
$\displaystyle=\oint_{\partial\mathcal{V}}\bigl{[}e_{m}({\bf U}-{\bf u})+{\bf
T}\cdot{\bf u}-p\,{\bf u}\bigr{]}\cdot{\bf\hat{n}}\,\mathrm{d}A-\dot{\Phi}$
(29)

This work-energy principle states that the temporal change of $E_{m}$ is
compensated by the flux of mechanical energy into $\mathcal{V}$ and the
external work done by the total stresses acting on $\partial\mathcal{V}$, due
to the material velocity ${\bf u}$, and the associated internal work of the
deviatoric (viscous) stresses, equal to the negative mechanical dissipation.
It equals zero for a stationary flow.

For many application, especially for stationary flows, it is valuable to
consider a cylindrical control volume $\mathcal{V}$ aligned with $z$-direction
and having the patch $\mathcal{A}_{s}$ of the free surface as its top
boundary, as indicated in Fig. 1. In the general, non-stationary situation its
lateral sides facing the bulk of the flow and $\mathcal{A}_{s}$ are both
moving. In this important case of a static control volume, its projection $B$
and thus $C$ stay stationary while $\mathcal{A}_{s}$ is moving though with the
normal component of the material flow (when the flow is not stationary). Thus
${\bf U}={\bf 0}$ on the lateral sides of $\mathcal{V}$ and ${\bf
U}_{C}\equiv{\bf 0}$ but ${\bf U}=u_{n}{\bf\hat{n}}_{s}$,
${\bf\hat{n}}={\bf\hat{n}}_{s}$, $\mathrm{d}A=\mathrm{d}A_{s}$ on its top face
$\mathcal{A}_{s}$. In this case, (29) is written as

$\displaystyle\frac{\mathrm{d}E_{m}}{\mathrm{d}t}=\frac{\mathrm{d}}{\mathrm{d}t}\int_{0}^{h(x,y,t)}\mathrm{d}z\int_{B}e_{m}\,\mathrm{d}x\mathrm{d}y$
$\displaystyle=-\\!\int_{\mathcal{A}_{s}}\Delta
p_{L}u_{n}\,\mathrm{d}A_{s}+\\!\int_{\partial\mathcal{V}_{r}}\left[{\bf
T}\cdot{\bf u}-(e_{m}+p-p_{0}){\bf
u}\right]\cdot{\bf\hat{n}}\,\mathrm{d}A-\dot{\Phi}$ (30)

where $\partial\mathcal{V}_{r}$ denotes the remaining boundary, i.e., the
boundary excluding the free surface, and where we have subtracted the constant
$p_{0}$ from the pressure everywhere since $\oint_{\partial\mathcal{V}}{\bf
u}\cdot{\bf\hat{n}}\,\mathrm{d}A=0$ by continuity and incompressibility.To
determine the total energy budget, we note that the Laplace pressure acts on
both sides of the interface separating $\mathcal{A}_{s}$ and $\mathcal{V}$.
But the powers it exerts on $\mathcal{A}_{s}$, see (25), and on $\mathcal{V}$,
see (30), differ in theirs signs and therefore cancel. Since we saw from the
comparison with (14) that the loop integral in (18) vanishes for ${\bf
U}_{C}\equiv{\bf 0}$, this yields the total energy budget in the form

$\frac{\mathrm{d}(E_{m}+E_{s})}{\mathrm{d}t}=\\!\int_{\partial\mathcal{V}_{r}}\left({\bf
T}\cdot{\bf u}-(e_{m}+p-p_{0}){\bf
u}\right)\cdot{\bf\hat{n}}\,\mathrm{d}A-\dot{\Phi}.$ (31)

For stationary flow $\mathrm{d}E_{s}/\mathrm{d}t$ vanishes as does $u_{n}$;
see (18). Then the energy equation balance (31), combining that of the bulk
and the surface flow, becomes indistinguishable from that just governing the
bulk flow, (30). Correspondingly, the only trace of capillarity then comes
from the pressure $p$, which must include the Laplace pressure to respect the
boundary condition (8)-(9). This is the standard approach, and it validity was
recently emphasized in [14] and [15]. If this seems strange, one should keep
in mind that surface area is not advected like mass or kinetic energy. In the
bulk of the flow, each fluid particle carries a particular mass and kinetic
energy with it, but it does not carry a particular surface area with it.

The additional energy term (1) for a static control volume in a steady flow
was originally postulated by [12] without explanation. In their follow-up
paper [13] they argue that it actually stems from the viscous advection terms
in (29). Writing the global stress balance (10) in the form (their eq. (3.4)
in our notation)

$\int_{\mathcal{A}_{s}}(p-p_{0}){\bf\hat{n}}_{s}\mathrm{d}A_{s}-\int_{\mathcal{A}_{s}}{\bf\hat{n}}_{s}\cdot{\bf
T}\,\mathrm{d}A_{s}+\gamma\\!\oint_{\mathcal{C}_{s}}{\bf\mathrm{d}l}\times{\bf\hat{n}}_{s}=0,$
(32)

they then take the inner product with $\bf u$ inside both integrals to get

$\int_{\mathcal{A}_{s}}(p-p_{0}){\bf
u}\cdot{\bf\hat{n}}_{s}\mathrm{d}A_{s}-\int_{\mathcal{A}_{s}}{\bf u}\cdot{\bf
T}\cdot{\bf\hat{n}}_{s}\,\mathrm{d}A_{s}+\gamma\\!\oint_{\mathcal{C}_{s}}{\bf
u}\cdot({\bf\mathrm{d}l}\times{\bf\hat{n}}_{s})=0$ (33)

and conclude, since ${\bf u}\cdot{\bf\hat{n}}_{s}$ (for a stationary flow),
that

$\int_{\mathcal{A}_{s}}{\bf u}\cdot{\bf
T}\cdot{\bf\hat{n}}_{s}\,\mathrm{d}A_{s}=-\gamma\\!\oint_{\mathcal{C}_{s}}{\bf
u}\cdot({\bf\mathrm{d}l}\times{\bf\hat{n}}_{s}),$ (34)

relating viscous and capillary effects and reinstating the viscous term that
was dropped from the free surface part of the energy equation (31). However,
multiplying by ${\bf u}$ inside the integral, which might seem a legal
operation since the area of integration is arbitrary, is actually illegal
since one of the integrals is a loop integral. Rewriting the loop integral in
(32) as a surface integral, using the identity (11), we get (as in (10)–(11))

$\int_{\mathcal{A}_{s}}\big{[}(p-p_{0}){\bf\hat{n}}_{s}-{\bf
T}\cdot{\bf\hat{n}}_{s}-\gamma(\nabla\cdot{\bf\hat{n}}_{s}){\bf\hat{n}}_{s}\big{]}\mathrm{d}A_{s}=0$
(35)

and we can now take ${\mathcal{A}_{s}}$ as any shape we want and conclude –
since the first and last terms are explicitly orthogonal to ${\bf u}$ – that

$\int_{\mathcal{A}_{s}}{\bf u}\cdot{\bf
T}\cdot{\bf\hat{n}}_{s}\,\mathrm{d}A_{s}=0$ (36)

Indeed the integrand ${\bf u}_{s}\cdot{\bf T}\cdot{\bf\hat{n}}_{s}$ (which we
removed from the free surface to get (31)) must be zero everywhere on the
interface, which is simply the shear part of the usual dynamic boundary
condition (8) on a free surface, and does not involve the surface tension.

As we mentioned above, it is perhaps strange that the surface tension forces
${\bf f}_{\gamma}$ in (25) do no work when the control volume does not move,
since these forces are still acting on a moving fluid. The reason for this is
that the surface energy, being proportional to the surface area (13), does not
depend on the fluid velocity ${\bf u}$. In other words, the surface energy has
no kinetic energy term, which is due to the fact that the surface, in our
idealized continuum approach, is infinitely thin and thus has no mass. But how
can it then have an energy? The reason for this is found in the fact that the
surface energy, coming from the breaking of bonds between the liquid
molecules, is at least of the order of $kT$ pr. particle or $RT$ pr. mol, i.e.
around 2 kJ/mol ($k$, $R$, and $T$ denote the Boltzmann and the universal gas
constant, and the absolute temperature, respectively), whereas the kinetic
energy of matter, say water, moving at 1 m/s is around 0.02 J/mol. The ratio
of these two energies of around 105 allows us to neglect the kinetic energy of
the organized motion represented by the velocity field ${\bf u}$ compared to
the surface energy for the same amount of matter – as opposed to the kinetic
energy of the random thermal motion, which is also of the order of $kT$ pr.
particle.

This work was co-funded by the project COMET InTribology, FFG-No. 872176
(project coordinator: AC2T research GmbH, Austria). Tomas Bohr is grateful to
Anders Andersen and Alexis Duchesne for many valuable discussions and to Jens
Eggers for useful comments.

The authors contributed equally to this work.

## References

  * [1] J.W.S. Rayleigh, F.R.S., On the instability of jets, P. Lond. Math. Soc. s1-10, (1), 4 (1878).
  * [2] J.W.S. Rayleigh, F.R.S., On the stability, or instability, of certain fluid motions, P. R. Soc. London 29 (196)–(199), 79 (1879).
  * [3] L.D. Landau and B.V. Levich, Dragging of a liquid by a moving plate, Acta Physicochim. USSR 17 (1–2), 42 (1942).
  * [4] G.I. Taylor, The dynamics of thin sheets of fluid. III. Disintegration of fluid sheets, Proc. R. Soc. A–Math. Phys. 253 (1274), 313 (1959).
  * [5] F.E.C. Culick, Comments on a Ruptured Soap Film, J. Appl. Phys. 31 (6), 1128 (1960).
  * [6] F.P. Bretherton, The motion of long bubbles in tubes, J. Fluid. Mech. 10 (2), 166 (1961).
  * [7] J. Eggers, Nonlinear dynamics and breakup of free-surface flows, Rev. Mod. Phys. 69 (3), 865 (1997).
  * [8] D. Bonn, J. Eggers, J. Indekeu, J. Meunier, and E. Rolley, Wetting and spreading, Rev. Mod. Phys. 81 (2), 739 (2009).
  * [9] P.-G. de Gennes, F. Brochard-Wyart, and D. Quéré, Capillarity and Wetting Phenomena (Drops, Bubbles, Pearls, Waves), Springer, (2004).
  * [10] L.D. Landau and E.M. Lifshitz, Fluid Dynamics, 2nd ed., Course of Theoretical Physics, vol. 6, Pergamon (1987).
  * [11] G.K. Batchelor, An Introduction to Fluid Dynamics, Cambridge Mathematical Library (2000).
  * [12] R.K. Bhagat, N.K. Jha, P.F. Linden, and D.I. Wilson, On the origin of the circular hydraulic jump in a thin liquid film, J. Fluid Mech. 851, R5 (2018).
  * [13] R.K. Bhagat and P.F. Linden, The circular capillary jump, J. Fluid Mech. 851, A25 (2020).
  * [14] A. Duchesne, A. Andersen, and T. Bohr, Surface tension and the origin of the circular hydraulic jump in a thin liquid film, Phys. Rev. Fluids 4 (8), 084001 (2019).
  * [15] B. Scheichl, Centred Splash of a Vertical Jet on a Horizontal Rotating Disc: Recent Findings and Resolving Controversies Over the Hydraulic Jump, Boundary and Interior Layers: Computational and Asymptotic Methods (BAIL 2018), Lecture Notes in Computational Science and Engineering, Springer book series, vol. 135, chap. 13, p. 201 (2020).
  * [16] C. Pozrikidis, Introduction to Theoretical and Computational Fluid Dynamics, 2nd ed., Oxford, 2011\.
  * [17] H. Flanders, Differentiation under the integral Sign, Am. Math. Mon. 80 (6), 615 (1973).
  * [18] B. Chow, P. Lu, and L. Ni, Hamilton’s Ricci Flow, Graduate Studies in Mathematics, vol. 77, AMS (2006).
  * [19] H.A. Stone, A simple derivation of the time‐dependent convective‐diffusion equation for surfactant transport along a deforming interface, Phys. Fluids A–Fluid 2 (1), 111 (1990).
  * [20] L.G. Leal, Advanced Transport Phenomena: Fluid Mechanics and Convective Transport Processes, Cambridge Series in Mechanical Engineering, vol. 7, 2007.
  * [21] B. Mohajer and R. Li, Circular hydraulic jump on finite surfaces with capillary limit, Phys. Fluids 27 (11), 117102 (2015).
  * [22] R.I. Bowles and F.T. Smith, The standing hydraulic jump: theory, computations and comparisons with experiments, J. Fluid Mech. 242, 145 (1992).