(LABEL:)

# Robust estimation of exposure ratios in multi-exposure image stacks

Param Hanji1, Rafał K. Mantiuk2
Department of Computer Science and Technology, University of Cambridge
Email<EMAIL_ADDRESS>
###### Abstract

Merging multi-exposure image stacks into a high dynamic range (HDR) image
requires knowledge of accurate exposure times. When exposure times are
inaccurate, for example, when they are extracted from a camera’s EXIF
metadata, the reconstructed HDR images reveal banding artifacts at smooth
gradients. To remedy this, we propose to estimate exposure ratios directly
from the input images. We derive the exposure time estimation as an
optimization problem, in which pixels are selected from pairs of exposures to
minimize estimation error caused by camera noise. When pixel values are
represented in the logarithmic domain, the problem can be solved efficiently
using a linear solver. We demonstrate that the estimation can be easily made
robust to pixel misalignment caused by camera or object motion by collecting
pixels from multiple spatial tiles. The proposed automatic exposure estimation
and alignment eliminates banding artifacts in popular datasets and is
essential for applications that require physically accurate reconstructions,
such as measuring the modulation transfer function of a display. The code for
the method is available.

###### Index Terms:

High dynamic range imaging; camera noise model, statistical estimation, multi-
exposure fusion

††footnotetext: This work is published in the Transactions of Computational
Imaging, 2023. © 2023 IEEE. Personal use of this material is permitted.
Permission from IEEE must be obtained for all other uses, in any current or
future media, including reprinting/republishing this material for advertising
or promotional purposes, creating new collective works, for resale or
redistribution to servers or lists, or reuse of any copyrighted component of
this work in other works. Digital Object Identifier 10.1109/TCI.2023.3301338

##  I Introduction

When the dynamic range of a scene exceeds the operating range of a standard
digital sensor, one can overcome this limitation by capturing a stack of
images with different camera settings: modulating exposure times [13], sensor
gain [16, 22] or by capturing and averaging a burst of images [20]. Then, the
captured exposure stack can be merged into a single image to both expand the
dynamic range and reduce noise.

Regardless of the approach, regions of the reconstructed high dynamic range
(HDR) images that contain smooth intensity gradients often end up with banding
artifacts, such as those shown in Figure 1. Banding artifacts are highly
visible in both images and videos, and they have been the subject of extensive
research, particularly in the field of video streaming. Banding tends to be
more noticeable compared to other video compression artifacts [44, 40, 27].
The reason for banding in merged exposure stacks is a mismatch between the
actual capture parameters and those reported by the camera, typically in the
exchangeable image file format (EXIF) header. Such inaccuracies could be
caused by

  * •

limited accuracy of the (mechanical) aperture and shutter

  * •

wrongly reported EXIF data due to rounding (e.g. $\nicefrac{{1}}{{60}}$
exposure time could be reported instead of $\nicefrac{{1}}{{64}}$)

  * •

changes in scene illumination, for example, due to flickering lights or
overcast skies with moving clouds

We stress that the artifacts shown in Figure 1 are not due to spatial
misalignment because we carefully selected image stacks without camera or
object motion. Further, since we merged linear RAW images, the artifacts are
not due to incorrect camera response function (CRF) inversion.

The inaccuracies of exposure time found in the EXIF metadata are substantial,
reaching 40%, as shown in Figure 2. To prepare this figure, we computed the
relative error between the reported EXIF exposure times and the times
estimated by the proposed method. The two plots show similar error
distribution for two datasets of multi-exposure HDR image stacks [10, 18]. The
errors of such magnitude can easily introduce banding artifacts, as shown in
Figure 1.

Figure 1: The first row depicts patches from images of overcast skies taken
from large HDR datasets [10, 18]. Misalignment between the reported camera
settings results in banding artifacts (pointed at by the red arrows) at the
boundary between two exposures. This problem is not limited to the sky and
affects other regions too, as shown by patches in the second row. We carefully
picked exposure stacks consisting of RAW spatially-aligned images that do not
require CRF inversion or motion compensation. Figure 2: Histograms of relative
errors for two popular HDR datasets. The computed sample standard deviations
of relative errors are $9\%$ for the Fairchild survey [10] and $14\%$ for the
SI-HDR dataset [18].

In this work, we propose to estimate exposure ratios directly from a stack of
images, which lets us eliminate the banding artifacts from the merged HDR
images. We formulate the problem to be efficiently solved as a sparse linear
system while accounting for heteroskedastic camera noise with known or unknown
noise parameters. Thus, our exposure estimation can be used to improve ground
truth reconstructions of several existing stack-based HDR datasets. Further,
algorithms that enhance merged reconstructions by iterative optimization [29],
as well as deep neural networks that use HDR images for other related tasks
[46, 6] stand to be improved by our more accurate exposures.

Correctly estimating exposures is also relevant for HDR merging with
deghosting, which involves detecting and registering pixels that belong to
moving objects [41, 32]. Erroneous capture parameters make this task more
challenging, as pixel differences could be caused by either object movement or
incorrect exposure time.

Although banding artifacts are visible only in smooth image regions in Figure
1, inaccurate exposures affect all pixels in the image. Incorrect pixel
intensities pose a serious complication for applications that utilize cameras
in place of expensive light measurement instruments. Multi-image exposure
stacks often serve as substitutes for more accurate measurements from
instruments such as spectrophotometers. In Section IV-D, we highlight
potential discrepancies in the measured modulation transfer function (MTF) of
an HDR display when the exposure stack of a slant-edge [36] is merged with
incorrectly reported capture parameters.

From a practical standpoint, our method is suitable for large images (8k or
more) from modern cameras. We achieve this by solving a smaller system of
equations by collecting pixel-pairs with the lowest relative variances
according to the noise model we consider. Contrary to existing methods that
utilize pixel pairs from consecutive exposures in the stack, we introduce a
greedy algorithm that provides optimal pixel pairs to ensure that all
exposures are well-estimated. Our reduced linear system is a union of the
highest weighted spanning trees induced on an _exposure multigraph_ ,
resulting in a balanced linear system. Additionally, our proposed method is
spatially balanced because we split the input image into tiles and collect
pixels from all tiles. This improves the robustness of our estimator to
ghosting caused by camera or object motion.

Here is a brief overview of the paper that summarizes our contributions. We
highlight that camera metadata may be unreliable and motivate the need to
estimate corrective per-exposure ratios to obtain artifact-free HDR
reconstruction. We show that such exposure ratios can be estimated by solving
a large linear system of equations (Section III-C), where each equation
connects pixel intensities in the logarithmic domain. To deal with
underexposed pixels in shorter exposures, we model heteroskedastic camera
noise with inverse-variance weights (Section III-D). Then, in Section III-F,
we show how to reduce the system of equations for faster equations without
sacrificing the estimation quality. We finally validate our exposure
estimation framework in Section IV, both on synthetic and real image captures.
††Code for the method has been integrated into the software for noise-optimal
HDR merging (_HDRutils_) and can be found at
https://github.com/gfxdisp/HDRutils.

##  II Related work

The problem of inaccurate capture parameters was identified in very early
works in HDR merging [7, 30]. However, most of these focused on the
challenging task of inverting the CRF under the assumption of film or sensor
_reciprocity_ [38]. Mitsunaga and Nayar [30] used a polynomial model to
jointly estimate the CRF and exposure times. Then, Grossberg and Nayar [15]
demonstrated how to recover the brightness transfer function (BTF), a function
that describes how the brightness transfers from one image to another. They
recovered the BTF from image histograms and stipulated that it can be used to
estimate exposure ratios if the CRF is known. More recently, Rodríguez at al.
have shown that the previously mentioned methods rely on incorrect assumptions
[35] about the independence of color channels and linearity of exposures,
resulting in estimation errors and hue shifts. We avoid those problems by
directly operating on demosaiced RAW images.

Based on [15], Cerman and Hlavac [3] assumed a linear CRF by relying on RAW
pixel values. They computed the brightness-transferred image histograms and
used them to weight a system of equations in the linear pixel domain.

In contrast to their work, we solve a weighted linear system in the
logarithmic domain to estimate exposure ratios. The weights in our system of
equations are derived from a popular statistical camera noise model and ensure
a noise-optimal solution. Our approach systematically handles heteroskedastic
camera measurements and provides accurate estimates even when many pixels are
under-exposed or affected by noise.

###  II-A Camera noise model

A key contribution of our work is the use of a parametric noise model to weigh
some pixel correspondences more than others. Our weights are based on the
popular Poisson-normal statistical camera noise model [12, 2, 28] that has
signal-dependent and static (or signal-independent) components. Many earlier
works have used approximations of similar noise models for noise-optimal HDR
reconstructions [21, 14, 19]. Although deep generative networks [1, 4] model
spatially-varying components of real camera noise better, the Poisson-normal
statistical model and its normal approximation are better suited for our
problem as they offer a convenient algebraic form under the assumption that
noise is independent at each pixel.

###  II-B HDR datasets

Merged HDR images of many multi-exposure datasets [10, 25, 24, 17, 32, 18] can
be improved with accurate exposure estimation. We observed the banding
artifacts depicted in Figure 1 in various images from these datasets. Recent
HDR deep learning reconstruction methods for tasks like HDR deghosting [24,
34, 47, 33, 5] and inverse tone-mapping [8, 37, 43] utilize these as well as
other multi-exposure datasets for testing and evaluation. All these works are
likely to produce better results when trained with accurate ground truth
information due to our work on better exposure alignment.

##  III Methodology

After introducing the camera model and some terminology, we describe how to
estimate exposure ratios directly from the input image stack by solving a
weighted linear system. We then derive noise-optimal weights for the system
based on the widely-used Poisson-normal camera noise model. Finally, we
discuss practical considerations for making the system agnostic to sensitive
noise parameters and improve robustness to ghosting caused by camera or object
motion.

###  III-A Camera model

To digitally represent an HDR scene, we start by capturing $N$ images with
varying exposure times or gains (ISO). Although earlier works included CRF
estimation in their image formation pipelines [7, 30, 15], we skip this step
since modern cameras provide access to linear RAW pixel values. We model the
captured RAW pixels as samples from independent random variables, which are
linearly related to the scene radiance, and denote them as:

$Y_{i}(p),\quad i=1{\ldots}N,\quad p=1\ldots M\,,$ (1)

where $i$ is the exposure index with $N$ total exposures, and $p$ is the pixel
index with $M$ total pixels. Throughout this work, we use upper case for
random variables and lower case for observed values. Before merging
(averaging) the exposure stack, we need to convert them to relative radiance
units by compensating for exposure time $t_{i}$, gain $g_{i}$, and aperture
f-number $a_{i}$:

$X_{i}(p)=\frac{Y_{i}(p)}{t_{i}\,g_{i}\,\pi\,\left(\frac{f}{2\,a_{i}}\right)^{2}}=\frac{Y_{i}(p)}{d_{i}}\,,$
(2)

where $f$ is the focal length. Typically, the scaling constant,
$d_{i}=t_{i}\,g_{i}\,\pi\,\left(\frac{f}{2\,a_{i}}\right)^{2}$, is computed
directly from the camera EXIF header. The problem is that $t_{i}$, $g_{i}$,
and $a_{i}$ could be incorrect due to the limited accuracy of the mechanical
shutter or the rounding of exposure values. Incorrect values of $d_{i}$ will
lead to inaccurate HDR reconstructions Thus, samples from $X_{i}(p)$, the
exposure-compensated or absolute estimates in the captured images, may
represent biased measurements of the true scene radiance. Merging images using
the inaccurately reported parameters results in the banding artifacts depicted
in Figure 1 and Figure 4. In all the images, the artifacts appear when a
longer exposure image saturates.

Figure 3: A linear gradient spanning $13$ stops is captured with exposure
times $[0.5,2,4]$ (in seconds) by simulating the _Sony ILCE-7R_ at ISO $3200$.
Each capture is quantized to $8$ bits for easy visualization. The left plot
depicts exposure-compensated pixels that represent samples from $X(p)$
(according to Eq. eq:exp-compens) with inaccurate exposures (red, yellow, and
brown dashed lines). The reconstructed gradient (blue line) is jagged around
$\phi=2^{11}$ and $\phi=2^{12}$ (green circles) due to misaligned exposures
caused by biased exposure values. We can reconstruct the smooth gradient by
correctly aligning exposure ratios, as shown in the right plot.

###  III-B Banding due to inaccurate exposure

To better understand the reason for banding, consider the one-dimensional
linear gradient depicted in Figure 3 (left). The noisy measurements (dashed
lines), obtained by simulating captures using calibrated noise parameters of
the _Sony ILCE-7R_ , become misaligned when scaled with inaccurate exposure
times. Merging such images results in a jagged reconstruction (solid blue
line), causing banding in an otherwise smooth output. Notice that although the
reconstruction deviates from the ground truth for almost all pixels, artifacts
will be visible only at transition points when an image in the stack saturates
as highlighted by the green circles.

While Figure 3 demonstrates the problem with simple averaging, banding is
further exaggerated when sophisticated algorithms based on physical noise
models [14, 21, 19] are used. This is because pixels of longer exposures are
more reliable (due to smaller noise variance), but they saturate at lower
physical values. Thus, these estimators weigh longer exposure pixels more, and
there is a sharp change upon saturation of any image in the stack.

If we are able to estimate relative exposure ratios w.r.t one of the images
correctly, we can align all exposures to produce the reconstruction in Figure
3 (right). Note that the exposure-aligned reconstruction may still not
coincide with the ground truth since we do not have the correct baseline.
However, accurately estimating relative ratios is sufficient to align the
capture parameters and eliminate banding.

###  III-C Exposure estimation in real images

To prevent banding and obtain physically accurate pixel values, we align the
exposures of all images in the stack by estimating all scaling constants
$d_{i}$. However, computing them directly from input pixels is impossible
since the correct absolute measurements representing observations of $X(p)$
from Eq. eq:exp-compens are unknown. We thus eliminate these unknowns by
estimating the ratio of exposures between any two images in the stack instead.
For a scene with constant illumination and no scene motion, this ratio should
be the same for all pixels from a given pair of images. Its expected value is:

$d_{ij}=\mathbb{E}_{p}\left[\frac{Y_{i}(p)}{Y_{j}(p)}\right]\,,$ (3)

where $i$ and $j$ index different images in the exposure stack, and
$\mathbb{E}_{p}$ indicates that the expected value is computed over all
pixels. To allow for fast computation using linear solvers, we operate on
logarithmic values. Thus, let

$e_{ij}=\log{d_{ij}}\quad\mathrm{and}\,L_{i}(p)=\log{Y_{i}(p)}\,.$ (4)

Although we cannot write a closed-form expression for the density function of
$L_{i}(p)$ because of the $\log$ transformation, it is possible to approximate
the expected value of any transformed random variable using its Taylor
expansion. Our results, detailed in Eq. eq:appendix-ev in the Appendix, are
only applicable to normally distributed random variables. It is thus
imperative to use only well-exposed pixels (we show how to do this in Section
III-F), since the Poisson photon noise component of such pixels is well-
approximated by a normal distribution. We can obtain an approximate expression
for the expected value by applying the result for $\log{Y(p)}$:

$\begin{split}e_{ij}&=\log\mathbb{E}_{p}\left[\frac{Y_{i}(p)}{Y_{j}(p)}\right]\approx\mathbb{E}_{p}\left[\log{\frac{Y_{i}(p)}{Y_{j}(p)}}\right]\\\
&=\mathbb{E}_{p}\left[L_{i}(p)-L_{j}(p)\right]\,.\end{split}$ (5)

The equality approximately holds for the operating range of pixel values, as
we will detail in Section III-D, Eq. eq:log-moments. This allows us to compute
the expected value over all pixels by setting up and solving a linear system.
For the most reliable estimate, we should utilize information from all the
available exposures (all possible values of $i$ and $j$). Unlike previous work
[3], we thus consider not only ratios between neighboring exposures but
between all pairs. This results in a large yet sparse linear system:

$\displaystyle\sqrt{\bm{W}}{\tiny\begin{bmatrix}1&-1&0&\cdots&0\\\
1&-1&0&\cdots&0\\\ \vdots&\vdots&\vdots&\ddots&\vdots\\\ 1&0&-1&\cdots&0\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\ 0&0&0&\cdots&-1\\\
\end{bmatrix}\begin{bmatrix}e_{1}\\\ e_{2}\\\ \vdots\\\
e_{N}\end{bmatrix}}=\sqrt{\bm{W}}{\tiny\begin{bmatrix}L_{1}(1)-L_{2}(1)\\\
L_{1}(2)-L_{2}(2)\\\ \vdots\\\ L_{1}(1)-L_{3}(1)\\\ \vdots\\\
L_{N-1}(M)-L_{N}(M)\\\ \end{bmatrix}}\,,$ (6)

or more compactly,

$\sqrt{\bm{W}}\bm{O}\bm{e}=\sqrt{\bm{W}}\bm{m}\,,$ (7)

where $\bm{W}$ is a diagonal weight matrix denoting the relative importance of
each row of the system. In practice, we found that the weighted system does
not always provide a good solution for shorter exposures. Since a rough
estimate of exposure values $\bm{e}_{0}$ is available in the image metadata,
we introduce a Tikhonov penalty with weight $\lambda$ and solve to get:

$\begin{split}\bm{\hat{e}}_{\mathrm{WLS}}&=\operatorname*{arg\,min}_{\bm{e}}\,\left\lVert\sqrt{\bm{W}}(\bm{O}\bm{e}-\bm{m})\right\rVert^{2}_{2}+\lambda\lVert\bm{e}-\bm{e}_{0}\rVert^{2}_{2}\\\
&=(\bm{O}^{\mathsf{T}}\bm{W}\bm{O}+\lambda\bm{I})^{-1}(\bm{O}^{\mathsf{T}}\bm{W}\bm{m}+\lambda\bm{e}_{0})\,.\end{split}$
(8)

(a) Parameters from EXIF

(b) Estimated parameters

Figure 4:  The left column shows exposures (gamma-encoded for visualization,
$\gamma=2.2$) of the HDR image reconstructed with EXIF parameters, while the
right column shows exposures (encoded with $\gamma=2.2$) of the HDR image
reconstructed by solving the linear system represented by Eq. eq:exp-sparse-
lin with $\bm{W}=\bm{I}$. Banding artifacts are visible in dark (top row) and
bright (bottom row) regions of the reconstruction using EXIF metadata (a).
Aligning exposures according to Eq. eq:lsqr fixes the problem (b).

To demonstrate the effectiveness of exposure alignment, we captured an
exposure stack of a simple HDR scene consisting of a bright light shining at
an angle, which produces an inverse-square fall-off in intensity. Since most
pixels are well-exposed, as shown in Figure 4, we could assume that they are
equally reliable and set $\bm{W}=\bm{I}$. This works well for the carefully
controlled scene and eliminates banding artifacts that would have appeared
when merging with EXIF parameters. However, the constant noise assumption
breaks down for real-world HDR stacks since the noise in camera pixels is
heteroskedastic [11] and thus, different pixels provide different amounts of
information.

###  III-D Heteroskedastic pixels

To determine noise-optimal weights for images of real cameras, we need to
derive an expression for the variance of each row of Eq. eq:exp-sparse-lin.
Then, we populate $\bm{W}$ with inverse-variance weights (i.e., each weight is
given by the reciprocal of variance) to obtain the weighted least-square
estimate $\bm{\hat{e}_{\mathrm{WLS}}}$. This is equivalent to maximum
likelihood estimation (MLE) under the assumption that the system of equations
models additive, normally distributed noise, where the variance is different
for different rows. Since each entry of the output vector $\bm{m}$ is a
difference of two random variables, the inverse-variance weight for each row
of the system is given by:

$w_{k,k}=\frac{1}{\mathbb{V}[L_{i}(p)-L_{j}(p)]}=\frac{1}{\mathbb{V}[L_{i}(p)]+\mathbb{V}[L_{j}(p)]}\,.$
(9)

Here, $k$ indexes each row of Eq. eq:exp-sparse-lin since $\bm{W}$ is
diagonal. For instance, the second row corresponds to $k=2,i=1,j=2$.

In order to get an expression for the denominator in Eq. eq:w-inv-var, we
refer to detailed studies of the noise characteristics of cameras [2, 28],
which indicate that real camera noise follows a compound Poisson-normal
distribution. For the working range of pixels in most images, this can be
approximated by zero-mean additive noise that follows a normal distribution.
Temporarily dropping the exposure index for brevity of notation, the variance
at each pixel thus consists of a signal-dependent component as well as a
static component and is equal to,

$\sigma_{Y}^{2}(p)=\mathbb{V}[Y(p)]=\alpha\,\mathbb{E}[Y(p)]+\beta\,,$ (10)

where $\alpha$ and $\beta$ are camera-specific noise parameters, which also
depend on the sensor’s gain.

Since $L(p)$ is a random variable obtained by applying the $\log$
transformation to $Y(p)$, its density function does not have an exact
expression. This is because the domain of $Y(p)$ includes negative values for
which the $\log$ function is undefined. As we will show in Section III-F, we
operate on a small subset of available pixels, selected for the lowest
relative variances, and whose intensities tend to be much greater than $0$. We
show, in Eq. eq:appendix-var in the Appendix, how to approximate the variance
when any random variable is transformed by an invertible function. Here, we
apply the result for $Y(p)$ when it is transformed by the $\log$ function,

$\begin{split}&\mathbb{E}[L(p)]\approx\log{\mu_{Y}(p)}\,,\\\
&\mathbb{V}[L(p)]\approx\frac{\sigma^{2}_{Y}(p)}{\mu^{2}_{Y}(p)}\,,\end{split}$
(11)

where $\sigma^{2}_{Y}(p)$ is the variance of the noise from Eq. eq:noise-model
and $\mu^{2}_{Y}(p)$ is the expected pixel value, which we approximate by
$\mathbb{E}[Y(p)]=\mu_{Y}(p)\approx y(p)$. We reiterate that we are able to
use this result because we work with well-exposed pixels, that are
approximately normal. After substituting the computed variance in Eq.
eq:w-inv-var, the diagonal weights become:

$w_{k,k}=\left(\frac{\alpha y_{i}(p)+\beta}{y_{i}^{2}(p)}+\frac{\alpha
y_{j}(p)+\beta}{y_{j}^{2}(p)}\right)^{-1}\,.$ (12)

###  III-E Camera noise calibration

While the inverse-variance weights in Eq. eq:weights help us compute noise-
optimal estimates of the exposure ratios, a fundamental limitation is their
dependence on calibrated noise parameters $\alpha$ and $\beta$. In the
considered noise model [12], these are camera and gain specific and may not be
available at the time of HDR merging. Moreover, the quality of HDR
reconstructions is highly sensitive to accurate noise parameters [2],
motivating the need for methods that do not rely on them.

If noise parameters are unavailable or inaccurate, it is still possible to
solve the weighted linear system by assuming that the static noise parameter
$\beta$ is $0$ as we work with well-exposed pixels (see Section III-F). Since
the entries of $\bm{W}$ determine the relative importance of the different
rows of Eq. eq:exp-sparse-lin, we can eliminate the common signal-dependent
constant, $\alpha$ as well. The new weights which no longer depend on
calibration-sensitive parameters are:

$w^{\prime}_{k,k}=\left(\frac{1}{y_{i}(p)}+\frac{1}{y_{j}(p)}\right)^{-1}\,.$
(13)

In practice, we can get away with using camera-independent $\bm{W^{\prime}}$
instead of camera-specific $\bm{W}$ from Eq. eq:weights.

###  III-F Reducing the linear system

The linear system given by Eq. eq:exp-sparse-lin contains up to
$M\binom{N}{2}$ equations corresponding to all pixels ($M$) and all pairs of
exposures in the stack ($N$). Solving such a large system may be impossible
for large images or deep exposure stacks due to computation and memory
limitations. However, this system is strongly overdetermined as only $N-1$
exposure ratios need to be estimated. Therefore, we only need a small
percentage of equations to solve Eq. eq:lsqr. Another reason for reducing the
system is to eliminate underexposed pixels pairs because they do not satisfy
the assumptions made in the derivations of previous sections. For example,
dark pixels may be negative or poorly approximated by a normal distribution.

A logical reduction strategy is to select pixels pairs with the highest
weights since they are least affected by noise. However, such a strategy is
heavily biased towards longer exposures because pixel pairs that include these
images will have the highest weights. The shorter exposures will then be
poorly represented and, thus, poorly estimated. Another issue is related to
the spatial location of bright objects in the scene, such as the Sun or other
light sources. If we select a small fraction of pixels (say 5% of the pixels
per image), all of them are likely to be concentrated in one portion of the
image, corresponding to these bright objects. If those objects happen to be in
motion, the exposure estimation will fail.

We propose two orthogonal design choices to balance the system of equations
and handle both these biases in estimated exposure ratios.

####  III-F1 Spatial balance: Tiling

To ensure that the linear system contains samples from all parts of a scene,
we split the input image stack into $t\times t$ pixel tiles ($16\times 16$ in
our experiments). We can then select a fixed number of pixel pairs from each
tile and pixels from a few bright objects will not dominate the system of
equations.

The tiled processing provides a convenient way to vectorize the construction
of the reduced linear system. Several tiles can be processed in parallel for
faster execution with multi-core or multi-threaded systems.

Our noise-based solution helps provide a robust estimation even though some
tiles may contain noisy pixels corresponding to dimly-lit parts of the scene.
Variance-optimal balancing of exposures is crucial for such tiles, as we will
show in the next section.

####  III-F2 Exposures balance: Spanning trees

Figure 5: Consider the exposure multigraph in (a) where different colored
vertices and edges represent different pixels. The MST may connect pairwise
exposures (b) or connect the longest exposure to all others (c). Algorithm 1
will provide different solutions for different inputs. For example, the MST
for the blue horizontal tile near the Sun in (d) is the pairwise-connected
spanning tree connected at pixels marked by the two triangles in plot (e).
However, for the orange tile, the MST is given by the longest-exposure
spanning tree at the pixel marked by the star in plot (e).

The next objective is guaranteeing that all $N-1$ exposure ratios are
correctly estimated. Within each spatial tile, we need to include pixels from
all exposure pairs, including noisy short exposures despite their relatively
smaller weights.

Consider the _exposure graph_ : an undirected weighted multigraph (a graph
that contains more than one edge between two vertices as shown in Figure 5a),
whose vertices represent the $N$ exposures and edges link pairs of co-located
pixels from two exposures. Each edge corresponds to one row of Eq. eq:exp-
sparse-lin, with its weight given by Eq. eq:weights or Eq. eq:weights-prime.
The different colored vertices and edges in Figure 5a represent different
pixel locations. Reducing the linear system from Eq. eq:exp-sparse-lin is then
equivalent to removing edges from this dense multigraph.

The optimal subset contains $(N-1)k$ edges with the largest weights that
connect all exposures in a balanced manner. Such a subset can be found by
computing the $k$ highest-weighted spanning trees of the multigraph, where the
weight of a spanning tree is the sum of the weights of its edges. The previous
work directly linked pixels from neighboring exposures [3], resulting in
pairwise connectivity as shown in Figure 5b. However, this solution is sub-
optimal for some inputs (such as the orange tile in Figure 5d) because using
edges linked to the longest exposures, as shown in Figure 5c), results in
higher weights and, therefore better estimates.

An optimal solution would sequentially extract $k$ MSTs using _Kruskal’s_ or
_Prim’s_ algorithm. Better algorithms [26, 9] extract the $k$ highest spanning
trees more efficiently. However, explicitly creating a large multigraph and
computing MSTs is both memory and computationally expensive. Below we show
that the optimal solution can be found by a simpler greedy algorithm.

##### Greedy MST solution

input : $y_{1},\ldots,y_{N}$ (Images sorted by exposure time)

output : MST (Maximimum spanning Tree)

1 Function  _greedyMST(_$y[\;]$_)_

2 $N\leftarrow length(y)$

3 $\textrm{MST}\leftarrow[\;]$

4 for  _$i\leftarrow 1$ to $N-1$_ do

/* Iterate over all exposures starting from the shortest */

5 $mask\leftarrow$ isValid(_$y_{i}$_) and isValid(_$y_{i+1}$_)

$p_{*}\leftarrow$ maxWeight(_$y_{i}[mask],y_{i+1}[mask]$_)

// Location of the highest weighted edge connecting images $y_{i}$ and
$y_{i+1}$

6 for  _$j\leftarrow N$ to $i+1$_ do

7 if  _isValid(_$y_{j}(p_{*})$_)_ then

8
$\textrm{MST}.$addEdge(_$i\mathrel{\bullet\mkern-3.0mu{-}\mkern-3.0mu\bullet}j$_)

9 break

10

11

12

13 return _MST_

Algorithm 1 Greedy MST algorithm

We start by extracting all valid pixels — those that are unsaturated and
sufficiently above the noise floor. We iterate through all exposures from the
shortest to the longest. For each exposure $i$, we identify $p_{*}$, the pixel
location that contains the highest-weighted edge between exposures $i$ and
$i+1$. This is typically the brightest pixel that is not saturated in exposure
$i+1$. Then, we find the longest exposure $j>i$ in which pixel $p_{*}$ is not
saturated. We add an edge between $i$ and $j$. By selecting the longest
exposure, we ensure that the weight of the edge is maximized.

This procedure, summarized in Algorithm 1, can be repeated $k$ times to
extract the $k$ highest-weighted spanning trees without explicitly
constructing a graph. If all longer exposures
$y_{N}(p_{*}),y_{N-1}(p_{*}),\ldots y_{i+2}(p_{*})$ are invalid due to
saturation, the solution reduces to pairwise connectivity depicted in Figure
5b.

Note that while a pair of pixels forming an edge must have the same position,
each edge in the MST can come from a different pixel position. Consider the
1-dimensional blue tile in the tonemapped image in Figure 5d and its
horizontal scanline in Figure 5e. The MST is given by two edges: an edge
between the medium and long exposures and an edge between the medium and short
exposures (marked by triangles in the plot). Note that these are the brightest
unsaturated pixels in the long and medium exposures. The MST for the blue tile
is thus the pairwise connected spanning tree. The limited dynamic range of the
orange tile means that a single pixel, marked by the star in Figure 5e, aligns
all three exposures for the orange tile, resulting in the spanning tree in
which the longest exposure is linked to all other exposures.

###  III-G Handling outliers and pixel misalignment

(a) Regular $\bm{\hat{e}_{\mathrm{WLS}}}$

(b) $\bm{\hat{e}_{\mathrm{WLS}}}$ with outlier removal

Figure 6: For an exposure stack containing some spatially misaligned pixels,
solving the weighted linear system Eq. eq:wls results in inaccurate estimates.
In this scene, camera motion causes pixels of the building to be misaligned
due to camera motion, producing artifacts near the windows in Figure 6a.
However, pixels from other areas, such as the sky, can still be utilized for
exposure estimation. Figure 6b demonstrates the advantage of using the
outlier-removal procedure described in Section III-G to improve robustness.

The weighted least-squares solution Eq. eq:lsqr using rows of the reduced
system will provide an accurate estimate only if there is no movement in the
scene and all pixels are well aligned. If there is motion across exposures,
for example, due to camera shake or object motion, those pairs of pixels will
affect our estimates of exposure ratios, as shown in Figure 6. Notice the
artifacts surrounding bright regions such as the windows and red lights in
Figure 6a due to camera motion. Here, we estimated the exposures using the
original stack and then registered the images with homography alignment [39]
before merging. Thus, any visible artifacts are due to incorrect exposure
alignment and not due to ghosting.

The spatial tiling described in Section III-F1 ensures that motion in a few
bright objects does not adversely affect the estimated exposures, i.e., our
procedure is robust to localized object motion. In case of widespread pixel
misalignment (e.g., due to camera motion), we can still provide reasonable
exposure estimates if the scene contains uniform image patches. We identify
_usable_ tiles by collecting $k$ MSTs within a tile and solving the smaller
per-tile system of equations separately. If the predictions of a tile deviate
too much from the EXIF values, we treat the tile as an outlier and do not
include the corresponding MSTs in the final system.

The added overhead of solving per-tile linear systems is negligible because
each system contains only $k$ equations ($k=50$ in our experiments). This
outlier removal is much faster than an iterative algorithm such as the one
described in [3].

##  IV Results and applications

The primary application of our method is merging exposure stacks with varying
exposure times or gains. We first validate our results on synthetically
generated stacks, for which we know the ground truth, and then show
qualitative comparisons on real captures. Finally, in Section IV-D, we show
how our method can improve the estimate of the MTF of a display. All our
results use $\lambda=10$ for the Tikhonov penalty term described in Eq.
eq:lsqr. For merging exposures stacks after exposure estimation, we used a
noise-aware HDR estimator [19]. However, our results hold for other simpler
methods too.

###  IV-A Validation on synthetic dataset

Figure 7: Accuracies of exposure estimation methods with increasing ISO
computed by simulating captures on an existing HDR dataset [10]. Error bars
indicate $95\%$ confidence intervals computed across all images. The blue line
represents corrupted exposures according to Eq. eq:res-error, while the
baseline, in green, was obtained by directly estimating exposures by solving a
system based on Eq. eq:expectation-ratio.

First, we rely on synthetic exposure stacks to compare the accuracy of
different methods for exposure estimation. We simulated HDR exposure stacks
using the noise parameters of the _Canon PowerShot S100_ with exposure times
[$\nicefrac{{1}}{{64}}$, $\nicefrac{{1}}{{8}}$, 1, 8] (in seconds). All images
were quantized to 14 bits to match the bit-depth of the camera, and the
experiment was repeated for ISO settings between 100 and 800. We simulated
noise according to Eq. eq:noise-model, with noise parameters listed in Table
I.††Code for noise simulation is available at
https://github.com/gfxdisp/HDRutils/tree/main/HDRutils/noise_modeling.

TABLE I: Noise parameters measured from the _Canon PowerShot S100_ at various ISOs. These parameters have been calculated after normalizing the camera sensor values to a range of 0 to 1. ISO | $\alpha_{\textrm{R}}$ | $\alpha_{\textrm{G}}$ | $\alpha_{\textrm{B}}$ | $\beta_{\textrm{R}}$ | $\beta_{\textrm{G}}$ | $\beta_{\textrm{B}}$
---|---|---|---|---|---|---
100 | 2.46e-5 | 1.67e-5 | 7.41e-5 | 3.58e-8 | 2.13e-8 | 1.28e-7
200 | 4.57e-5 | 3.02e-5 | 1.32e-4 | 9.89e-8 | 6.07e-8 | 2.66e-7
400 | 9.12e-5 | 5.95e-5 | 2.59e-4 | 2.21e-7 | 1.72e-7 | 5.61e-7
800 | 1.85e-4 | 1.19e-4 | 5.26e-4 | 4.94e-7 | 4.28e-7 | 1.14e-6

Source HDR images were taken from the Fairchild photographic survey [10]
containing $105$ scenes at a resolution of $4312\times 2868$. All simulated
captures are linearly related to the HDR image, i.e., we do not apply a CRF
since it is likely to introduce artifacts.

Before merging the images, we corrupted the exposure times by introducing a
small amount of normally distributed noise:

$e^{\prime}_{i}=e_{i}+\eta_{i}\quad\mathrm{where}\,\,\eta_{i}\sim\mathcal{N}(0,0.15\,e_{i})\,.$
(14)

The relative standard deviation factor $0.15$ was selected to match real
camera EXIF errors plotted in Figure 2.

Figure 7 plots the relative root-mean-squared error (RMSE) (in percent) of the
exposure ratios, for different ISO levels. Note that we can compute such
relative errors only for synthetic datasets. We include results of a simple
baseline by solving a linear system that realises Eq. eq:expectation-ratio
without relying on a noise model. Note that the baseline includes other
components described in Section III, such as spatial balance via tiling and
exposure balance via MSTs. The blue line shows the error of noisy exposures
according to Eq. eq:res-error. The synthetic exposure stacks represent a
challenging estimation problem with exposures three stops apart. As a result,
the simple baseline (green) and the histogram-based BTF method [3] (orange)
are unable to provide good results due to the adverse impact of pixels in
shorter exposures. A key limitation of both methods is not explicitly modeling
camera noise, which we address in our formulation.

Figure 7 depicts results for the two versions of our weighted least-squares
solution (Eq. eq:wls), showing both pairwise connectivity (in red) and the
greedy MST heuristic (in purple). Both our solutions use weights that model
the noise characteristics of the camera. Thus, they result in better estimates
with smaller errors. Overall, the greedy MST heuristic is the best-performing
estimator of exposure ratios.

Execution times: Histogram matching used in [3] is computationally expensive
for high-resolution images, resulting in an average execution time of 2.11
seconds. Further, we observed that the expensive iterative procedure for
removing outliers is needed for good estimates increasing time to 8.54
seconds. Our tiled reduction (Section III-F1) results in a significantly
faster average execution time of 0.265 seconds for pairwise connected
exposures and 0.29 seconds for the greedy MST solution. The reported times
were computed on an _Intel i7-8700 CPU_ for images of resolution $4312\times
2868$.

###  IV-B Performance on real captures

Figure 8: Zoomed in patches of light sources, appropriately exposed to
highlight banding close to the source due to inaccurately reported EXIF
exposures (first row). The baseline method (second row) that does not model
camera noise is unable to estimate shorter exposures, resulting in artifact-
ridden reconstructions. Similarly, only some exposures can be aligned with
histogram weights based on BTF (third row). By modeling camera noise, we
simultaneously align all exposures (last row).

Next, we show the advantage of our proposed exposure estimation over naive
merging using EXIF data and compare it with the histogram-based BTF approach
[3]. The RAW images have a bit-depth of 14 and are linearly related to the
scene radiance. We did not apply CRF or tone-mapping for any image shown in
this section. Thus, all artifacts visible are due to incorrect exposures. For
all results in this section, we use calibration-free weights given by Eq.
eq:weights-prime.

Banding artifacts are visible for scenes containing a smooth gradient at the
transition point when one of the input exposures saturates. This frequently
occurs very close to bright light sources such as the Sun during the day or
street lights at night, as indicated by red arrows in Figure 8 (a). Pixels
close to light sources tend to be unsaturated only in shorter exposures where
other parts of the image are strongly affected by noise. Aligning exposures
under these conditions is challenging, resulting in the poor performance of
the baseline as well as the histogram-based BTF weights [3]. The baseline
method obtained by directly solving Eq. eq:expectation-ratio completely fails
to recover exposure ratios. When using the BTF weights, the banding artifacts
also persist at the same locations or are introduced in other locations as
shown in Figure 8 (c). By accounting for heteroskedastic noise with inverse-
variance weights in Figure 8 (d), we can correctly align all the exposures to
produce banding-free results.

Even in the absence of point light sources, banding can appear at natural
smooth gradients or due to defocus blur, as shown in Figure 9. Large regions
of the images are well-exposed in such scenes, and most methods work
reasonably well. However, the baseline and BTF weights can still sometimes
fail to recover the exposure ratio for shorter exposures (for example, see the
red patch in the third column and green and blue patches in the second column
of the first scene). Our noise-model motivated approach consistently aligns
all exposures to produce banding-free reconstructions across the scenes.

Figure 9: Fixing banding in smooth image gradients when images contain
sufficient well-exposed pixels to populate the linear system. Our inverse-
variance solution (last column) is more robust and produces better
reconstructions.

###  IV-C Deghosting for scenes with motion

(a) Aligned
exposures

(b) Corrupted
exposures

(c) Reference

Figure 10: Comparison of HDR reconstructions of an attention-based deep
network [42] with aligned (first column) and corrupted (middle column)
exposures. These insets are from the HDR deghosting dataset [24], which the
deghosting network was trained on.

When multi-image stacks with misaligned exposures are fused, the resulting HDR
reconstructions often exhibit accentuated ghosting artifacts. Even recent deep
learning-based methods [24, 45, 42] are unable to account for inaccurate
exposure ratios. Figure 10b illustrates the artifacts that arise when
exposures are corrupted according to Eq. eq:res-error. It is important to note
that these artifacts are in addition to the ghosting artifacts visible in
Figure 10a. Our proposed algorithm, which includes tile-based outlier
detection to account for scene motion, effectively aligns exposures.

###  IV-D Measuring display MTF

Figure 11: Edge profile of an HDR display recovered from a slant edge captured
with an exposure stack. Errors in EXIF metadata lead to deviation from the
actual edge spread function, affecting MTF estimation. On inspecting the
captured images, we observed the banding artifacts close to the captured
slant-edge when the image was merged with EXIF exposures.

Digital cameras are often used as inexpensive light measuring instruments, for
example, for measuring the spatial characteristic of an electronic display,
such as spatial uniformity of a modulation transfer function (MTF) [48].
Inaccuracies in pixel values due to exposure misalignment may lead to errors
in the estimated display characteristics. Here, we demonstrate the need for
exposure alignment for measuring MTF of HDR displays.

The MTF measurement involves taking an image of a slanted edge [36], where one
side has its pixels set to 0 and the other side to the maximum pixel value.
Modern displays can reach very high contrast levels, so such a slanted-edge
image can be faithfully captured only using an exposure stack. The
inaccuracies of exposure times can easily introduce bias to our measurements.

Figure 11 shows the edge profile computed with (blue) and without (orange)
exposure estimation. The difference in reconstructions is due to inaccurate
EXIF exposure values, visible as banding in the merged image. This translates
to the bias at low-radiance pixels in Figure 11 and shows that errors in
exposure times can easily lead to biases in the final MTF estimates.

##  V Conclusions

We propose to estimate relative exposures directly from images of an HDR
exposure stack instead of relying on inaccurate camera EXIF metadata. The
problem can be formulated as a weighted linear system of equations with
Tikhonov regularization. Our formulation considers the noise characteristics
of a camera to derive weights. We show that good performance on a large number
of scenes is possible without the need for camera- and gain-specific noise
parameters.

We also describe an efficient approach based on the union of the highest
weighted spanning trees of the exposure multigraph to reduce the size of the
system for large images. When applied to multi-image HDR stacks, our method
eliminates banding artifacts at smooth image gradients close to light sources
or due to defocus blur. The exposure-aligned images are closer to physical
quantities, making our work essential for using a camera as a light-
measurement instrument.

## Acknowledgments

This project has received funding from the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation programme
(grant agreement N∘ 725253–EyeCode).

[Moments of a transformed random variable] Here, we use the Taylor expansion
to derive expressions for the first and second moments (expected value and
variance respectively) of a random variable under an arbitrary invertible
function. In Section III-D we apply these results to obtain expressions for
the $\log$ of a normally distributed random variable.

To get an approximate variance for a nonlinear function of a random variable,
we describe the Delta method [31, 23] that uses a Taylor expansion. Let $Z$ be
an asymptotically normal random variable with a known mean $\mu_{Z}$ and
variance $\sigma_{Z}^{2}$, and let $f$ be an invertible and differentiable
function. To get expressions for the expected value and variance of $f(Z)$ we
use its first order Taylor expansion about the mean $\mu_{Z}$,

$\begin{split}f(Z)=\,&f(\mu_{Z})+f^{\prime}(\mu_{Z})(Z-\mu_{Z})\\\
&+\textrm{higher-order terms}\,.\end{split}$ (15)

Ignoring the diminishing contributions of higher-order terms, the expected
value is,

$\begin{split}\mathbb{E}[f(Z)]&\approx\mathbb{E}[f(\mu_{Z})]+\mathbb{E}[f^{\prime}(\mu_{Z})(Z-\mu_{Z})]\\\
&=\mathbb{E}[f(\mu_{Z})]+f^{\prime}(\mu_{Z})(\mathbb{E}[Z]-\mathbb{E}[\mu_{Z})])\\\
&=f(\mu_{Z})+f^{\prime}(\mu_{Z})(\mu_{Z}-\mu_{Z})\\\
&=f(\mu_{Z})\,,\end{split}$ (16)

and the variance is

$\begin{split}\mathbb{V}[f(Z)]&\approx\mathbb{V}[f(\mu_{Z})]+\mathbb{V}[f^{\prime}(\mu_{Z})(Z-\mu_{Z})]\\\
&=\mathbb{V}[f(\mu_{Z})]+f^{\prime}(\mu_{Z})^{2}\,\mathbb{V}[Z-\mu_{Z}]\\\
&=\mathbb{V}[f(\mu_{Z})]+f^{\prime}(\mu_{Z})^{2}\,(\mathbb{V}[Z]+\mathbb{V}[\mu_{Z}])\\\
&=0+f^{\prime}(\mu_{Z})^{2}\,(\sigma_{Z}^{2}+0)\\\
&=f^{\prime}(\mu_{Z})^{2}\,\sigma_{Z}^{2}\,.\end{split}$ (17)

## References

  * [1] A. Abdelhamed, M. A. Brubaker, and M. S. Brown, “Noise Flow: Noise Modeling with Conditional Normalizing Flows,” in _International Conference on Computer Vision (ICCV)_ , 2019.
  * [2] C. Aguerrebere, J. Delon, Y. Gousseau, and P. Musé, “Study of the digital camera acquisition process and statistical modeling of the sensor raw data,” Laboratoire Traitement et Communication de l’Information, Instituto de Ingeniería Eléctrica, Tech. Rep., Sep. 2012. [Online]. Available: https://hal.archives-ouvertes.fr/hal-00733538
  * [3] L. Cerman and V. Hlavac, “Exposure time estimation for high dynamic range imaging with hand held camera,” in _Proc. of Computer Vision Winter Workshop, Czech Republic_. Citeseer, 2006\.
  * [4] K.-C. Chang, R. Wang, H.-J. Lin, Y.-L. Liu, C.-P. Chen, Y.-L. Chang, and H.-T. Chen, “Learning camera-aware noise models,” in _Proceedings of European Conference on Computer Vision (ECCV)_ , 2020.
  * [5] J. Chen, Z. Yang, T. N. Chan, H. Li, J. Hou, and L.-P. Chau, “Attention-guided progressive neural texture fusion for high dynamic range image restoration,” _IEEE Transactions on Image Processing_ , vol. 31, pp. 2661–2672, 2022.
  * [6] Y. Chen, G. Jiang, M. Yu, Y. Yang, and Y.-S. Ho, “Learning stereo high dynamic range imaging from a pair of cameras with different exposure parameters,” _IEEE Transactions on Computational Imaging_ , vol. 6, pp. 1044–1058, 2020\.
  * [7] P. E. Debevec and J. Malik, “Recovering high dynamic range radiance maps from photographs,” in _Proceedings of the 24th Annual Conference on Computer Graphics and Interactive Techniques_ , ser. SIGGRAPH ’97. USA: ACM Press/Addison-Wesley Publishing Co., 1997, p. 369–378. [Online]. Available: https://doi.org/10.1145/258734.258884
  * [8] G. Eilertsen, J. Kronander, G. Denes, R. Mantiuk, and J. Unger, “Hdr image reconstruction from a single exposure using deep cnns,” _ACM Transactions on Graphics (TOG)_ , vol. 36, no. 6, 2017.
  * [9] D. Eppstein, “Finding the k smallest spanning trees,” in _Scandinavian Workshop on Algorithm Theory_. Springer, 1990, pp. 38–47.
  * [10] M. D. Fairchild, “The hdr photographic survey,” in _Color and imaging conference_ , vol. 2007. Society for Imaging Science and Technology, 2007, pp. 233–238.
  * [11] A. Foi, “Clipped noisy images: Heteroskedastic modeling and practical denoising,” _Signal Processing_ , vol. 89, no. 12, pp. 2609–2629, 2009, special Section: Visual Information Analysis for Security. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0165168409001996
  * [12] A. Foi, M. Trimeche, V. Katkovnik, and K. Egiazarian, “Practical poissonian-gaussian noise modeling and fitting for single-image raw-data,” _IEEE Transactions on Image Processing_ , vol. 17, no. 10, pp. 1737–1754, 2008.
  * [13] O. Gallo and P. Sen, “Stack-Based Algorithms for HDR Capture and Reconstruction,” in _High Dynamic Range Video_. Elsevier, 2016, pp. 85–119.
  * [14] M. Granados, B. Ajdin, M. Wand, C. Theobalt, H.-P. Seidel, and H. P. Lensch, “Optimal hdr reconstruction with linear digital cameras,” in _2010 IEEE Computer Society Conference on Computer Vision and Pattern Recognition_. IEEE, 2010, pp. 215–222.
  * [15] M. D. Grossberg and S. K. Nayar, “What can be known about the radiometric response from images?” in _Computer Vision — ECCV 2002_ , A. Heyden, G. Sparr, M. Nielsen, and P. Johansen, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg, 2002, pp. 189–205.
  * [16] S. Hajisharif, J. Kronander, and J. Unger, “Adaptive dualiso hdr reconstruction,” _EURASIP Journal on Image and Video Processing_ , vol. 2015, no. 1, p. 41, 2015.
  * [17] P. Hanji, M. Z. Alam, N. Giuliani, H. Chen, and R. K. Mantiuk, “Hdr4cv: High dynamic range dataset with adversarial illumination for testing computer vision methods,” _Journal of Imaging Science and Technology_ , 2021. [Online]. Available: http://www.cl.cam.ac.uk/research/rainbow/projects/hdr4cv-dataset/
  * [18] P. Hanji, R. K. Mantiuk, G. Eilertsen, S. Hajisharif, and J. Unger, “Comparison of single image hdr reconstruction methods — the caveats of quality assessment,” in _Special Interest Group on Computer Graphics and Interactive Techniques Conference Proceedings (SIGGRAPH ’22 Conference Proceedings)_ , 2022. [Online]. Available: https://www.cl.cam.ac.uk/research/rainbow/projects/sihdr_benchmark/
  * [19] P. Hanji, F. Zhong, and R. K. Mantiuk, “Noise-aware merging of high dynamic range image stacks without camera calibration,” in _Advances in Image Manipulation (ECCV workshop)_. Springer, 2020, pp. 376–391. [Online]. Available: https://doi.org/10.1007/978-3-030-67070-2_23
  * [20] S. W. Hasinoff, D. Sharlet, R. Geiss, A. Adams, J. T. Barron, F. Kainz, J. Chen, and M. Levoy, “Burst photography for high dynamic range and low-light imaging on mobile cameras,” _ACM Transactions on Graphics (TOG)_ , vol. 35, no. 6, p. 192, 2016.
  * [21] S. Hasinoff, F. Durand, and W. Freeman, “Noise-optimal capture for high dynamic range photography,” in _CVPR_. IEEE, 2010, pp. 553–560.
  * [22] F. Heide, M. Steinberger, Y.-T. Tsai, M. Rouf, D. Pajak, D. Reddy, O. Gallo, J. Liu, W. Heidrich, K. Egiazarian _et al._ , “Flexisp: A flexible camera image processing framework,” _ACM Transactions on Graphics (TOG)_ , vol. 33, no. 6, pp. 1–13, 2014.
  * [23] J. M. V. Hoef, “Who invented the delta method?” _The American Statistician_ , vol. 66, no. 2, pp. 124–127, 2012. [Online]. Available: https://doi.org/10.1080/00031305.2012.687494
  * [24] N. K. Kalantari and R. Ramamoorthi, “Deep high dynamic range imaging of dynamic scenes,” _ACM Transactions on Graphics (Proceedings of SIGGRAPH 2017)_ , vol. 36, no. 4, 2017.
  * [25] K. Karaduzovic-Hadziabdic, J. H. Telalovic, and R. Mantiuk, “Subjective and Objective Evaluation of Multi-exposure High Dynamic Range Image Deghosting Methods,” in _Eurographics 2016 - Short Papers_ , 2016.
  * [26] N. Katoh, T. Ibaraki, and H. Mine, “An algorithm for finding k minimum spanning trees,” _SIAM Journal on Computing_ , vol. 10, no. 2, pp. 247–255, 1981.
  * [27] M. Kim, M. Azimi, and R. K. Mantiuk, “Perceptually motivated model for predicting banding artefacts in high-dynamic range images,” _Color and Imaging Conference_ , vol. 2020, no. 28, pp. 42–48, nov 2020. [Online]. Available: https://www.ingentaconnect.com/content/10.2352/issn.2169-2629.2020.28.8
  * [28] M. Konnik and J. Welsh, “High-level numerical simulations of noise in ccd and cmos photosensors: review and tutorial,” _arXiv preprint arXiv:1412.4031_ , 2014.
  * [29] K. Ma, Z. Duanmu, H. Yeganeh, and Z. Wang, “Multi-exposure image fusion by optimizing a structural similarity index,” _IEEE Transactions on Computational Imaging_ , vol. 4, no. 1, pp. 60–72, 2018.
  * [30] T. Mitsunaga and S. Nayar, “Radiometric self calibration,” in _Proceedings. 1999 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (Cat. No PR00149)_ , vol. 1, 1999, pp. 374–380 Vol. 1\.
  * [31] G. W. Oehlert, “A note on the delta method,” _The American Statistician_ , vol. 46, no. 1, pp. 27–29, 1992. [Online]. Available: https://www.tandfonline.com/doi/abs/10.1080/00031305.1992.10475842
  * [32] E. Perez-Pellitero, S. Catley-Chandar, A. Leonardis, and R. Timofte, “Ntire 2021 challenge on high dynamic range imaging: Dataset, methods and results,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR) Workshops_ , June 2021, pp. 691–700.
  * [33] K. R. Prabhakar, S. Agrawal, and R. V. Babu, “Self-gated memory recurrent network for efficient scalable hdr deghosting,” _IEEE Transactions on Computational Imaging_ , vol. 7, pp. 1228–1239, 2021.
  * [34] Z. Pu, P. Guo, M. S. Asif, and Z. Ma, “Robust high dynamic range (hdr) imaging with complex motion and parallax,” in _Proceedings of the Asian Conference on Computer Vision (ACCV)_ , November 2020.
  * [35] R. G. Rodríguez, J. Vazquez-Corral, and M. Bertalmío, “Issues with common assumptions about the camera pipeline and their impact in hdr imaging from multiple exposures,” _https://doi.org/10.1137/19M1250248_ , vol. 12, pp. 1627–1642, 10 2019. [Online]. Available: https://epubs.siam.org/doi/10.1137/19M1250248
  * [36] E. Samei, M. J. Flynn, and D. A. Reimann, “A method for measuring the presampled mtf of digital radiographic systems using an edge test device,” _Medical physics_ , vol. 25, no. 1, pp. 102–113, 1998.
  * [37] M. S. Santos, R. Tsang, and N. Khademi Kalantari, “Single image hdr reconstruction using a cnn with masked features and perceptual loss,” _ACM Transactions on Graphics_ , vol. 39, no. 4, 7 2020.
  * [38] T. Tani, _Photographic sensitivity: theory and mechanisms_. Oxford University Press on Demand, 1995, no. 8.
  * [39] A. Tomaszewska and R. Mantiuk, “Image registration for multi-exposure high dynamic range image acquisition,” in _Proceedings of the International Conference in Central Europe on Computer Graphics, Visualization and Computer Vision_ , 01 2007.
  * [40] Z. Tu, J. Lin, Y. Wang, B. Adsumilli, and A. C. Bovik, “Bband index: A no-reference banding artifact predictor,” in _ICASSP 2020 - 2020 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , 2020, pp. 2712–2716.
  * [41] O. T. Tursun, A. O. Akyüz, A. Erdem, and E. Erdem, “The state of the art in hdr deghosting: a survey and evaluation,” in _Computer Graphics Forum_ , vol. 34, no. 2. Wiley Online Library, 2015, pp. 683–707.
  * [42] T. V.V, “High dynamic range image synthesis via attention non-local network,” 2020\. [Online]. Available: github.com/tuvovan/ANL-HDRI
  * [43] L. Wang and K.-J. Yoon, “Deep learning for hdr imaging: State-of-the-art and future trends,” _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 2021.
  * [44] Y. Wang, S.-U. Kum, C. Chen, and A. Kokaram, “A perceptual visibility metric for banding artifacts,” in _2016 IEEE International Conference on Image Processing (ICIP)_ , 2016, pp. 2067–2071.
  * [45] Q. Yan, L. Zhang, Y. Liu, Y. Zhu, J. Sun, Q. Shi, and Y. Zhang, “Deep hdr imaging via a non-local network,” _IEEE Transactions on Image Processing_ , vol. 29, pp. 4308–4322, 2020.
  * [46] X. Yang, K. Xu, Y. Song, Q. Zhang, X. Wei, and R. W. Lau, “Image correction via deep reciprocating hdr transformation,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2018, pp. 1798–1807.
  * [47] Q. Ye, J. Xiao, K.-m. Lam, and T. Okatani, _Progressive and Selective Fusion Network for High Dynamic Range Imaging_. New York, NY, USA: Association for Computing Machinery, 2021, p. 5290–5297. [Online]. Available: https://doi.org/10.1145/3474085.3475651
  * [48] X. Zhang, T. Kashti, D. Kella, T. Frank, D. Shaked, R. Ulichney, M. Fischer, and J. P. Allebach, “Measuring the modulation transfer function of image capture devices: what do the numbers really mean?” in _Image Quality and System Performance IX_ , vol. 8293. SPIE, 2012, pp. 64–74.

  *[HDR]: high dynamic range
  *[EXIF]: exchangeable image file format
  *[CRF]: camera response function
  *[BTF]: brightness transfer function
  *[MLE]: maximum likelihood estimation
  *[MST]: maximum spanning trees
  *[RMSE]: root-mean-squared error


The Belle Collaboration

#
Search for a heavy neutrino in $\tau$ decays at Belle

D. Liventsev  0000-0003-3416-0056 I. Adachi  0000-0003-2287-0173 H. Aihara
0000-0002-1907-5964 S. Al Said  0000-0002-4895-3869 D. M. Asner
0000-0002-1586-5790 H. Atmacan  0000-0003-2435-501X R. Ayad
0000-0003-3466-9290 V. Babu  0000-0003-0419-6912 Sw. Banerjee
0000-0001-8852-2409 M. Bauer  0000-0002-0953-7387 P. Behera
0000-0002-1527-2266 K. Belous  0000-0003-0014-2589 J. Bennett
0000-0002-5440-2668 M. Bessner  0000-0003-1776-0439 T. Bilka
0000-0003-1449-6986 D. Biswas  0000-0002-7543-3471 D. Bodrov
0000-0001-5279-4787 G. Bonvicini  0000-0003-4861-7918 J. Borah
0000-0003-2990-1913 A. Bozek  0000-0002-5915-1319 M. Bračko
0000-0002-2495-0524 P. Branchini  0000-0002-2270-9673 T. E. Browder
0000-0001-7357-9007 A. Budano  0000-0002-0856-1131 M. Campajola
0000-0003-2518-7134 D. Červenkov  0000-0002-1865-741X M.-C. Chang
0000-0002-8650-6058 A. Chen  0000-0002-8544-9274 B. G. Cheon
0000-0002-8803-4429 H. E. Cho  0000-0002-7008-3759 K. Cho  0000-0003-1705-7399
S.-J. Cho  0000-0002-1673-5664 S.-K. Choi  0000-0003-2747-8277 Y. Choi
0000-0003-3499-7948 S. Choudhury  0000-0001-9841-0216 D. Cinabro
0000-0001-7347-6585 S. Das  0000-0001-6857-966X G. De Nardo
0000-0002-2047-9675 G. De Pietro  0000-0001-8442-107X R. Dhamija
0000-0001-7052-3163 F. Di Capua  0000-0001-9076-5936 Z. Doležal
0000-0002-5662-3675 T. V. Dong  0000-0003-3043-1939 D. Dossett
0000-0002-5670-5582 D. Epifanov  0000-0001-8656-2693 T. Ferber
0000-0002-6849-0427 D. Ferlewicz  0000-0002-4374-1234 B. G. Fulsom
0000-0002-5862-9739 R. Garg  0000-0002-7406-4707 V. Gaur  0000-0002-8880-6134
A. Giri  0000-0002-8895-0128 P. Goldenzweig  0000-0001-8785-847X E. Graziani
0000-0001-8602-5652 T. Gu  0000-0002-1470-6536 K. Gudkova  0000-0002-5858-3187
C. Hadjivasiliou  0000-0002-2234-0001 X. Han  0000-0003-1656-9413 T. Hara
0000-0002-4321-0417 K. Hayasaka  0000-0002-6347-433X H. Hayashii
0000-0002-5138-5903 M. T. Hedges  0000-0001-6504-1872 D. Herrmann
0000-0001-9772-9989 M. Hernández Villanueva  0000-0002-6322-5587 C.-L. Hsu
0000-0002-1641-430X K. Inami  0000-0003-2765-7072 N. Ipsita
0000-0002-2927-3366 A. Ishikawa  0000-0002-3561-5633 R. Itoh
0000-0003-1590-0266 M. Iwasaki  0000-0002-9402-7559 W. W. Jacobs
0000-0002-9996-6336 E.-J. Jang  0000-0002-1935-9887 S. Jia
0000-0001-8176-8545 Y. Jin  0000-0002-7323-0830 A. B. Kaliyar
0000-0002-2211-619X K. H. Kang  0000-0002-6816-0751 T. Kawasaki
0000-0002-4089-5238 C. Kiesling  0000-0002-2209-535X C. H. Kim
0000-0002-5743-7698 D. Y. Kim  0000-0001-8125-9070 K.-H. Kim
0000-0002-4659-1112 Y.-K. Kim  0000-0002-9695-8103 P. Kodyš
0000-0002-8644-2349 A. Korobov  0000-0001-5959-8172 S. Korpar
0000-0003-0971-0968 E. Kovalenko  0000-0001-8084-1931 P. Križan
0000-0002-4967-7675 P. Krokovny  0000-0002-1236-4667 M. Kumar
0000-0002-6627-9708 R. Kumar  0000-0002-6277-2626 K. Kumara
0000-0003-1572-5365 Y.-J. Kwon  0000-0001-9448-5691 T. Lam
0000-0001-9128-6806 J. S. Lange  0000-0003-0234-0474 S. C. Lee
0000-0002-9835-1006 P. Lewis  0000-0002-5991-622X L. K. Li
0000-0002-7366-1307 Y. Li  0000-0002-4413-6247 L. Li Gioi  0000-0003-2024-5649
J. Libby  0000-0002-1219-3247 K. Lieret  0000-0003-2792-7511 Y.-R. Lin
0000-0003-0864-6693 T. Luo  0000-0001-5139-5784 T. Matsuda
0000-0003-4673-570X D. Matvienko  0000-0002-2698-5448 F. Meier
0000-0002-6088-0412 M. Merola  0000-0002-7082-8108 F. Metzner
0000-0002-0128-264X K. Miyabayashi  0000-0003-4352-734X R. Mizuk
0000-0002-2209-6969 G. B. Mohanty  0000-0001-6850-7666 R. Mussa
0000-0002-0294-9071 I. Nakamura  0000-0002-7640-5456 T. Nakano
0000-0003-3157-5328 M. Nakao  0000-0001-8424-7075 Z. Natkaniec
0000-0003-0486-9291 A. Natochii  0000-0002-1076-814X L. Nayak
0000-0002-7739-914X M. Nayak  0000-0002-2572-4692 N. K. Nisar
0000-0001-9562-1253 S. Nishida  0000-0001-6373-2346 S. Ogawa
0000-0002-7310-5079 H. Ono  0000-0003-4486-0064 P. Oskin  0000-0002-7524-0936
P. Pakhlov  0000-0001-7426-4824 G. Pakhlova  0000-0001-7518-3022 S. Pardi
0000-0001-7994-0537 H. Park  0000-0001-6087-2052 J. Park  0000-0001-6520-0028
A. Passeri  0000-0003-4864-3411 S. Patra  0000-0002-4114-1091 S. Paul
0000-0002-8813-0437 T. K. Pedlar  0000-0001-9839-7373 R. Pestotnik
0000-0003-1804-9470 L. E. Piilonen  0000-0001-6836-0748 T. Podobnik
0000-0002-6131-819X E. Prencipe  0000-0002-9465-2493 M. T. Prim
0000-0002-1407-7450 A. Rostomyan  0000-0003-1839-8152 N. Rout
0000-0002-4310-3638 G. Russo  0000-0001-5823-4393 S. Sandilya
0000-0002-4199-4369 A. Sangal  0000-0001-5853-349X L. Santelj
0000-0003-3904-2956 V. Savinov  0000-0002-9184-2830 G. Schnell
0000-0002-7336-3246 J. Schueler  0000-0002-2722-6953 C. Schwanda
0000-0003-4844-5028 Y. Seino  0000-0002-8378-4255 K. Senyo
0000-0002-1615-9118 M. E. Sevior  0000-0002-4824-101X M. Shapkin
0000-0002-4098-9592 C. Sharma  0000-0002-1312-0429 C. P. Shen
0000-0002-9012-4618 J.-G. Shiu  0000-0002-8478-5639 B. Shwartz
0000-0002-1456-1496 F. Simon  0000-0002-5978-0289 A. Sokolov
0000-0002-9420-0091 E. Solovieva  0000-0002-5735-4059 M. Starič
0000-0001-8751-5944 M. Sumihama  0000-0002-8954-0585 T. Sumiyoshi
0000-0002-0486-3896 W. Sutcliffe  0000-0002-9795-3582 M. Takizawa
0000-0001-8225-3973 U. Tamponi  0000-0001-6651-0706 K. Tanida
0000-0002-8255-3746 F. Tenchini  0000-0003-3469-9377 M. Uchida
0000-0003-4904-6168 S. Uehara  0000-0001-7377-5016 S. Uno  0000-0002-3401-0480
Y. Usov  0000-0003-3144-2920 R. van Tonder  0000-0002-7448-4816 G. Varner
0000-0002-0302-8151 K. E. Varvell  0000-0003-1017-1295 A. Vinokurova
0000-0003-4220-8056 M.-Z. Wang  0000-0002-0979-8341 X. L. Wang
0000-0001-5805-1255 M. Watanabe  0000-0001-6917-6694 S. Watanuki
0000-0002-5241-6628 O. Werbycka  0000-0002-0614-8773 E. Won
0000-0002-4245-7442 X. Xu  0000-0001-5096-1182 B. D. Yabsley
0000-0002-2680-0474 W. Yan  0000-0003-0713-0871 S. B. Yang
0000-0002-9543-7971 J. Yelton  0000-0001-8840-3346 J. H. Yin
0000-0002-1479-9349 C. Z. Yuan  0000-0002-1652-6686 L. Yuan
0000-0002-6719-5397 Z. P. Zhang  0000-0001-6140-2044 V. Zhilich
0000-0002-0907-5565 V. Zhukova  0000-0002-8253-641X

###### Abstract

We report on a search for a heavy Majorana neutrino in the decays
$\tau^{-}\to\pi^{-}\nu_{h}$, $\nu_{h}\to\pi^{\pm}\ell^{\mp}$, $\ell=e,\mu$.
The results are obtained using the full data sample of $988\,\mathrm{fb}^{-1}$
collected with the Belle detector at the KEKB asymmetric energy $e^{+}e^{-}$
collider, which contains $912\times 10^{6}$ $\tau\tau$ pairs. We observe no
significant signal and set 95% CL upper limits on the couplings of the heavy
right-handed neutrinos to the conventional SM left-handed neutrinos in the
mass range $0.2-1.6\,\textrm{GeV}/c^{2}$. This is the first study of a mixed
couplings of heavy neutrinos to $\tau$ leptons and light-flavor leptons.

Tau decay, Heavy neutral lepton

###### pacs:

12.60.-i,13.35.Dx,14.60.Pq

In the Standard model (SM), neutrinos are strictly massless since there are no
right-handed neutrino components. However, experimental data on neutrino
oscillations conclusively show that neutrinos are massive [1], though neutrino
mass measurements show that their masses are very small [2]. One approach to
resolve this disagreement is to include right-handed neutrinos, also known as
sterile neutrinos, heavy neutrinos, or heavy neutral leptons (HNL), into the
model. Such particles do not participate in any of the weak, strong, and
electromagnetic interactions; if we exclude gravitation, the only way they
interact with matter is via mixing with left-handed neutrinos. Singlet right-
handed neutrinos may also have Majorana mass, naturally explaining the
smallness of the observed neutrino masses via the so-called “see-saw”
mechanism [3]. One example of the models realizing such a mechanism is
$\nu$MSM [4]. It introduces three right-handed singlet HNLs, so that every
left-handed particle gets its right-handed counterpart, and manages to explain
neutrino oscillations, dark matter existence, and baryogenesis with the same
set of parameters. HNLs also appear in other extensions of the SM; see a
review in Ref. [5].

In general, neutrino flavor eigenstates need not to coincide with the mass
eigenstates but may be related through a unitary transformation, similar to
that in the quark sector:

$\nu_{\alpha}=\sum_{i}{U_{\alpha
i}\nu_{i}},\quad\alpha=e,\mu,\tau,...,\,i=1,2,3,4,...$ (1)

where Greek (Latin) indexes denote flavor (mass) eigenstates. The coupling of
the HNLs to charged or neutral currents of flavor $\alpha$ is characterized by
the quantities $U_{\alpha 4}$, $U_{\alpha 5}$ etc, which we denote for
convenience as $U_{\alpha}$. A generic HNL is denoted by $\nu_{h}$. Its
production and decay diagrams are shown in Fig. 1. Existing experimental
results are reviewed and discussed in Ref. [5].

In our previous analysis [6] we searched for the decays of HNLs produced in
$B$ decays. No signal was found and upper limits on $|U_{e}|^{2}$,
$|U_{\mu}|^{2}$ and $|U_{e}||U_{\mu}|$ as functions of the mass of the HNL
were set.

In this analysis, we reconstruct $\tau^{-}\to\pi^{-}_{p}\nu_{h}$ decays, where
$\pi_{p}$ refers to the “prompt pion” and the HNL decays into a pion-lepton
pair in the detector volume: $\nu_{h}\to\pi^{\pm}\ell^{\mp}$, $\ell=e,\mu$
(the charge-conjugate decay mode being included throughout this Letter). Both
combinations of $\tau$ and $\ell$ charges are retained for further analysis.
In the final state, we have two pions and a lepton:
$\pi^{-}_{p}\pi^{\pm}\ell^{\mp}$.

Following [7, 8], we interpret the result in terms of the minimal realistic
model with two quasi-degenerate HNLs with close masses and couplings and not
trivial $U_{\alpha}$. When two HNL masses are not exactly the same, HNL
oscillations occur and we consider two extreme cases: the “Dirac-like limit”,
where only lepton-number conserving final states are allowed, and the
“Majorana-like limit”, where lepton-number violating final states are also
allowed with the same branching fractions. The ratio of different $U_{\alpha}$
is determined from the neutrino oscillation data. In the normal hierarchy case
(NH), the relative mixing coefficients
$x_{\alpha}\equiv|U_{\alpha}|^{2}/|U|^{2}$,
$|U|^{2}=\sum_{\alpha}|U_{\alpha}|^{2}$ $(\alpha=e,\mu,\tau)$ are taken to be
$x_{e}=0.06$, $x_{\mu}=0.48$ and $x_{\tau}=0.46$; for the inverted hierarchy
case (IH), we use the values $x_{\alpha}=1/3$ $(\alpha=e,\mu,\tau)$ [7].

A distinctive feature of the HNL is its long lifetime. We can estimate it as
$c\tau\sim|U|^{-2}M(\nu_{h})^{-5}$ [9]: for $M(\nu_{h})=1\,\textrm{GeV}/c^{2}$
and $|U|^{2}=10^{-4}$, the lifetime is $c\tau\sim 20\,\textrm{m}$; thus, a
$\pi\ell$ pair forms a vertex displaced from the interaction point (IP). BaBar
[10] and Belle [11] previously searched for $\tau\to\ell hh^{\prime}$ decays;
however, both analyses required all tracks to originate in the vicinity of the
IP. For the long-lived HNL, this greatly reduces the reconstruction
efficiency. In contrast, we do not impose such a requirement on the HNL
daughters.

Results presented here are based on all available Belle data, including on-
resonance, off-resonance and energy scans. The collision energy is around
10.58 GeV. The total integrated luminosity is $988\,\mathrm{fb}^{-1}$ [12] and
the total number of $\tau^{+}\tau^{-}$ pairs is calculated using direct
production cross sections [13] and $\Upsilon(NS)$ branching fractions to be
$N_{\tau\tau}=(912\pm 13)\times 10^{6}$, where the error arises from the
luminosity measurement.

The Belle detector is a large-solid-angle magnetic spectrometer that consists
of a silicon vertex detector (SVD), a 50-layer central drift chamber (CDC), an
array of aerogel threshold Cherenkov counters (ACC), a barrel-like arrangement
of time-of-flight scintillation counters (TOF), and an electromagnetic
calorimeter comprised of CsI(Tl) crystals (ECL) located inside a super-
conducting solenoid coil that provides a 1.5 T magnetic field. An iron flux-
return located outside of the coil is instrumented to detect $K_{L}^{0}$
mesons and to identify muons (KLM). The detector is described in detail
elsewhere [14].

To study backgrounds, we use the following Monte Carlo (MC) simulated samples:
$e^{+}e^{-}\to q\bar{q}$ $(q=u,d,s,c,b)$, $e^{+}e^{-}$, $\mu^{+}\mu^{-}$,
$\tau^{+}\tau^{-}$ and two-photon processes ($e^{+}e^{-}\to
e^{+}e^{-}\ell^{+}\ell^{-}$, $e^{+}e^{-}q\bar{q}$). Background processes are
generated by EvtGen [15], BHLUMI [18], KKMC [16], KKMC ($\tau\tau$ production)
and TAUOLA ($\tau$ decay) [17], and AAFH [19], respectively. We use signal MC
samples generated with different HNL masses of
$M(\nu_{h})=0.2\,\textrm{GeV}/c^{2}$ to $1.6\,\textrm{GeV}/c^{2}$ (with
$0.1\,\textrm{GeV}/c^{2}$ step) and life times of $c\tau=0.2,0.5,$ and
$1.0\,\textrm{m}$ to study the response of the detector and to determine its
acceptance and signal efficiency dependence on the neutrino mass and the
distance of the decay point from the IP. This efficiency does not depend on
$c\tau$ ($|U|^{2}$). Signal MC samples are $\tau^{+}\tau^{-}$ pairs where one
of the $\tau$ leptons decays according to the modes under study and the other
decays generically. Signal events are generated using EvtGen; radiative
corrections are included using PHOTOS [20]. HNLs are produced and decayed
uniformly in phase space. GEANT3 [21] is used to model the detector response.

Electrons are identified using the energy and shower profile in the ECL, the
light yield in the ACC and the specific ionization energy loss in the CDC
($dE/dx$). This information is used to form an electron ($\mathcal{L}_{e}$)
and non-electron ($\mathcal{L}_{\bar{e}}$) likelihood; these are combined into
a likelihood ratio
$\mathcal{P}_{e}=\mathcal{L}_{e}/(\mathcal{L}_{e}+\mathcal{L}_{\bar{e}})$
[22]. Muons are distinguished from other charged tracks by their range and hit
profiles in the KLM. This information is utilized in a likelihood ratio
approach [23] similar to the one used for the electron identification.

Charged tracks with laboratory momentum greater than $0.5\,\textrm{GeV}/c$ and
electron likelihood ratio $\mathcal{P}_{e}>0.9$ or muon ratio
$\mathcal{P}_{\mu}>0.9$ are treated as leptons. These requirements correctly
identify leptons with an efficiency of approximately 95% and a
misidentification rate of less than 2%. All charged tracks not identified as
leptons and satisfying the electron veto ${\cal P}_{e}<0.5$ are treated as
pions.

HNL candidates are formed from a pion $\pi$ and a lepton $\ell$ of the
opposite sign. The pion and lepton are then fitted to a common vertex. HNL
candidates are combined with a prompt pion $\pi_{p}$. The second vertex fit of
the HNL candidate and the prompt pion is performed with a vertex constraint to
the IP, which is determined run-by-run using charged tracks. The $\chi^{2}$ of
both vertex fits is required to be $\chi^{2}/ndf<25$, where $ndf$ is the
number of degrees of freedom. Kinematics of the particles are updated after
the fits are performed. For the prompt pion, we require the closest distance
to IP along the detector symmetry axis ($dz$) to be $|dz|<5\,\textrm{cm}$ and
in the transverse plane to be $dr<1\,\textrm{cm}$.

Since the signal $\tau$ lepton is fully reconstructed, we can utilize the
kinematic constraint of the known initial four-momentum of the colliding
$e^{+}e^{-}$ pair to define $\Delta E\equiv E(\pi_{p}\pi\ell)-E_{\textrm{cm}}$
in the center-of-mass (c.m.) frame of a $\tau$ candidate relative to the beam
energy $E_{\textrm{cm}}$. In $\tau$ decays, $\Delta E$ and $M(\pi_{p}\pi\ell)$
are highly correlated; therefore, we use an elliptically shaped requirement,
that encompasses $\sim$ 95% efficiency as computed from the signal simulation.

In the rest of the event, we select tracks with $dr<1\,\textrm{cm}$,
$|dz|<5\,\textrm{cm}$ and a transverse momentum $p_{t}>0.1\,\textrm{GeV}/c$.
We classify clusters in the ECL not associated with charged tracks as photons
and require $E(\gamma)>0.05\,\textrm{GeV}$ in the barrel
($32.2^{\circ}<\theta<128.7^{\circ}$), $E(\gamma)>0.1\,\textrm{GeV}$ in the
forward endcap ($12.4^{\circ}<\theta<31.4^{\circ}$) and
$E(\gamma)>0.15\,\textrm{GeV}$ in the backward endcap
($130.7^{\circ}<\theta<155.1^{\circ}$). Events are separated into two
hemispheres by the plane perpendicular to the thrust axis $\vec{n}_{T}$ [24],
defined to maximize the thrust magnitude value

$V_{T}=\frac{\sum|\vec{p}_{i}^{*}\cdot\vec{n}_{T}|}{\sum|\vec{p}_{i}^{*}|},$
(2)

where $\vec{p}_{i}^{*}$ are momenta of the selected tracks, photons and $\tau$
daughters in the c.m. frame. We require the signal hemisphere to contain no
additional tracks besides $\tau$ daughters, and the opposite side to contain
one or three tracks with a total charge opposite that of the prompt pion.

We select well-vertexed HNL candidates using $d\phi$, the angle between the
momentum vector and decay-vertex vector of the HNL candidate;
$dz_{\textrm{vtx}}$, the distance between the daughter tracks at their closest
approach in the direction parallel to the beam; and $dr$ for each track.
Requirements vary depending on the presence of SVD hits on the tracks and on
the HNL candidate flight length. These are summarized in Table 1. The four
event types in the Table are I: both neutrino daughter tracks have recorded
hits in SVD, II: one of the neutrino daughter tracks has recorded hits in SVD,
III: none of the neutrino daughter tracks have recorded hits in SVD, with
$M(\nu_{h})\geq 0.8\,\textrm{GeV}/c^{2}$, and IV: no SVD hits and
$M(\nu_{h})<0.8\,\textrm{GeV}/c^{2}$. There is a large contamination from the
conversion photons $\gamma\to ee$ in the last category.

Table 1: Summary of the vertex requirements. Type |  $d\phi$, rad |  $z_{\mathrm{vtx}}$, cm |  $dr$, cm
---|---|---|---
I | $<0.02$ | $<0.06$ | $>0.07$
II | $<0.024$ | $<1.5$ | $>0.08$
III | $<0.16$ | $<3.0$ | $>0.1$
IV | $<0.16$ | $<3.0$ | $>1.0$

Figure 2 shows the efficiency of HNL reconstruction with all requirements
applied as a function of the reconstructed travel distance $l$ for several
mass hypotheses.

The number of neutrinos detected by this method is (in units where
$\hbar=c=1$) [25]

$\displaystyle n(\nu_{h})$ $\displaystyle=2N_{\tau\tau}\
\mathcal{B}(\tau\to\pi\nu_{h})\ \mathcal{B}(\nu_{h}\to\pi\ell)\
\frac{m\Gamma}{p}\int\exp\Big{(}-\frac{m\Gamma l}{p}\Big{)}\varepsilon(m,l)dl$
$\displaystyle=|U_{\tau}|^{2}|U_{\ell}|^{2}\ 2N_{\tau\tau}\ f_{1}(m)\
f_{2}(m)\ \frac{m}{p}\int\exp\Big{(}-\frac{m\Gamma
l}{p}\Big{)}\varepsilon(m,l)dl,$ (3)

where $N_{\tau\tau}$ is the number of $\tau\tau$ pairs,
$\mathcal{B}(\tau\to\pi\nu_{h})$ is the branching fraction for $\nu_{h}$
production, $\mathcal{B}(\nu_{h}\to\pi\ell)$ is the branching fraction of the
reconstructed decay, $m$, $p$ and $\Gamma=\Gamma(m,U)$ are the mass, momentum
and full width of the HNL, respectively, and $\varepsilon(m,l)$ is the
reconstruction efficiency of the HNL of mass $m$ decaying at a distance $l$
from the IP. The momentum $p$ is approximated by the mean value for a given
mass, determined from the signal MC simulation. To factor out the
$|U_{\ell}|^{2}$ dependence, we define functions $f_{1,2}(m)$ as
$|U_{\tau}|^{2}f_{1}(m)\equiv\mathcal{B}(\tau\to\pi\nu_{h})$ and
$|U_{\ell}|^{2}f_{2}(m)\equiv\Gamma(\nu_{h}\to\pi\ell)=\mathcal{B}(\nu_{h}\to\pi\ell)\Gamma$,
where $\ell$ denotes the flavor ($e,\ \mu$) of the charged lepton produced in
the $\nu_{h}$ decay. Integration is performed over the full volume used to
reconstruct the HNL vertex. The expressions for
$\mathcal{B}(\tau\to\pi\nu_{h})$, $\Gamma(\nu_{h}\to\pi\ell)$ and the full
neutrino width $\Gamma$ are taken from Ref. [26] and require only general
assumptions (i.e., they are not specific to the $\nu$MSM model). In the
Majorana case $\Gamma(\nu_{h}\to\pi\ell)$ is twice that in the Dirac case.
Given number of observed events, we solve Eq.  Search for a heavy neutrino in
$\tau$ decays at Belle for the variable $|U|^{2}$ using the relative mixing
coefficients $x_{\alpha}$ defined above.

The systematic uncertainties in number of events calculated according to Eq.
Search for a heavy neutrino in $\tau$ decays at Belle are enumerated in Table
2. We estimate the systematic uncertainties of event selection criteria from
the differences in their efficiencies obtained in data and MC simulation.
Since all particles used in the systematic uncertainty study decay relatively
close to the IP compared to the expectation for an HNL, we require where
possible that the decay vertices be farther than 4 cm from the IP in the
transverse plane to put more weight on large decay lengths. To estimate the
systematic uncertainty due to tracking, we compare the number of fully and
partially reconstructed $D^{*+}$ decays in the decay chain $D^{*+}\to
D^{0}\pi^{+}$, $D^{0}\to K_{S}^{0}\pi^{+}\pi^{-}$,
$K_{S}^{0}\to\pi^{+}\pi^{-}$, where in the latter case one of the pions from
the $K_{S}^{0}$ is explicitly left unreconstructed. To estimate the systematic
uncertainty of the lepton identification, we reconstruct
$J/\psi\to\ell^{+}\ell^{-}$, $\ell=e,\mu$ events, where one of the daughter
particles is identified as an electron or muon. The difference of the
identification efficiency of the other daughter between data and MC simulation
is treated as a systematic uncertainty. To estimate the systematic uncertainty
of the vertex quality requirements we apply them to $K^{0}_{S}$ decays, which
have a topology similar to HNL decays. Signal events were generated using
EvtGen, which is not optimized for $\tau$ decays. To estimate the effect of
this we prepared two samples of $\tau\tau$ events — one generated with EvtGen
and the other with KKMC and TAUOLA — then reconstructed $\tau\to\ell\nu\nu$
decays, where $\ell=e,\mu$, using the same tagging as for the signal events
and compared reconstruction efficiency in both cases. The phase space model
may not give the correct angle distribution of the HNL. We vary it by
reweighting generated events and treat the change as a systematic uncertainty.
The calculation uncertainty comes from the efficiency and momentum
approximations in Eq.  Search for a heavy neutrino in $\tau$ decays at Belle
and was estimated by comparing predicted and observed number of events in
different subsets of the signal MC simulation. Systematic uncertainties
induced by the fitting procedure were found by varying the signal resolution
and background shape within their errors. The theoretical uncertainty arises
from uncertainties in the constants used in Eq.  Search for a heavy neutrino
in $\tau$ decays at Belle. Correlations between different systematic
uncertainties are found to be small and are ignored. All systematic
uncertainties are summed in quadrature, leading to total systematic
uncertainties of 16% and 12% for the $\pi\pi e$ and $\pi\pi\mu$ modes,
respectively. The largest contributions are lepton identification (12% and 6%
for the electron and muon identification, respectively) and vertex quality
requirements (5.3%).

Table 2: Summary of systematic uncertainties in number of events calculated according to Eq.  Search for a heavy neutrino in $\tau$ decays at Belle. Requirement | Systematic uncertainty, %
---|---
Tracking | 1.2
${\cal P}_{e}(\ell)$ | 12
${\cal P}_{\mu}(\ell)$ | 6
Vertex quality | 5.3
$N_{\tau\tau}$ | 1.4
Generator | 2
Angle distribution | 5
Calculation | 4
Signal resolution | 2.5
Background shape | 4
Theoretical | 0.35
Total, $e/\mu$ | 16/12

We study the dependence of the coupling constant $|U|^{2}$ on the HNL mass
using simultaneous fit to $\pi\pi e$ and $\pi\pi\mu$ modes taking into account
the relative mixing coefficients $x_{\alpha}$ defined above.

The $\Delta E$ vs $M(\pi\pi\ell)$ distributions with all requirements but
$\Delta E$ and $M(\pi\pi\ell)$ imposed are shown in Fig. 3. The mass
distributions after application of all reconstruction requirements are shown
in Fig. 4 for the same-charge $\tau$ and $\ell$ (“Dirac-like limit”) and both
same- and opposite-charge combinations (“Majorana-like limit”). From the
background MC simulation study, we expect to see a wide peak around $\sim
0.2\,\textrm{GeV}/c^{2}$ from the conversion process $\gamma\to e^{+}e^{-}$ in
the $\pi\pi e$ mode and a narrow peak from the $K_{S}^{0}\to\pi^{+}\pi^{-}$
process at $\sim 0.48\,\textrm{GeV}/c^{2}$ in the $\pi\pi\mu$ mode. Since the
conversion peak is wide, we can distinguish a narrow signal peak from the HNL
decay under it, but since the $K_{S}^{0}$ peak is narrow, we exclude the
$K_{S}^{0}$ region at $0.464-0.494\,\textrm{GeV}/c^{2}$ from consideration,
which corresponds to $\pm 2\sigma$ of the peak width.

The HNL mass is unknown and we search for it in the kinematically accessible
region for the mass; for the decays under study, this lies between
$M_{\pi}+M_{\ell}$ and $M_{\tau}-M_{\pi}$. We perform a series of binned
likelihood fits to the mass distributions using the sum of a Gaussian signal
function and background (described below) varying the mass hypothesis in each
fit. The neutrino mass is set to the center of a histogram bin which has a
width of $2\,\textrm{MeV}/c^{2}$. The signal-shape parameters used in the fits
to data are fixed to those obtained by fitting simulated events. The width
evolves linearly from $\sim 3\,\textrm{MeV}/c^{2}$ for
$M(\nu_{h})=0.2\,\textrm{GeV}/c^{2}$ to $\sim 10\,\textrm{MeV}/c^{2}$ for
$M(\nu_{h})=1.6\,\textrm{GeV}/c^{2}$. The background is described by the sum
of a constant and the conversion peak described above in the $\pi\pi e$ subset
and by a constant and the $K_{S}^{0}$ peak in the $\pi\pi\mu$ subset. The
functions for the peaking components are defined as smoothed histograms from
the background MC simulation. Yields of all components are free parameters of
the fit.

The statistical significance of the HNL signal is defined as
$\sqrt{-2\ln{L_{0}}/{L}}$, where $L_{0}$ and $L$ are the likelihoods returned
by the fit with the signal yield fixed at zero and at the fitted value,
respectively. The maximum local significance in four fits does not exceed
$2.5\,\sigma$, and we set upper limits on $|U|^{2}$ at 95% CL in the “Dirac-
like limit” and the “Majorana-like limit” for the two neutrino-mass hierarchy
scenarios bin by bin using pyhf package [27, 28]. The resulting upper limits
on the coupling constants at 95% CL are shown in Fig. 5. Comparison of all
four upper limits on one plot may be found in the Supplimentary Materials
[25].

In conclusion, we search for a heavy neutrino in $\tau$ decays and observe no
significant signal. This is the first study of a mixed couplings of heavy
neutrinos to $\tau$ leptons and light-flavor leptons. Upper limits on the
mixing of HNLs to the active neutrinos in the mass range
$0.2-1.6\,\textrm{GeV}/c^{2}$ are set. The maximum sensitivities are achieved
around $1.0\,\textrm{GeV}/c^{2}$ and the corresponding upper limits at 95% CL
for $|U|^{2}$ are $1.4\times 10^{-4}$ ($1.5\times 10^{-4}$) in the “Dirac-like
limit” for the normal (inverted) hierarchy and $1.0\times 10^{-4}$ ($1.1\times
10^{-4}$) in the “Majorana-like limit” for the normal (inverted) hierarchy.

This work, based on data collected using the Belle detector, which was
operated until June 2010, was supported by the Ministry of Education, Culture,
Sports, Science, and Technology (MEXT) of Japan, the Japan Society for the
Promotion of Science (JSPS), and the Tau-Lepton Physics Research Center of
Nagoya University; the Australian Research Council including grants
DP180102629, DP170102389, DP170102204, DE220100462, DP150103061, FT130100303;
Austrian Federal Ministry of Education, Science and Research (FWF) and FWF
Austrian Science Fund No. P 31361-N36; the National Natural Science Foundation
of China under Contracts No. 11675166, No. 11705209; No. 11975076; No.
12135005; No. 12175041; No. 12161141008; Key Research Program of Frontier
Sciences, Chinese Academy of Sciences (CAS), Grant No. QYZDJ-SSW-SLH011; the
Ministry of Education, Youth and Sports of the Czech Republic under Contract
No. LTT17020; the Czech Science Foundation Grant No. 22-18469S; Horizon 2020
ERC Advanced Grant No. 884719 and ERC Starting Grant No. 947006 “InterLeptons”
(European Union); the Carl Zeiss Foundation, the Deutsche
Forschungsgemeinschaft, the Excellence Cluster Universe, and the
VolkswagenStiftung; the Department of Atomic Energy (Project Identification
No. RTI 4002) and the Department of Science and Technology of India; the
Istituto Nazionale di Fisica Nucleare of Italy; National Research Foundation
(NRF) of Korea Grant Nos. 2016R1D1A1B02012900, 2018R1A2B3003643,
2018R1A6A1A06024970, RS202200197659, 2019R1I1A3A01058933, 2021R1A6A1A03043957,
2021R1F1A1060423, 2021R1F1A1064008, 2022R1A2C1003993; Radiation Science
Research Institute, Foreign Large-size Research Facility Application
Supporting project, the Global Science Experimental Data Hub Center of the
Korea Institute of Science and Technology Information and KREONET/GLORIAD; the
Polish Ministry of Science and Higher Education and the National Science
Center; the Ministry of Science and Higher Education of the Russian
Federation, Agreement 14.W03.31.0026, and the HSE University Basic Research
Program, Moscow; University of Tabuk research grants S-1440-0321, S-0256-1438,
and S-0280-1439 (Saudi Arabia); the Slovenian Research Agency Grant Nos.
J1-9124 and P1-0135; Ikerbasque, Basque Foundation for Science, Spain; the
Swiss National Science Foundation; the Ministry of Education and the Ministry
of Science and Technology of Taiwan; and the United States Department of
Energy and the National Science Foundation. These acknowledgements are not to
be interpreted as an endorsement of any statement made by any of our
institutes, funding agencies, governments, or their representatives. We thank
the KEKB group for the excellent operation of the accelerator; the KEK
cryogenics group for the efficient operation of the solenoid; and the KEK
computer group and the Pacific Northwest National Laboratory (PNNL)
Environmental Molecular Sciences Laboratory (EMSL) computing group for strong
computing support; and the National Institute of Informatics, and Science
Information NETwork 6 (SINET6) for valuable network support.

## References

  * [1] R. L. Workman et al. (Particle Data Group), Prog. Theor. Exp. Phys. 2022, Issue 8, 083C01 (2022).
  * [2] M. Aker et al. (KATRIN Collaboration), Nat. Phys. 18, 160 (2022) [arXiv:2105.08533].
  * [3] P. Minkowski, Phys. Lett. B 67, 421 (1977);
T. Yanagida, in Proc. of the Workshop on Grand Unified Theory and Baryon
Number of the Universe, KEK, Japan, 1979;
M. Gell-Mann, P. Ramond and R. Slansky, (1980), Print-80-0576 (CERN);
R.N. Mohapatra and G. Senjanovic, Phys. Rev. Lett. 44, 912 (1980).

  * [4] T. Asaka, S. Blanchet and M. Shaposhnikov, Phys. Lett. B 631, 151 (2005) [arXiv:hep-ph/0503065].
  * [5] F. F. Deppisch, P. S. Bhupal Dev, A. Pilaftsis, New J. Phys. 17, 075019 (2015) [arXiv:1502.06541];
P. D. Bolton, F. F. Deppisch, P. S. Bhupal Dev, J. High Energ. Phys. 03, 170
(2020) [arXiv:1912.03058].

  * [6] D. Liventsev et al. (Belle Collaboration), Phys. Rev. D. 87, 071102 (2013) [arXiv:1301.1105].
  * [7] J.-L. Tastet, O. Ruchayskiy, and I. Timiryasov, J. High Energ. Phys. 12 182 (2021), [arXiv:2107.12980].
  * [8] ATLAS collaboration, arXiv:2204.11988.
  * [9] M. Gronau, C. N. Leung and J. L. Rosner, Phys. Rev. D 29, 2539 (1984).
  * [10] B. Aubert et al. (BaBaR Collaboration), Phys. Rev. Lett. 95, 191801 (2005) [arXiv:hep-ex/0506066].
  * [11] Y. Miyazaki et al. (Belle Collaboration), Phys. Lett. B 719, 346 (2013) [arXiv:1206.5595].
  * [12] J. Brodzicka et al. (Belle Collaboration), Prog. Theor. Exp. Phys. 2012, Issue 1, 4D001 (2012) [arXiv:1212.5342].
  * [13] S. Banerjee, B. Pietrzyk, J. M. Roney, and Z. Was, Phys. Rev. D 77, 054012 (2008).
  * [14] A. Abashian et al. (Belle Collab.), Nucl. Instr. and Meth. A 479, 117 (2002).
  * [15] D. J. Lange, Nucl. Instr. and Meth. A 462, 152 (2001).
  * [16] S. Jadach, B. Ward, and Z. Wa̧s, Comput. Phys. Commun. 130, 260 (2000).
  * [17] N. Davidson et al., Comput. Phys. Commun. 183, 821 (2012).
  * [18] S. Jadach et al., Comput. Phys. Commun. 70, 305 (1992).
  * [19] F. A. Berends et al., Comput. Phys. Commun. 40, 285 (1986).
  * [20] E. Barberio and Z. Was, Comput. Phys. Commun. 79, 291 (1994).
  * [21] R. Brun et al., CERN Report No. DD/EE/84-1 (1984).
  * [22] K. Hanagaki et al., Nucl. Instr. and Meth. A 485, 490 (2002).
  * [23] A. Abashian et al., Nucl. Instr. and Meth. A 491, 69 (2002).
  * [24] S. Brandt, C. Peyrou, R. Sosnowski, and A. Wroblewski, Phys. Lett. 12, 57 (1964);
E. Farhi, Phys. Rev. Lett. 39, 1587 (1977).

  * [25] See Supplimentary material (Appendix A).
  * [26] K. Bondarenko, A. Boyarsky, D. Gorbunov, and O. Ruchayskiy, J. High Ener. Phys. 11, (2018) 032, [arXiv:1805.08567].
  * [27] L. Heinrich, M. Feickert, and G. Stark, “pyhf:v0.7.2”, https://github.com/scikit-hep/pyhf/releases/tag/v0.7.2, doi: 10.5281/zenodo.1169739.
  * [28] L. Heinrich et al., J. of Open Source Software 6(58), 2823 (2021), https://doi.org/10.21105/joss.02823.

---
Figure 1: Feynman diagrams for HNL production (top) and decay (bottom). Figure 2: Dependence of the HNL reconstruction efficiency on the neutrino travel distance $l$ for different neutrino masses $M(\nu_{h})$. Efficiency is almost identical for $e$ and $\mu$. |
---|---
a) | b)
Figure 3: $\Delta E$ vs $M(\pi\pi\ell)$ distributions with all requirements
but $\Delta E$ and $M(\pi\pi\ell)$ imposed for $\pi\pi e$ (a) and $\pi\pi\mu$
(b) in data. The signal region is shown as a red ellipse. Figure 4: Final
distributions of $M(\nu_{h})$ for $\pi\pi e$ and $\pi\pi\mu$ reconstruction
modes in data. The filled histograms are for candidates with opposite-charge
$\tau$ and $\ell$, while the open histograms are for candidates with same-
charge combinations. The curves are the fits with the signal yield fixed at
zero. Figure 5: Upper limits at 95% CL on $|U|^{2}$. The upper (lower) plot is
for the normal (inverted) hierarchy. The solid (dashed) curve shows the result
in the “Majorana-like limit” (the “Dirac-like limit”). The excluded area is
above the curves.

##  Appendix A Derivation of Eq. 3

If a particle with a mass $m$ and width $\Gamma$ has a momentum $p$, then the
probability that it travels distance $l$ or greater is (in units where
$\hbar=c=1$)

$P(l)=\exp\Big{(}-\frac{m\Gamma l}{p}\Big{)},$ (4)

thus probability to decay in a segment $dl$ at distance $l$ is

$dP(l)=\frac{m\Gamma}{p}\exp\Big{(}-\frac{m\Gamma l}{p}\Big{)}dl.$ (5)

The number of neutrinos detected in the Belle detector is

$\displaystyle n(\nu_{h})$ $\displaystyle=N_{0}\int\varepsilon(m,l)dP(l)$
$\displaystyle=2N_{\tau\tau}\ \mathcal{B}(\tau\to\pi\nu_{h})\
\mathcal{B}(\nu_{h}\to\pi\ell)\ \frac{m\Gamma}{p}\int\exp\Big{(}-\frac{m\Gamma
l}{p}\Big{)}\varepsilon(m,l)dl$ $\displaystyle=|U_{\tau}|^{2}|U_{\ell}|^{2}\
2N_{\tau\tau}\ f_{1}(m)\ f_{2}(m)\ \frac{m}{p}\int\exp\Big{(}-\frac{m\Gamma
l}{p}\Big{)}\varepsilon(m,l)dl$ $\displaystyle=|U|^{2}\ x_{\tau}x_{\ell}\
2N_{\tau\tau}\ f_{1}(m)\ f_{2}(m)\ \frac{m}{p}\int\exp\Big{(}-\frac{m\Gamma
l}{p}\Big{)}\varepsilon(m,l)dl,$ (6)

where $N_{0}=2N_{\tau\tau}\ \mathcal{B}(\tau\to\pi\nu_{h})\
\mathcal{B}(\nu_{h}\to\pi\ell)$ is the total number of $\pi\ell$ pairs
produced in the $\tau\to\pi\nu_{h}$, $\nu_{h}\to\pi\ell$ decay chain and the
rest of definitions is the same as in the paper.

The $|U_{\tau}|^{2}$ coupling comes from the branching fraction
$\mathcal{B}(\tau\to\pi\nu_{h})$, and the $|U_{\ell}|^{2}$ coupling comes from
the partial width
$\Gamma(\nu_{h}\to\pi\ell)=\mathcal{B}(\nu_{h}\to\pi\ell)\Gamma$,
$|U_{\alpha}|^{2}=\ x_{\alpha}|U|^{2}$. The full width $\Gamma$ is calculated
by summing over a number of two- and three-body decays taking into account the
relative mixing coefficients $x_{\alpha}$.

Figure 6: Upper limits at 95% CL on $|U|^{2}$. The solid, dashed, dotted and
dash-dotted curves show the result in the “Majorana-like limit” for normal
($x_{e}=0.06$, $x_{\mu}=0.48$ and $x_{\tau}=0.46$) and inverted
($x_{\alpha}=1/3,\alpha=e,\mu,\tau$) hierarchy and in the “Dirac-like limit”
for normal and inverted hierarchy, respectively. The excluded area is above
the curves.