ISAC

    Integrated Sensing and Communication
MIMO

    Multiple-input multiple-output
SNR

    signal-to-noise ratio
mmWave

    millimeter wave
OFDM

    orthogonal frequency division multiplexing
LoS

    line-of-sight
ULA

    uniform linear array
RF

    radio frequency
AoA

    angle of arrival
AoD

    angle of departure
HDA

    hybrid digital analog
BS

    Base Station
UE

    User Equipment
KF

    Kalman Filter

# Simultaneous Communication and Tracking in Arbitrary Trajectories via Beam-
Space Processing

Fernando Pedraza1, Saeid K. Dehkordi1, Mari Kobayashi2, Giuseppe Caire1
1Technical University of Berlin, Germany
2Technical University of Munich, Germany
Emails: {f.pedrazanieto, s.khalilidehkordi<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In this paper, we develop a beam tracking scheme for an orthogonal frequency
division multiplexing (OFDM) Integrated Sensing and Communication (ISAC)
system with a hybrid digital analog (HDA) architecture operating in the
millimeter wave (mmWave) band. Our tracking method consists of an estimation
step inspired by radar signal processing techniques, and a prediction step
based on simple kinematic equations. The hybrid architecture exploits the
predicted state information to focus only on the directions of interest,
trading off beamforming gain, hardware complexity and multistream processing
capabilities. Our extensive simulations in arbitrary trajectories show that
the proposed method can outperform state of the art beam tracking methods in
terms of prediction accuracy and consequently achievable communication rate,
and is fully capable of dealing with highly non-linear dynamic motion
patterns.

###### Index Terms:

Integrated Sensing and Communication, Beam Tracking, OFDM, Hybrid Digital
Analog.

##  I Introduction

ISAC has emerged as a key enabler for 5G and beyond wireless systems to deal
with the challenging requirement in terms of spectral efficiency,
localization, as well as power consumption [1]. In particular, the waveform
design for ISAC, sharing common hardware and spectrum, has been extensively
studied in the literature (see e.g. [2] and references therein). In this
context, OFDM has been widely investigated as an ISAC waveform, in virtue of
its availability in wireless communication systems and the capability to
achieve accurate radar estimation performance [3, 4].

Most of the envisioned ISAC applications are expected to operate in mmWave
[4]. A key feature of the mmWave channel that motivates this design is its
sparsity in the beamspace domain [5], which connects the scattering conditions
with the geometry of the environment. Due to this property, accurate channel
state information is needed in order to efficiently operate in the mmWave
channel. Since acquiring this information leads to overhead, methods for
dynamically inferring the channel state given past observations are currently
an active area of research, often under the name of beam tracking [6, 7, 8, 9,
10, 11]. Another key issue in mmWave are the large bandwidths involved,
resulting in the need for very high sampling rates and a subsequent raise in
receiver complexity and energy consumption. This problem is exacerbated when
large antenna arrays are considered and the received signal is sampled and
demodulated at each antenna output.

Most of the previous approaches to beam tracking [7, 8, 9, 10] consider
channel sounding strategies, where pilot signals are necessary to track the
channel, therefore experiencing overhead. ISAC systems eliminate this problem
by considering a radar receiver co-located with the transmitter, such that the
transmitted waveform can carry information while being perfectly known upon
reception. However, many of the current works in this framework suffer from
some practical limitations. These include (i) the use of fully digital
receiver architectures incurring high energy consumption [6, 12, 11], (ii) the
restriction to very simple trajectories such as straight roads [6, 13, 14],
and (iii) the assumption of perfect matched filtering, thus neglecting the
effect of discrete sampling in time and frequency [6, 11].

In this paper, we extend our previous work on beam refinement with HDA
architectures [15] to a tracking scenario by combining it with recent ideas
from the literature. Our method outperforms other existing approaches by
showing that well designed low dimensional observations are sufficient to
achieve almost optimal tracking performance, therefore alleviating the
hardware complexities in mmWave. Furthermore, we evaluate the proposed method
in highly non-linear environments without the commonly used assumption of
perfect matched filtering. Our simulation results demonstrate that it yields
excellent performance when accurate mobility models are not available.

##  II System Model

We consider an ISAC scenario where a Base Station (BS) with $N_{\rm tx}$
antennas and $N_{\rm tx}^{\rm rf}$ RF chains, co-located with a radar receiver
equipped with $N_{\rm rx}$ antennas and $N_{\rm rx}^{\rm rf}$ RF chains, tries
to communicate with $K$ users, while simultaneously tracking their positions
over time. We let the system operate in the mmWave band. Our model assumes
that the BS has coarse knowledge of the angular locations of the $K$ users at
each measurement epoch. Such knowledge can be acquired by a beam alignment
method (see e.g. [16]) if a new user enters the coverage region of the BS, or
predicted from previous observations for tracked users, as explained in
Section III. This section focuses on a single measurement step, therefore we
drop the measurement index to simplify notation.

###  II-A Channel model

Assuming the $K$ users are well separated in space, we adopt the widely used
mmWave radar channel model (see e.g. [4]) where where the received echoes
arrive from the line-of-sight (LoS) direction corresponding to each user, i.e.

$\displaystyle{\bf H}(t,\tau)=\sum_{k=1}^{K}h_{k}{\bf b}(\phi_{k}){\bf
a}^{\scriptscriptstyle\mathsf{H}}(\phi_{k})\delta(\tau-\tau_{k})e^{j2\pi\nu_{k}t}\,,$
(1)

where $h_{k}$, $\tau_{k}$, $\nu_{k}$ and $\phi_{k}$ are respectively the
complex channel coefficient, delay, Doppler and angle of arrival (AoA) of the
$k$-th user. Notice that, due to the co-location of BS and radar receiver, the
angles of departure and arrival coincide. For simplicity, we focus on uniform
linear arrays (ULA), such that their array response vectors have elements
given by

$\displaystyle[{\bf a}(\phi_{k})]_{i}=e^{j\pi(i-1)\sin(\phi_{k})},\quad i\in
0,\dots,N_{\rm tx}-1$ (2) $\displaystyle[{\bf
b}(\phi_{k})]_{i}=e^{j\pi(i-1)\sin(\phi_{k})},\quad i\in 0,\dots,N_{\rm rx}-1$
(3)

According to the radar equation [17], the channel coefficient $h_{k}$
satisfies

$\displaystyle\lvert h_{k}\rvert^{2}=\frac{\lambda^{2}\sigma_{{\rm
rcs},k}}{(4\pi)^{3}d^{4}_{k}},$ (4)

where $\lambda$ is the wavelength at the central (carrier) frequency, and
$\sigma_{{\rm rcs},k}$ and $d_{k}$ are respectively the radar cross section
(RCS) and range of the $k$-th user at a given epoch.

###  II-B Signal model

We consider OFDM as our modulation scheme since it is one of the standardized
waveforms for mmWave systems and due to its good applicability as an ISAC
waveform [4]. In particular, we use the OFDM pulse shape

$\displaystyle p_{n,m}(t)={\sf rect}\left(\frac{t-nT_{\rm 0}}{T_{\rm
0}}\right)e^{j2\pi m\Delta f(t-T_{\rm cp}-nT_{\rm 0})}\,,$ (5)

where ${\sf rect}(x)$ is a function taking value $1$ when $0\leq x\leq 1$ and
$0$ elsewhere, $\Delta f$ is the subcarrier spacing, $T_{\rm cp}$ is the
cyclic prefix duration, and $T_{\rm
0}\stackrel{{\scriptstyle\Delta}}{{=}}1/\Delta f+T_{\rm cp}$ is the total
symbol duration including cyclic prefix. Considering $N$ symbols and $M$
subcarriers, the transmitted OFDM frame is given by

$\displaystyle{\bf s}(t)=\frac{1}{\sqrt{K}}\sum_{k=1}^{K}{\bf
f}(\hat{\phi}_{k})\sum_{n=0}^{N-1}\sum_{m=0}^{M-1}\zeta_{k}[n,m]p_{n,m}(t)$
(6)

where $\\{\hat{\phi}_{k}\\}$ is a set of predicted AoAs, ${\bf
f}:[-\pi/2,\pi,2]\mapsto\mbox{\bb C}^{N_{\rm tx}}$ is a beamforming function
that generates a unit norm beam pointing towards the intended user, and
$\zeta_{k}[n,m]$ is the $n$-th information symbol intended to user $k$ sent
over subcarrier $m$. For simplicity, we considered $K\leq N_{\rm tx}^{\rm rf}$
in (6). As it is common in the literature [18], we will make the assumption

$\displaystyle\Delta f\gg\nu_{\rm max},$ (7)

where $\nu_{\rm max}$ is the maximum Doppler frequency to be expected in the
channel.

The received backscattered signal in the antenna plane and in the absence of
noise is given by

$\displaystyle{\bf y}(t)$
$\displaystyle=\frac{1}{\sqrt{K}}\sum_{k=1}^{K}h_{k}{\bf b}(\phi_{k}){\bf
a}^{\scriptscriptstyle\mathsf{H}}(\phi_{k}){\bf
s}(t-\tau_{k})e^{j2\pi\nu_{k}t}$ (8)
$\displaystyle=\frac{1}{\sqrt{K}}\sum_{k=1}^{K}h_{k}{\bf
b}(\phi_{k})\sum_{k^{\prime}=1}^{K}{\bf
a}^{\scriptscriptstyle\mathsf{H}}(\phi_{k}){\bf f}(\hat{\phi}_{k^{\prime}})$
$\displaystyle\qquad\qquad\sum_{n=0}^{N-1}\sum_{m=0}^{M-1}\zeta_{k^{\prime}}[n,m]p_{n,m}(t-\tau_{k})e^{j2\pi\nu_{k}t}$
(9) $\displaystyle\approx\frac{1}{\sqrt{K}}\sum_{k=1}^{K}h_{k}{\bf
b}(\phi_{k}){\bf a}^{\scriptscriptstyle\mathsf{H}}(\phi_{k}){\bf
f}(\hat{\phi}_{k})$
$\displaystyle\qquad\qquad\sum_{n=0}^{N-1}\sum_{m=0}^{M-1}\zeta_{k}[n,m]p_{n,m}(t-\tau_{k})e^{j2\pi\nu_{k}t},$
(10)

where the last step follows from the approximation $\lvert{\bf
a}^{\scriptscriptstyle\mathsf{H}}(\phi_{k}){\bf
f}(\hat{\phi}_{k^{\prime}})\rvert\approx 0$ for $k^{\prime}\neq k$, which is
accurate in massive MIMO systems when the users are spatially separated and
the predictions are close to the true value [19].

Aiming to reduce hardware complexity and energy consumption at the radar
receiver, we process the received signal ${\bf y}(t)$ by a reduction matrix
before sampling. In order to be able to produce super-resolution angle
estimates, we let the radar receiver focus on a single user $k^{*}$ for each
OFDM frame, and estimate different users sequentially in time. We achieve this
by tuning the reduction matrix ${\bf U}_{k^{*}}\in\mbox{\bb C}^{N_{\rm
rx}\times N_{\rm rx}^{\rm rf}}$ as will be described later. Then, after
standard OFDM processing (see e.g. [3]) and including noise, the sampled
signal when focusing on user $k^{*}$ is given by

$\displaystyle{\bf y}_{k^{*}}[n,m]$ $\displaystyle={\bf
U}_{k^{*}}^{\scriptscriptstyle\mathsf{H}}\bigg{(}\frac{1}{\sqrt{K}}\sum_{k=1}^{K}h_{k}e^{j2\pi(nT_{\rm
0}\nu_{k}-m\Delta f\tau_{k})}{\bf b}(\phi_{k})$
$\displaystyle\qquad\qquad\;\;{\bf
a}^{\scriptscriptstyle\mathsf{H}}(\phi_{k}){\bf
f}(\hat{\phi}_{k})\zeta_{k}[n,m]+{\bf w}[n,m]\bigg{)}$ (11)
$\displaystyle={\bf
U}_{k^{*}}^{\scriptscriptstyle\mathsf{H}}\left(\sum_{k=1}^{K}h_{k}g_{t,k}{\bf
b}(\phi_{k})\tilde{\zeta}_{k}[n,m]+{\bf w}[n,m]\right)$ (12)
$\displaystyle\approx{\bf
U}_{k^{*}}^{\scriptscriptstyle\mathsf{H}}h_{k^{*}}g_{t,k^{*}}{\bf
b}(\phi_{k^{*}})\tilde{\zeta}_{k^{*}}[n,m]+{\bf
U}_{k^{*}}^{\scriptscriptstyle\mathsf{H}}{\bf w}[n,m],$ (13)

where in (12) we defined $g_{{\rm
tx},k}\stackrel{{\scriptstyle\Delta}}{{=}}\frac{1}{\sqrt{K}}{\bf
a}^{\scriptscriptstyle\mathsf{H}}(\phi_{k}){\bf f}(\hat{\phi}_{k})$ and
$\tilde{\zeta}_{k}[n,m]\stackrel{{\scriptstyle\Delta}}{{=}}\zeta_{k}[n,m]e^{j2\pi(nT_{\rm
0}\nu_{k}-m\Delta f\tau_{k})}$, ${\bf w}[n,m]\in\mbox{\bb C}^{N_{\rm rx}}$ is
white Gaussian noise with variance $\sigma_{n}^{2}$, and the approximation in
(13) follows from designing ${\bf U}_{k^{*}}$ such that $\|{\bf
U}_{k^{*}}^{\scriptscriptstyle\mathsf{H}}{\bf b}(\phi_{k})\|\approx 0$ for
$k\neq k^{*}$.

##  III Proposed Tracking Scheme

In this section, we describe the two main components of our tracking scheme.
First, we estimate parameters of interest from the sampled signal ${\bf
y}_{k^{*}}[n,m]$. For our tracking scheme, we are interested in the ranges
$\\{d_{l,k}\\}$ and angles of arrival $\\{\phi_{l,k}\\}$ of the different
users, where $l$ indexes measurements sequentially in time. However, due to
the nature of our estimator, we also produce estimates of the Doppler
frequencies $\\{\nu_{l,k}\\}$, which could be used for other purposes outside
of the scope of this paper. Then, we update the set of predicted angles
$\\{\phi_{l,k}\\}$ based on the history of estimations. The overall procedure
is summarized in Algorithm 1.

Algorithm 1 Beam Tracking

Refresh period $\Delta T$, predicted angles $\\{\hat{\phi}_{0,k}\\}$ at $t=0$.

$l\leftarrow 0$

loop

for $k^{*}=1,\dots,K$ do111Iteration $k^{*}$ of the inner loop corresponds to
$t\in[l\Delta T+(k^{*}-1)NT_{\rm 0},l\Delta T+k^{*}NT_{\rm 0})$.

Transmit beamformed OFDM frames towards

all directions $\\{\hat{\phi}_{l,k}\\}$ in the predicted set

as shown in Section II.

Obtain estimate $(\check{d}_{l,k^{*}},\check{\phi}_{l,k^{*}})$ as shown in

Section III-A.

Predict $\hat{\phi}_{l+1,k^{*}}$ as shown in Section III-B.

end for

$l\leftarrow l+1$

Radar receiver idle until $t=l\Delta T$

end loop

###  III-A Radar Parameter Estimation

For this section, we will focus on estimating the parameters of user $k^{*}$
at the $l$-th measurement round and will drop the indices to simplify
notation. The measurement model is extensively described in [15], therefore
only a brief summary will be provided here.

We define the reduction matrix ${\bf U}$ as

$\displaystyle{\bf U}={\bf D}(\hat{\phi})\hbox{\boldmath$\Psi$},$ (14)

where ${\bf
D}(\hat{\phi})\stackrel{{\scriptstyle\Delta}}{{=}}{\hbox{diag}}(1,e^{j\pi\sin(\hat{\phi})},\dots,e^{j\pi(N_{\rm
rx}-1)\sin(\hat{\phi})})$ is a tunable network of phase shifters and $\Psi$ is
a fixed network of beamformers pointing towards the broadside direction (i.e.
$\phi=0$) , satisfying ${\bf U}^{\scriptscriptstyle\mathsf{H}}{\bf
U}=\frac{1}{N_{\rm rx}^{\rm rf}}{\bf I}_{N_{\rm rx}^{\rm rf}}$, where ${\bf
I}_{N}$ is the identity matrix of rank $N$. In particular, we let the columns
of $\Psi$ be obtained as the first $N_{\rm rx}^{\rm rf}$ Slepian sequences of
length $N_{\rm rx}$ and any user defined time-bandwidth product. A brief
description of Slepian sequences can be found in [15] and more extensive
treatment in [20].

Using this design, we can obtain the sample covariance matrix from $N\times M$
samples

$\displaystyle\hat{{\bf R}}$
$\displaystyle=\frac{1}{NM}\sum_{n=0}^{N-1}\sum_{m=0}^{M-1}{\bf y}[n,m]{\bf
y}^{\scriptscriptstyle\mathsf{H}}[n,m]$ $\displaystyle\approx P_{\rm tx}\lvert
h\rvert^{2}\lvert g_{\rm
tx}\rvert^{2}\hbox{\boldmath$\Psi$}^{\scriptscriptstyle\mathsf{H}}{\bf
b}(\phi^{\prime}){\bf
b}^{\scriptscriptstyle\mathsf{H}}(\phi^{\prime})\hbox{\boldmath$\Psi$}+\sigma_{n}^{2}{\bf
I}_{N_{\rm rx}^{\rm rf}},$ (15)

where $P_{\rm tx}\stackrel{{\scriptstyle\Delta}}{{=}}\mbox{\bb
E}\left[\lvert\zeta[n,m]\rvert^{2}\right]$ is the average power of the
transmitted signal and
$\phi^{\prime}\stackrel{{\scriptstyle\Delta}}{{=}}\sin^{-1}(\sin(\phi)-\sin(\hat{\phi}))$.
Given the structure of $\hat{{\bf R}}$, we can process it via spectral methods
such as MUSIC [20] in order to obtain $\check{\phi}$ as a super-resolution
estimate of $\phi$.

Finally, we can use the signal model in (13) together with our estimate of
$\phi$ to obtain a least squares estimate of $\tau$ and $\nu$. To do so, we
solve

$\displaystyle(\check{\tau},\check{\nu})=\underset{(\tau,\nu)}{{\hbox{arg}}\min}\;\;\sum_{n,m}\left\|{\bf
y}[n,m]-hg_{\rm tx}{\bf U}^{\scriptscriptstyle\mathsf{H}}{\bf
b}(\check{\phi})\tilde{\zeta}[n,m]\right\|_{2}^{2}\;\;,$ (16)

which, as shown in [15], is equivalent to solving

$\displaystyle(\check{\tau},\check{\nu})=\underset{(\tau,\nu)}{{\hbox{arg}}\max}\;\;\left|\sum_{n,m}y^{\prime}[n,m]\zeta^{*}[n,m]e^{-j2\pi(nT_{\rm
0}\nu-m\Delta f\tau)}\right|,$ (17)

where $y^{\prime}[n,m]\stackrel{{\scriptstyle\Delta}}{{=}}\frac{{\bf
b}^{\scriptscriptstyle\mathsf{H}}(\check{\phi}){\bf U}{\bf y}[n,m]}{{\bf
b}^{\scriptscriptstyle\mathsf{H}}(\check{\phi}){\bf U}{\bf
U}^{\scriptscriptstyle\mathsf{H}}{\bf b}(\check{\phi})}$. Problem (17) can be
efficiently solved by applying FFT and finding a peak in a 2D grid. We can
estimate the range $d$ as

$\displaystyle\check{d}=\frac{c\check{\tau}}{2},$ (18)

where $c$ is the speed of light.

###  III-B Next Angle Prediction

In this section, we summarize the key points of the beam prediction step in
[11] and propose a combination of their simple tracking equations with the
radar receiver introduced in Section III-A that yields an overall improved
performance. By focusing still on a single user, we will introduce the time
index $l$ to capture the prediction based on the past observations.

The approach in [11] aims to predict the coordinates of the user of interest
given the kinematic state in the last three measurement epochs. In particular,
focusing on the $x$ coordinate,

$\displaystyle\begin{cases}&x_{l+1}-x_{l}=v_{x,l}\Delta T+a_{x,l}\Delta
T^{2}/2,\\\ &x_{l}-x_{l-1}=v_{x,l-1}\Delta T+a_{x,l-1}\Delta T^{2}/2,\\\
&x_{l-1}-x_{l-2}=v_{x,l-2}\Delta T+a_{x,l-2}\Delta T^{2}/2,\\\
&v_{x,l}-v_{x,l-1}=a_{x,l-1}\Delta T,\\\ &v_{x,l-1}-v_{x,l-2}=a_{x,l-2}\Delta
T\end{cases}$ (19)

where $x_{l}$, $v_{x,l}$ and $a_{x,l}$ are respectively the $x$-th coordinate,
the velocity and the acceleration in the $x$-th coordinate at the $l$-th
measurement step, and $\Delta T$ is the interval between measurements
corresponding to the same user. In this work, we choose $\Delta T\gg T_{\rm
0}$ for two reasons. First, since different users are estimated sequentially
in time, the number of users this system can accommodate is upper bounded by
$\lfloor\Delta T/T_{\rm 0}\rfloor$, where $\lfloor x\rfloor$ indicates the
largest integer not greater than $x$. Second, $T_{\rm 0}$ can be too short for
any non-negligible motion to happen. Therefore, using a much longer $\Delta T$
does not significantly affect performance, while greatly simplifying hardware
operation by reducing the rate at which beams must be reconfigured.

By setting the assumption that the acceleration does not significantly change
in three measurement steps (i.e. $a_{x,l-2}\approx a_{x,l-1}\approx a_{x,l}$),
the system of equations in (19) can be solved for $x_{l+1}$ yielding

$\displaystyle\hat{x}_{l+1}=3x_{l}-3x_{l-1}+x_{l-2}\approx
3\check{x}_{l}-3\check{x}_{l-1}+\check{x}_{l-2},$ (20)

where, keeping our notation consistent, $\check{x}$ represents an estimated
value (e.g. by following the processing in Section III-A)222Notice that
estimates of $x$ and $y$ can be directly obtained from estimates of $d$ and
$\phi$ using basic trigonometry., and $\hat{x}$ represents a predicted value.

The derivation holds verbatim for the $y$-th coordinate, so we can obtain a
prediction for the next AoA from the relation between the $x$ and $y$
coordinates as

$\displaystyle\hat{\phi}_{l+1}=\tan^{-1}\left(\frac{3\check{x}_{l}-3\check{x}_{l-1}+\check{x}_{l-2}}{3\check{y}_{l}-3\check{y}_{l-1}+\check{y}_{l-2}}\right).$
(21)

##  IV Numerical Results

We compare our approach with the one proposed in [11] that assumes a fully
digital receiver (i.e. $N_{\rm rx}^{\rm rf}=N_{\rm rx}$) and considers
estimation via local linear approximation [11, Section III.B]. We remark here
that our method, by operating in the beam-space domain, overcomes two major
issues associated with fully digital architectures. First, it greatly reduces
hardware complexity since the wideband signal is sampled by only $N_{\rm rf}$
A/D converters, where we can in general have $N_{\rm rf}\ll N_{\rm a}$ as will
be shown in this section. Second, it should be noted that the raw samples at
the antenna outputs are extremely noisy due to the large pathloss of the two-
way mmWave channel. Our method projects the high dimensional signal into a
suitable low dimensional subspace where most of the information content is
preserved but noise from uninteresting directions is rejected, thus resulting
in a much higher SNR per sample. Slepian sequences provide a systematic way to
define such a subspace by procuring an orthonormal basis upon which signals
coming from directions of interest have maximal projection, whereas those
coming from any other angular region are approximately orthogonal. We also
note that we had to assume perfect matched filtering for the method in [11] to
perform adequately, whereas our method is robust to small processing errors
(e.g. those caused by discrete sampling in time and frequency). A detailed
description of what this assumption entails can be found in Appendix A. For
the sake of completeness, we also include performance evaluation for the
method in [11] when no perfect matched filtering is assumed. Furthermore, for
ease of evaluation, in the presented results we focus on a single user
scenario (i.e. $K=1$). Table I summarizes the parameters used in our
simulations.

TABLE I: System parameters $N=64$ | $M=512$
---|---
$N_{\rm tx}=N_{\rm rx}=64$ | $N_{\rm tx}^{\rm rf}=N_{\rm rx}^{\rm rf}=4$
$f_{\rm c}=$ $60\text{\,}\mathrm{GHz}$ |  $\Delta f=$ $1\text{\,}\mathrm{MHz}$
$d_{\rm max}=$ $100\text{\,}\mathrm{m}$ |  $v_{\rm max}=$ $30\text{\,}\mathrm{m}\mathrm{/}\mathrm{s}$
$P_{\rm TX}$ = $50\text{\,}\mathrm{mW}$ |  $\sigma_{\rm rcs}=$ $20\text{\,}\mathrm{dBsm}$
Noise PSD: $N_{0}$ = $2\text{\times}{10}^{-21}\text{\,}\mathrm{W}\mathrm{/}\mathrm{Hz}$ |  $\Delta T=$ $100\text{\,}\mathrm{ms}$
Figure 1: Main road scenario and trajectory considered.

050100150200250Epoch number0246810Achievable spectral efficiency
[bps/Hz]Oracle[11][11], perfect MFProposed Approach

Figure 2: Achievable spectral efficiency as a function of the time step for
the selected trajectory.

$-$100$-$75$-$50$-$250255075100x [m]020406080100y [m]Appraoch in
[11]$-$100$-$75$-$50$-$250255075100x [m]020406080100y [m]Proposed
Approach012345678Spectral efficiency [bps/Hz]

Figure 3: Achievable spectral efficiency as a function of the position in the
road for each of the two trackers considered.

To verify our results in arbitrary road geometries, we consider the road
structure shown in Figure1. In particular, for our first result, we fix the
path shown in blue and study the performance of the two trackers in terms of
achievable spectral efficiency as a function of the time step, averaged over
500 Monte Carlo simulations. The achievable spectral efficiency assumes a
single antenna user, and is computed as

$\displaystyle\eta_{l}\approx\mbox{\bb
E}\left[\log_{2}\left(1+\left(\frac{\lambda}{4\pi
d_{l}}\right)^{2}\frac{P_{\rm tx}|g_{{\rm tx},l}|^{2}}{N_{0}(M\Delta
f)}\right)\right],$ (22)

where the approximation follows from using an empirical estimation of the
expectation and $N_{0}$ is the noise power spectral density at the user. The
results are presented in Figure2, and show how our method is able to quickly
recover from the non-smooth changes of direction while the approach in [11]
loses track of the user and is not able to find it again.

In order to average results over different mobility patterns, we generate
trajectories where a user moves within the road structure defined in Figure1,
but now is able to take any of the paths, and where the speed at each time
step is sampled from a random process. Since now there is no association
between time step and kinematic state, we illustrate performance as a function
of the position by means of a heat map. This is shown in Figure3, where our
method clearly outperforms the baseline in [11] when averaging over an
ensemble of different trajectories. In particular, it can be observed that the
approach in [11] fails most of the time to track the user in the locations
closest to the BS, where the angles change at a faster rate, even though the
path loss there is small. On the other hand, the performance of our proposed
method depends only on the distance (or equivalently the path loss) to the
user, but is robust to fast changes in the trajectory.

Figure 4: Average beamforming gain as a function of the number of RF chains
for different numbers of antenna elements.

Finally, we study the effect of the number of RF chains in the tracking
performance. To do so, we generate multiple trajectories and average the
beamforming gain obtained at each point of each path for different antenna
front end configurations. This is presented in Figure4, where it can be seen
that as few as 4 RF chains are enough to achieve close to optimal performance.
This result suggests that projecting the high dimensional signal at the
antenna plane onto a much lower dimensional space preserves most of the
information necessary for estimation, provided that the considered subspace is
adequately selected. The figure also shows how the gap between the achieved
results and the upper bound becomes larger when the number of antennas
increases. This can be explained by the fact that larger arrays generate
narrower beams, resulting on a higher sensitivity to small pointing errors.

##  V Conclusions

In this paper, we have shown an efficient method to perform beam tracking in
OFDM ISAC systems, based on the combination of an advanced HDA receiver and a
simple tracking equation. Our close to optimal results indicate that large
savings in computation and hardware complexity can be obtained without
sacrificing performance. Moreover, the validation of our results in
complicated road geometries suggests that our method is general enough to
perform well in realistic mobility scenarios.

## References

  * [1] G. Flagship, “6G white paper on localization and sensing,” _University of Oulu, Finland_ , June 2020.
  * [2] D. Ma, N. Shlezinger, T. Huang, Y. Liu, and Y. C. Eldar, “Joint radar-communication strategies for autonomous vehicles: Combining two key automotive technologies,” _IEEE Signal Process. Mag._ , vol. 37, no. 4, pp. 85–97, 2020.
  * [3] C. Sturm and W. Wiesbeck, “Waveform design and signal processing aspects for fusion of wireless communications and radar sensing,” _Proc. IEEE_ , vol. 99, no. 7, pp. 1236–1259, 2011.
  * [4] J. A. Zhang, F. Liu, C. Masouros, R. W. H. J. au2, Z. Feng, L. Zheng, and A. Petropulu, “An overview of signal processing techniques for joint communication and radar sensing,” 2021.
  * [5] I. A. Hemadeh, K. Satyanarayana, M. El-Hajjar, and L. Hanzo, “Millimeter-wave communications: Physical channel models, design considerations, antenna constructions, and link-budget,” _IEEE Commun. Surveys Tuts._ , vol. 20, no. 2, pp. 870–913, 2018.
  * [6] W. Yuan, F. Liu, C. Masouros, J. Yuan, D. W. K. Ng, and N. González-Prelcic, “Bayesian predictive beamforming for vehicular networks: A low-overhead joint radar-communication approach,” _IEEE Transactions on Wireless Communications_ , vol. 20, no. 3, pp. 1442–1456, 2021.
  * [7] S. G. Larew and D. J. Love, “Adaptive beam tracking with the unscented kalman filter for millimeter wave communication,” _IEEE Signal Processing Letters_ , vol. 26, no. 11, pp. 1658–1662, 2019.
  * [8] V. Va, H. Vikalo, and R. W. Heath, “Beam tracking for mobile millimeter wave communication systems,” in _2016 IEEE Global Conference on Signal and Information Processing (GlobalSIP)_ , 2016, pp. 743–747.
  * [9] J. Lim, H.-M. Park, and D. Hong, “Beam tracking under highly nonlinear mobile millimeter-wave channel,” _IEEE Communications Letters_ , vol. 23, no. 3, pp. 450–453, 2019.
  * [10] S. K. Dehkordi, M. Kobayashi, and G. Caire, “Adaptive beam tracking based on recurrent neural networks for mmwave channels,” in _2021 IEEE 22nd International Workshop on Signal Processing Advances in Wireless Communications (SPAWC)_ , 2021, pp. 1–5.
  * [11] F. Liu and C. Masouros, “A tutorial on joint radar and communication transmission for vehicular networks—part iii: Predictive beamforming without state models,” _IEEE Commun. Lett._ , vol. 25, no. 2, pp. 332–336, 2021.
  * [12] Z. Du, F. Liu, W. Yuan, C. Masouros, Z. Zhang, S. Xia, and G. Caire, “Integrated sensing and communications for v2i networks: Dynamic predictive beamforming for extended vehicle targets,” 2021.
  * [13] S.-H. Hyun, J. Song, K. Kim, J.-H. Lee, and S.-C. Kim, “Adaptive beam design for v2i communications using vehicle tracking with extended kalman filter,” _IEEE Transactions on Vehicular Technology_ , vol. 71, no. 1, pp. 489–502, 2022.
  * [14] C. Liu, W. Yuan, S. Li, X. Liu, D. W. K. Ng, and Y. Li, “Learning-based predictive beamforming for integrated sensing and communication in vehicular networks,” 2021.
  * [15] F. Pedraza, M. Kobayashi, and G. Caire, “Beam refinement and user state acquisition via integrated sensing and communication with ofdm,” in _2021 IEEE 22nd International Workshop on Signal Processing Advances in Wireless Communications (SPAWC)_ , 2021, pp. 476–480.
  * [16] X. Song, S. Haghighatshoar, and G. Caire, “Efficient beam alignment for millimeter wave single-carrier systems with hybrid MIMO transceivers,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 3, pp. 1518–1533, 2019.
  * [17] M. A. Richards, _Fundamentals of radar signal processing_. Tata McGraw-Hill Education, 2005.
  * [18] L. Gaudio, M. Kobayashi, B. Bissinger, and G. Caire, “Performance analysis of joint radar and communication using ofdm and otfs,” in _2019 IEEE International Conference on Communications Workshops (ICC Workshops)_ , 2019, pp. 1–6.
  * [19] T. L. Marzetta, E. G. Larsson, H. Yang, and H. Q. Ngo, _Fundamentals of Massive MIMO_. Cambridge University Press, 2016.
  * [20] P. Stoica and R. Moses, _Spectral Analysis of Signals_. Pearson Prentice Hall, 2005.

##  Appendix A The perfect matched filtering assumption

In this appendix, we illustrate the limitation of the assumption on perfect
matched filtering frequently used in recent works. To that aim, let us focus
for simplicity in a single path channel model, such that the time-varying
channel matrix is given by

$\displaystyle{\bf H}(t,\tau)=h_{0}{\bf b}(\phi_{0}){\bf
a}^{\scriptscriptstyle\mathsf{H}}(\phi_{0})\delta(\tau-\tau_{0})e^{j2\pi\nu_{0}t},$
(23)

where $h_{0}$, $\phi_{0}$, $\tau_{0}$ and $\nu_{0}$ represent respectively the
channel coefficient, AoA, delay and Doppler associated with a single user. Let
us consider a generic transmitted signal ${\bf
s}(t)\stackrel{{\scriptstyle\Delta}}{{=}}{\bf f}s(t)$ time limited to the
interval $t\in[0,T]$, where ${\bf f}$ is a beamforming vector, such that the
received signal in the absence of noise is

$\displaystyle{\bf r}(t)=h_{0}g(\phi_{0}){\bf
b}(\phi_{0})s(t-\tau_{0})e^{j2\pi\nu_{0}t},$ (24)

where $g(\phi_{0})\stackrel{{\scriptstyle\Delta}}{{=}}{\bf
a}^{\scriptscriptstyle\mathsf{H}}(\phi_{0}){\bf f}$. The (ideal) matched
filtering process would then consist in processing the received signal as

$\displaystyle{\bf r}$ $\displaystyle=\int_{-\infty}^{\infty}{\bf
r}(t^{\prime})s^{*}(t^{\prime}-(t-T))e^{-j2\pi\nu
t^{\prime}}dt^{\prime}\bigg{|}_{\begin{subarray}{c}t=T+\tau_{0}\\\
\nu=\nu_{0}\end{subarray}}$ (25) $\displaystyle=h_{0}g(\phi_{0}){\bf
b}(\phi_{0})$
$\displaystyle\qquad\int_{-\infty}^{\infty}s(t^{\prime}-\tau_{0})s^{*}(t^{\prime}-(t-T))e^{-j2\pi(\nu_{0}-\nu)t^{\prime}}dt^{\prime}\bigg{|}_{\begin{subarray}{c}t=T+\tau_{0}\\\
\nu=\nu_{0}\end{subarray}}$ (26) $\displaystyle=h_{0}g(\phi_{0}){\bf
b}(\phi_{0})\int_{-\infty}^{\infty}s(t^{\prime}-\tau_{0})s^{*}(t^{\prime}-\tau_{0})dt^{\prime}$
(27) $\displaystyle=h_{0}{\cal E}_{s}g(\phi_{0}){\bf b}(\phi_{0}),$ (28)

where ${\cal
E}_{s}\stackrel{{\scriptstyle\Delta}}{{=}}\int_{-\infty}^{\infty}|s(t)|^{2}dt$
is the energy of signal $s(t)$. Recent works consider a noisy version of (28)
as a model for the sampled signal. However, as illustrated here, explicit
knowledge of the delay and Doppler of the channel was required to arrive to
this expression. This assumption becomes unrealistic in mobility settings
where delay and Doppler might be time varying magnitudes. In practice, most
systems sample a grid of $(\tau,\nu)$ tuples and find the pair maximizes the
norm of the output of the matched filter. Let us refer to such a pair as
$(\hat{\tau},\hat{\nu})$, where $(\hat{\tau},\hat{\nu})\neq(\tau_{0},\nu_{0})$
almost surely due to the finiteness of the grid. Then, the exact matched
filter output would be given by

$\displaystyle{\bf r}=h_{0}g(\phi_{0}){\bf
b}(\phi_{0})\int_{-\infty}^{\infty}s(t^{\prime}-\tau_{0})s^{*}(t^{\prime}-\hat{\tau})e^{j2\pi(\nu_{0}-\hat{\nu})t^{\prime}}dt^{\prime},$
(29)

which shows that the term ${\cal E}_{s}$ should be replaced with an expression
depending on
$\Delta\tau\stackrel{{\scriptstyle\Delta}}{{=}}\tau_{0}-\hat{\tau}$ and
$\Delta\nu\stackrel{{\scriptstyle\Delta}}{{=}}\nu_{0}-\hat{\nu}$. Notice also
that $\Delta\tau$ and $\Delta\nu$ can vary quickly in time since $\tau_{0}$
($\nu_{0}$) changes smoothly and $\hat{\tau}$ ($\hat{\nu}$) changes in
discrete steps. When using model (29) in lieu of (28), many of the works in
the literature become unusable. In particular, works that consider that the
sampled signal is described by a statistical distribution that is a perfectly
known function of the AoA and the channel coefficient fail to consider this
difficult to model phenomenon. Also, methods that take derivatives of the
model with respect to the parameters in order to linearize the functional
dependence, like the baseline used in this paper, have typically considered
the formulation in (28), resulting in simple expressions that are however
invalid once the perfect matched filtering assumption is lifted.

  *[OFDM]: orthogonal frequency division multiplexing
  *[ISAC]: Integrated Sensing and Communication
  *[HDA]: hybrid digital analog
  *[mmWave]: millimeter wave
  *[BS]: Base Station
  *[LoS]: line-of-sight


∎

e1e-mail<EMAIL_ADDRESS>e2e-mail<EMAIL_ADDRESS>e3e-mail:
<EMAIL_ADDRESS>e4e-mail<EMAIL_ADDRESS>
11institutetext: INFN Laboratori Nazionali di Frascati, Via Enrico Fermi, 40,
00044 Frascati, Italy

# Impacts of radiative corrections on measurements of lepton flavour
universality in ${B}\rightarrow{D}\ell{{\nu}_{\ell}}$ decays

Stefano Calíe1,addr1  Suzanne Klavere2,addr1  Marcello Rotondoe3,addr1
Barbara Sciasciae4,addr1

(Received: 28 May 2019 / Accepted: 27 August 2019)

###### Abstract

Radiative corrections to ${B}\rightarrow{D}\ell{{\nu}_{\ell}}$ decays may have
an impact on predictions and measurements of the lepton flavour universality
observables $\mathcal{R}({{D}^{+}})$ and $\mathcal{R}({{D}^{0}})$. In this
paper, a comparison between recent calculations of the effect of soft-photon
corrections on $\mathcal{R}({{D}^{+}})$ and $\mathcal{R}({{D}^{0}})$, and
corrections generated by the widely used package Photos is given. The impact
of long-distance Coulomb interactions, which are not simulated in Photos, is
discussed. Furthermore, the effect of high-energy photon emission is studied
through pseudo-experiments in an LHCb-like environment. It is found that over-
or underestimating these emissions can cause a bias on $\mathcal{R}({D})$ as
high as 7%. However, this bias depends on individual analyses, and future high
precision measurements require an accurate evaluation of these QED
corrections.

††journal: Eur. Phys. J. C

##  1 Introduction

The Standard Model (SM) assumes lepton universality (LU) implying that once
the mass difference is taken into account, all SM interactions treat the three
charged leptons identically. The mass difference results in a different phase
space between decays involving $\tau^{-}$ and the lighter $e^{-}$ and
$\mu^{-}$ leptons111Throughout this paper, the inclusion of charge-conjugate
processes is implied and natural units with $\hbar=c=1$ are used.. LU can be
tested by measuring the ratio of decay rates, ensuring that the Cabibbo-
Kobayashi-Maskawa matrix elements, as well as most of the form factors, cancel
in the ratio. This results in more accurate theoretical predictions and in the
cancellation of many experimental systematic uncertainties. One type of these
LU measurements is performed using semileptonic $B$ decays of the form
${b}\rightarrow{c}{\ell^{-}}{{\overline{\nu}}_{\ell}}$, commonly known as
measurements of $\mathcal{R}({H_{c}})$, defined as

$\displaystyle{\mathcal{R}({H_{c}})}$
$\displaystyle=\frac{{\mathcal{B}}({H_{b}}\rightarrow{H_{c}}{\tau^{-}}{{\overline{\nu}}_{\tau}})}{{\mathcal{B}}({H_{b}}\rightarrow{H_{c}}\ell^{-}{{\overline{\nu}}_{\ell}})}\,,$
(1)

where $H_{b}$ and $H_{c}$ are a $b$ and $c$ hadron, respectively, and $\ell$
is either an electron or muon.

Several measurements of $\mathcal{R}({H_{c}})$ have been performed by the
LHCb, Belle and BaBar experiments. For $\mathcal{R}({D})$, on which this paper
is focused, the predicted value Bigi:2016mdz ; Bernlochner:2017jka ;
Jaiswal:2017rve ; Aoki:2019cca  is

$\displaystyle{\mathcal{R}({D})}={\mathcal{R}({{D}^{+}})}={\mathcal{R}({{D}^{0}})}=0.299\pm
0.003\,,$ (2)

which assumes isospin symmetry. The average of the measured value of
$\mathcal{R}({D})$ is $0.349\pm 0.027\pm 0.015$ Lees:2012xj ; Huschle:2015rga
; Abdesselam:2019dgh , where the first uncertainty is statistical and the
second systematic. Even though $\mathcal{R}({D})$ differs from the SM
prediction by only 1.4$\sigma$, it is remarkable that the deviation from the
SM of the combined $\mathcal{R}({D})$ and $\mathcal{R}({{D}^{*}})$ observables
is 3.1$\sigma$ HFLAV16 .

Radiative corrections were long thought to be negligible at the level of
precision of measurements and predictions of $\mathcal{R}({D})$. Recently,
however, de Boer et al. deBoer:2018ipi  presented a new evaluation of the
long-distance electromagnetic (QED) contributions to ${{\kern
1.79993pt\overline{\kern-1.79993ptB}{}}{}^{0}}\rightarrow{{D}^{+}}{\ell^{-}}{{\overline{\nu}}_{\ell}}$
and ${{{B}^{-}}}\rightarrow{{D}^{0}}{\ell^{-}}{{\overline{\nu}}_{\ell}}$
decays, where ${\ell^{-}}={\mu^{-}},{\tau^{-}}$. They point out that these
soft-photon corrections are different for $\mu^{-}$ and $\tau^{-}$ decays,
such that they do not cancel in the ratios $\mathcal{R}({{D}^{+}})$ and
$\mathcal{R}({{D}^{0}})$. According to the authors of Ref. deBoer:2018ipi ,
the proper evaluation of the radiative corrections alters the SM predictions
of the $\mathcal{R}({D})$ and $\mathcal{R}({{D}^{*}})$ values and increases
their uncertainty. The current tension between the SM and measurements could
be weakened or strengthened if radiative corrections are not properly taken
into account.

All experiments measuring these types of LU are dependent on the simulation of
QED radiative corrections in decays of particles and resonances. The widely
used package to simulate these corrections is Photos Barberio:1993qi ;
Golonka:2005pn , which is used by all three experiments measuring
$\mathcal{R}({D})$ and $\mathcal{R}({{D}^{*}})$.

This paper starts by comparing the radiative corrections on
$\mathcal{R}({{D}^{+}})$ and $\mathcal{R}({{D}^{0}})$ from Ref. deBoer:2018ipi
with those simulated by Photos in Sect. 2. The sensitivity of measurements of
$\mathcal{R}({{D}^{+}})$ and $\mathcal{R}({{D}^{0}})$ to radiative corrections
in the $\mu^{-}$ and $\tau^{-}$ decay modes is studied with pseudo-experiments
in an LHCb-like environment, with different assumptions on the shape of the
total energy of the radiated photons. The method and the results of this study
are reported in Sect. 3. Conclusions and recommendations are summarised in
Sect. 4.

##  2 Radiative corrections in PHOTOS

Photos Golonka:2005pn ; Golonka:1379813  is a universal Monte Carlo algorithm
that simulates the effects of QED corrections in decays of particles and
resonances. It exploits the factorisation property of QED coming from the
exponentiation method used to improve the convergence of the perturbative
expansion. Any particle-decay process accompanied by bremsstrahlung photons
can be factorised into a tree term and bremsstrahlung factor. The latter
depends only on the four-momenta of those particles taking part in the decay,
and not on the underlying process. This approximation, which takes into
account both real and virtual corrections, converges to an exact expression in
the soft-photon region of phase space. It is worth noting that Photos does not
incorporate the emission of photons depending on the hadronic structure. These
so called structure-dependent (SD) photons impact the spin of the decay
particle, and may also interfere with bremsstrahlung photons. The effect of SD
photons depends on the specific decay under study and, as was the case for
kaon decays Bijnens:1994me , may not be negligible.

The latest versions of Photos include multi-photon emissions, and interference
between final-state photons. The validity of Photos has been tested
successfully by comparing its results to full calculations available in
various processes involving $W$, $Z$ and hadronic $B$ decays into scalar
mesons Golonka:2005pn ; Nanava:2006vv . Because of the universal treatment of
photon emission in Photos, its performances in specific processes should
always be checked, especially when high precision is desired or when signal
extraction is sensitive to detailed simulation of a phase space corner of the
decay.

The calculation by de Boer et al. in Ref. deBoer:2018ipi  is the first that
studies the impact of soft-photon corrections on $\mathcal{R}({{D}^{+}})$ and
$\mathcal{R}({{D}^{0}})$. It is valid in the regime in which the maximum
energy of the radiated photons is smaller than the lepton mass, which is the
muon mass in this case. This calculation includes more effects than Photos
does, in particular the interference between initial- and final-state photons,
and the Coulomb correction. The latter increases the decay rate of decays with
charged particles in the final state. It should be noted that the contribution
of the Coulomb correction is singular for null relative velocity between
final-state charged particles.

To compare QED corrections between Photos and Ref. deBoer:2018ipi , four
samples (${{\kern
1.79993pt\overline{\kern-1.79993ptB}{}}{}^{0}}\rightarrow{{D}^{+}}{\ell^{-}}{{\overline{\nu}}_{\ell}}$
and ${{{B}^{-}}}\rightarrow{{D}^{0}}{\ell^{-}}{{\overline{\nu}}_{\ell}}$,
where ${\ell^{-}}={\mu^{-}},{\tau^{-}}$) with three million $B$-meson decays
are generated by Pythia 8 Sjostrand:2006za ; Sjostrand:2007gs . The decays are
simulated by EvtGen Lange:2001uf , and the radiative corrections by Photos
v.3.56, with the “option with interference” switched on. QED corrections are
applied by Photos by modifying the charged track’s four-momentum in the event
record filled by EvtGen every time a photon is added.

The four-momentum of the total radiated photons, $\mbox{$p$}_{\gamma}$, is
defined as

$\mbox{$p$}_{\gamma}=\mbox{$p$}_{B}-\left(\mbox{$p$}_{{D}}+\mbox{$p$}_{{\ell^{-}}}+\mbox{$p$}_{{{\overline{\nu}}_{\ell}}}\right)\,,$
(3)

where $\mbox{$p$}_{B}$, $\mbox{$p$}_{{D}}$,
$\mbox{$p$}_{{{\overline{\nu}}_{\ell}}}$, and $\mbox{$p$}_{{\ell^{-}}}$ are
the four-momenta of the $B$, $D$, $\ell^{-}$ and ${\overline{\nu}}_{\ell}$
particles, respectively, taken from the event record updated by Photos. This
means that, in agreement with Ref. deBoer:2018ipi , the radiation of the $D$
decay products is not taken into account. The total energy of the radiated
photons, $E_{\gamma}$, is computed in the $B$ rest frame. As in Ref.
deBoer:2018ipi , the variable $E_{\rm{max}}$ is defined as the maximum value
that $E_{\gamma}$ is allowed to have to consider
${B}\rightarrow{D}\ell{\bar{\nu}}_{\ell}(\gamma)$ as signal.

Figure 1:  Distribution of the total energy of the radiated photons,
$E_{\gamma}$, up to 100 MeV for ${\kern
1.79993pt\overline{\kern-1.79993ptB}{}}{}^{0}$ $\rightarrow$ ${D}^{+}$
$\tau^{-}$ ${\overline{\nu}}_{\tau}$ and ${\kern
1.79993pt\overline{\kern-1.79993ptB}{}}{}^{0}$ $\rightarrow$ ${D}^{+}$
$\mu^{-}$ ${\overline{\nu}}_{\mu}$ decays as simulated by Photos. Once a value
of $E_{\rm{max}}$ is chosen, e.g. 40 MeV as in the plot, all events with
higher $E_{\gamma}$ values are discarded.

The QED correction, $\delta_{\rm{QED}}$, is given by the relative variation of
the branching ratio when events with total radiated energy greater than
$E_{\rm{max}}$ are discarded. This can be calculated as follows:

${\delta_{\rm{QED}}}=\frac{\int_{0}^{{E_{\rm{max}}}}N(E_{\gamma})dE_{\gamma}}{\int_{0}^{\infty}N(E_{\gamma})dE_{\gamma}}-1\,,$
(4)

where $N(E_{\gamma})$ is the distribution of events with $E_{\gamma}$. This
distribution is shown for ${\kern
1.79993pt\overline{\kern-1.79993ptB}{}}{}^{0}$ $\rightarrow$ ${D}^{+}$
$\tau^{-}$ ${\overline{\nu}}_{\tau}$ and ${\kern
1.79993pt\overline{\kern-1.79993ptB}{}}{}^{0}$ $\rightarrow$ ${D}^{+}$
$\mu^{-}$ ${\overline{\nu}}_{\mu}$ decays in Fig. 1. The considered energy
range is up to 100 MeV, which covers the majority of radiative photons, namely
98% of the $\mu^{-}$ decays and 99.7% for the $\tau^{-}$ decays generated by
Photos.

Comparisons between radiative corrections from Photos and Ref. deBoer:2018ipi
are shown in Fig. 2 for the ${{\kern
1.79993pt\overline{\kern-1.79993ptB}{}}{}^{0}}\rightarrow{{D}^{+}}{\ell^{-}}{{\overline{\nu}}_{\ell}}$
(left panel) and
${{{B}^{-}}}\rightarrow{{D}^{0}}{\ell^{-}}{{\overline{\nu}}_{\ell}}$ (middle
panel) branching fractions. These plots show differences of up to 2% for
${\kern 1.79993pt\overline{\kern-1.79993ptB}{}}{}^{0}$ decays, and $0.5-1\%$
for ${{B}^{-}}$ decays. Unfortunately, the effect does not cancel in the
ratios of branching fractions. This is clearly visible from the right panel of
Fig. 2 where radiative corrections on $\mathcal{R}({D})$,
$\delta_{QED}(\mathcal{R})$, are shown as a function of $E_{\rm{max}}$. Photos
predicts a QED correction that is 0.5% lower than the one in Ref.
deBoer:2018ipi  for $\mathcal{R}({{D}^{+}})$, while it is 0.5% higher than the
one in Ref. deBoer:2018ipi  for $\mathcal{R}({{D}^{0}})$.

Figure 2:  Radiative corrections to the branching ratios of
$\bar{B}^{0}\rightarrow D^{+}\ell^{-}\bar{\nu_{\ell}}$ (left) and
$B^{-}\rightarrow D^{0}\ell^{-}\bar{\nu_{\ell}}$ (middle) decays, as a
function of $E_{\rm{max}}$. The long-distance QED corrections to
$\mathcal{R}({{D}^{+}})$ (orange) and $\mathcal{R}({{D}^{0}})$ (violet) as a
function of $E_{\rm{max}}$ (right). The plots are obtained from simulated data
(solid lines, which include the statistical uncertainty) and from Ref.
deBoer:2018ipi  (dashed lines, filled with transparent colours when the
uncertainties are significant).

###  2.1 Coulomb correction

A significant part of the radiative corrections in Ref. deBoer:2018ipi
originates from Coulomb interactions, which are not included in Photos. Note
that the Coulomb correction is relevant for the ${D}^{+}$ mode, but not for
the ${D}^{0}$ mode. For a fermion-scalar (and fermion-fermion) pair, this
correction is given by

${\Omega_{\rm{C}}}=\frac{2\pi\alpha}{\beta_{D\ell}}\frac{1}{1-e^{-\frac{2\pi\alpha}{\beta_{D\ell}}}}{\mbox{,}}$
(5)

where $\alpha=1/137$ and $\beta_{D\ell}$ is the relative velocity between the
$D$ meson and the lepton, defined as

$\beta_{D\ell}=\left[1-\frac{4m_{D}^{2}m_{\ell}^{2}}{(s_{D\ell}-m_{D}^{2}-m_{\ell}^{2})^{2}}\right]^{1/2}\,,$
(6)

where $s_{D\ell}=(p_{D}+p_{\ell})^{2}$. A well-known approximation of the
Coulomb correction by Atwood and Marciano Atwood:1989em , yields
${\Omega_{\rm{C}}}=(1+\pi\alpha)\,\approx 1.023$ which occurs when
$\beta_{D\ell}\approx 1$. This is accurate for decays with light leptons, but
not for those with $\tau^{-}$ leptons. For the semitauonic mode, the typical
relative velocity is 0.5-0.9, resulting in a Coulomb correction between 2.5
and 5.0%.

QED corrections from Photos for the ${D}^{+}$ mode are also compared with
predictions not including the Coulomb correction from Ref. deBoer:2018ipi .
This reduces the difference of the corrections to the branching ratios between
Photos and the theoretical calculations to about 1% and brings the corrections
on $\mathcal{R}({{D}^{+}})$ in close agreement, as shown in Fig. 3 (left and
middle, respectively).

Figure 3:  Radiative corrections to the branching ratios of
$\bar{B}^{0}\rightarrow D^{+}\ell^{-}\bar{\nu_{\ell}}$ (left) and
$\mathcal{R}({{D}^{+}})$ (middle) in the case that no Coulomb correction is
applied. The plot on the right shows the ratio
${\delta_{\rm{QED}}}({\mathcal{R}({{D}^{+}})})/{\delta_{\rm{QED}}}({\mathcal{R}({{D}^{0}})})$.
The plots are obtained from simulated data (solid lines, which include the
statistical uncertainty) and from Ref. deBoer:2018ipi  (dashed lines, filled
with transparent colours when the uncertainties are significant).

Fig. 3 (right) shows the ratio of QED corrections on $\mathcal{R}({{D}^{+}})$
over those on $\mathcal{R}({{D}^{0}})$. It is worth noting that both Photos
and the calculation in Ref. deBoer:2018ipi  without Coulomb correction
conserve isospin symmetry ($\delta_{\rm{QED}}$ values for
$\mathcal{R}({{D}^{+}})$ and $\mathcal{R}({{D}^{0}})$ agree within the
errors), while the Coulomb correction introduces an isospin-breaking term.

##  3 Effects on LHCb-like analysis

The comparison in the previous section only holds for values of $E_{\gamma}$
up to 100 MeV. For higher energies, no calculations relevant to
$\mathcal{R}({D})$ are available222In Ref.Bernlochner:2010yd  a calculation of
the high-energy SD contribution to
${B}\rightarrow{D}\ell{{\overline{\nu}}_{\ell}}$ is reported. This is not
relevant to this study because of the missing lepton-mass dependent effects..
Nevertheless, Photos generates photons with energies larger than those treated
in Ref. deBoer:2018ipi , which is the range where the effect of SD photons
might be relevant.

To study the effect of under- or overestimating radiative corrections in
simulations used for measurements of $\mathcal{R}({D})$, a simplified analysis
is performed in an LHCb-like environment. Also here, the radiation emitted by
the decay products of the $D$ mesons is neglected because their contributions
largely cancel out in the ratio.

The strategy of this study consists of fitting a data sample with templates
describing the ${B}\rightarrow{D}{\mu^{-}}{{\overline{\nu}}_{\mu}}$ and
${B}\rightarrow{D}{\tau^{-}}{{\overline{\nu}}_{\tau}}$ components. The fits
are performed with templates built under the hypothesis that no radiation with
$E_{\gamma}$ above a certain value $E_{\rm{max}}$ is emitted. In particular,
five $E_{\rm{max}}$ values were chosen to cut on $E_{\gamma}$: 100, 300, 500,
800 and 1500 MeV. The bias on $\mathcal{R}({D})$, determined from these fits,
is an indication of the importance of the simulation of the $E_{\gamma}$
distribution in the high-energy region.

This analysis follows a strategy similar to the one used in Ref. Aaij:2015yra
, where $\mathcal{R}({{D}^{*}})$ is measured using a three-dimensional
templated fit. The data samples, referred to as pseudo-experiments in the
following text, are generated from a mixture of
${B}\rightarrow{D}{\mu^{-}}{{\overline{\nu}}_{\mu}}$ and
${B}\rightarrow{D}{\tau^{-}}{{\overline{\nu}}_{\tau}}$ decays, with radiative
corrections generated by Photos. Here $\mathcal{R}({D})$ is assumed to be 0.3
as predicted by the SM.

The variables used in the templated fit performed to extract the
${B}\rightarrow{D}{\mu^{-}}{{\overline{\nu}}_{\mu}}$ and
${B}\rightarrow{D}{\tau^{-}}{{\overline{\nu}}_{\tau}}$ yields from the pseudo-
experiments are: the muon energy computed in the $B$ meson rest frame,
$E_{\mu}$; the missing mass squared, ${m_{\rm
miss}^{2}}=(\mbox{$p$}_{{B}}-\mbox{$p$}_{{D}}-\mbox{$p$}_{\mu})^{2}$; and the
squared four-momentum transferred to the lepton system,
${q^{2}}=(\mbox{$p$}_{{B}}-\mbox{$p$}_{{D}})^{2}$. The variables are binned as
follows: four bins in $q^{2}$ in the range $-0.4<{q^{2}}<12.6$ GeV, 40 bins in
$m_{\rm miss}^{2}$ between $-2<{m_{\rm miss}^{2}}<10$ GeV2, and 30 bins in
muon energy in the range $100<{E_{\mu}}<2500$ MeV, consistent with Ref.
Aaij:2015yra . In this case study, only the signal
(${B}\rightarrow{D}{\tau^{-}}{{\overline{\nu}}_{\tau}}$) and normalisation
(${B}\rightarrow{D}{\mu^{-}}{{\overline{\nu}}_{\mu}}$) components are
considered, while all backgrounds are ignored.

Basic selection requirements are applied to mimic the acceptance of the LHCb
detector and its trigger following Ref. Ciezarek:2016lqu . Both production and
$B$-decay vertex positions are smeared to simulate the resolution of the LHCb
detector. The resolution on the production vertices is 13${\,\upmu\mathrm{m}}$
in $x$ and $y$, and 70${\,\upmu\mathrm{m}}$ in $z$ direction. For the
$B$-decay vertices a resolution of 20${\,\upmu\mathrm{m}}$ in $x$ and $y$, and
200${\,\upmu\mathrm{m}}$ in $z$ direction is used, after which the $B$
direction is computed. The $D$ mesons decay as
${{D}^{0}}\rightarrow{{K}^{-}}{{\pi}^{+}}$ and
${{D}^{+}}\rightarrow{{K}^{-}}{{\pi}^{+}}{{\pi}^{+}}$, the $\tau^{-}$ lepton
as ${\tau^{-}}\rightarrow{\mu^{-}}{{\overline{\nu}}_{\mu}}{{\nu}_{\tau}}$. The
muons and all decay products from the $D$ mesons are required to be in the
pseudorapidity range between $1.9$ and $4.9$. In addition, the momentum of
each of these particles is required to be larger than $5$ GeV, and its
component transverse to the beam direction must be larger than $250$ MeV. The
distance between the production and $B$-decay vertex should be at least 3 mm,
similar to the requirements applied in a typical trigger selection.

Due to the missing neutrino and unknown effective centre-of-mass energy of the
collision, the $B$-meson momentum cannot be reconstructed in an
$\mathcal{R}({H_{c}})$ analysis at LHCb. Therefore the momentum of the $B$
meson in the $z$ direction, $(\mbox{$p$}_{{B}})_{z}$, is approximated as
$(\mbox{$p$}_{{B}})_{z}=(m_{{B}}/m_{\rm vis})(\mbox{$p$}_{\rm vis})_{z}$,
where $m_{{B}}$ is the $B$ mass, and $m_{\rm vis}$ and $(\mbox{$p$}_{\rm
vis})_{z}$ are the momentum in the $z$ direction and the mass of the visible
decay products of the $B$ meson, respectively. This directly follows the
approach from Ref. Aaij:2015yra . After computing the $B$ momentum with the
above approximation and applying the selection criteria described in this
section, $q^{2}$, $m_{\rm miss}^{2}$ and $E_{\mu}$ are calculated. The
distributions for the signal and control samples are shown in Fig. 4. Even
using this simplified detector description, these distributions show the same
key features as those in Ref. Aaij:2015yra .

Figure 4: Shapes of the $q^{2}$, $m_{\rm miss}^{2}$, and $E_{\mu}$
distributions after applying basic selection requirements and the rest frame
approximation for the ${{\kern
1.79993pt\overline{\kern-1.79993ptB}{}}{}^{0}}\rightarrow{{D}^{+}}{\ell^{-}}{{\overline{\nu}}_{\ell}}$
decays where ${\ell^{-}}={\mu^{-}},{\tau^{-}}$.

When applying cuts on $E_{\gamma}$, the templates shapes change. This is most
clearly seen in the distributions of $m_{\rm miss}^{2}$, shown for the ${\kern
1.79993pt\overline{\kern-1.79993ptB}{}}{}^{0}$ decay in Fig. 5. Especially in
the $\mu^{-}$ decay mode the effect is large, altering the shape at high
values of $m_{\rm miss}^{2}$. Since this feature is not present in the
$\tau^{-}$ mode, this does not cancel when measuring the ratio
$\mathcal{R}({{D}^{+}})$. For completeness, changes in the shape of $q^{2}$
and $E_{\mu}$ for the $\mu^{-}$ mode are shown in A.

Figure 5:  Ratios of the cut $m_{\rm miss}^{2}$ distribution over the default
$m_{\rm miss}^{2}$ distributions for the $\bar{B}^{0}\rightarrow
D^{+}\ell^{-}\bar{\nu_{\ell}}$ decays, for various $E_{\rm{max}}$ cuts. On the
left for the $\mu^{-}$ and on the right for the $\tau^{-}$ mode.

The number of events generated to simulate data is determined from the
estimated number of events that LHCb gathered during their Run II data-taking
period. The estimate takes into account the $B$-production cross-section at 13
TeV, branching fractions, and assumes the average reconstruction efficiency is
the same as in Ref. Aaij:2015yra . This results in data samples of $1.0\times
10^{6}$ and $0.5\times 10^{5}$ for the ${{\kern
1.79993pt\overline{\kern-1.79993ptB}{}}{}^{0}}\rightarrow{{D}^{+}}{\ell^{-}}{{\overline{\nu}}_{\ell}}$
decays, and $4.4\times 10^{5}$ and $2.3\times 10^{4}$ for the
${{{B}^{-}}}\rightarrow{{D}^{0}}{\ell^{-}}{{\overline{\nu}}_{\ell}}$ decays,
where the first yield represents the $\mu^{-}$ sample, and the second the
$\tau^{-}$ sample. In an actual analysis, the efficiencies for
${B}\rightarrow{D}\ell{{\overline{\nu}}_{\ell}}$ decays are likely higher than
those for ${B}\rightarrow{D}^{*}\ell{{\overline{\nu}}_{\ell}}$ decays, where
${D}^{*}$ is reconstructed in the ${D}\pi$ decay mode.

The measured value of $\mathcal{R}({D})$ is determined from two components.
The first is the ratio of reconstruction efficiencies $\varepsilon_{\mu}$ and
$\varepsilon_{\tau}$ for the $\mu^{-}$ and $\tau^{-}$ samples, respectively,
which takes into account the selection requirements described earlier in this
section. The second component is the fraction of semitauonic decays in the
sample, $f_{\tau}$, determined from the three-dimensional template fit (the
absence of background events in the simulated samples implies that the
fraction of $\mu^{-}$ and $\tau^{-}$ components add up to one). These are
combined to measure $\mathcal{R}({D})$ as

$\displaystyle{\mathcal{R}({D})}=\frac{f_{\tau}}{1-f_{\tau}}\frac{\varepsilon_{\mu}}{\varepsilon_{\tau}}\,.$
(7)

The exercise of generating pseudo-experiments is repeated 10.000 times after
which the spread of the measured values of $\mathcal{R}({D})$ is taken as the
statistical uncertainty.

The resulting values of $\mathcal{R}({{D}^{+}})$ and $\mathcal{R}({{D}^{0}})$
as a function of $E_{\rm{max}}$ are shown in Fig. 6. From here it is clear
that there is a significant effect in underestimating the QED radiative
corrections which could be up to 0.02 for both $\mathcal{R}({{D}^{+}})$ and
$\mathcal{R}({{D}^{0}})$ values, corresponding to a relative bias of 7.5%. The
largest contribution to the observed bias is due to the fit fraction
$f_{\tau}$, which is strongly affected by the shapes of the $\mu^{-}$
templates. Instead, the ratio of efficiencies
${\varepsilon_{\mu}}/{\varepsilon_{\tau}}$ is only marginally dependent on
$E_{\rm{max}}$. However, this last statement holds only for this specific case
study. Different sets of selection cuts or different experimental environments
could indeed introduce a significant bias also in the ratio of efficiencies.
The observed bias can be understood when looking at the $m_{\rm miss}^{2}$
distribution in Fig. 4 and Fig. 5. When cutting on $E_{\gamma}$, part of the
tail of the $\mu^{-}$ distribution is removed, which is compensated by a
higher $\tau^{-}$ fraction in the fit.

Figure 6: Values of $\mathcal{R}({D})$ as a function of $E_{\rm{max}}$
measured in a simplified LHCb-like analysis of $\mathcal{R}({{D}^{+}})$ (left)
and $\mathcal{R}({{D}^{0}})$ (right). The error bars reflect the statistical
uncertainties from the generated MC samples. The blue bands correspond to fit
results obtained with no $E_{\rm{max}}$ cuts.

In an actual analysis there are radiative corrections in MC and radiated
photons in data. Therefore, it is useful to check the above approach using an
alternative strategy. In this case, the templates include all QED corrections
predicted by Photos while an $E_{\rm{max}}$ cut is applied on the pseudo-
experiments. This approach leads to an overestimate on the QED corrections,
resulting in a negative bias on the $\mathcal{R}({D})$ values. The results for
$\mathcal{R}({{D}^{+}})$ and $\mathcal{R}({{D}^{0}})$ as a function of
$E_{\rm{max}}$ are reported in A. The corrections are of the exact same size
as those in the baseline approach.

It is worth to note that, despite the fact that LHCb does not cut explicitly
on $E_{\gamma}$ in its analyses, indirect cuts on the total radiated energy
are applied through e.g. requirements on isolation variables or inefficient
reconstruction algorithms for low momentum particles. This could alter the
observed bias if $E_{\gamma}$ is not simulated correctly.

These studies show that radiative corrections play a crucial role in
$\mathcal{R}({D})$ measurements. Since part of these corrections are already
simulated in Photos, the above exercise shows the effect of a worst-case
scenario. Nevertheless, additional effects such as the Coulomb correction, as
detailed in next section, or the calculation for energies greater than 100 MeV
are becoming fundamental in view of the increased experimental precision
expected in the coming years. Also, these quantitative effects strongly depend
on explicit or implicit cuts on radiative photons and must be carefully
evaluated for each analysis measuring $\mathcal{R}({D})$.

###  3.1 Coulomb correction

Beyond affecting the SM prediction of $\mathcal{R}({{D}^{+}})$, the Coulomb
correction impacts the experimental results by changing the shape of the fit
templates. This is evaluated in the ${{\kern
1.79993pt\overline{\kern-1.79993ptB}{}}{}^{0}}\rightarrow{{D}^{+}}{\ell^{-}}{{\overline{\nu}}_{\ell}}$
decay by weighting each event by the term $\Omega_{\rm{C}}$.333The
contribution due to the non-factorizable loop corrections to the tree-level
differential decay rate, called $\tilde{\Gamma}^{{{D}^{+}}}$ in
Ref.deBoer:2018ipi , is small and not implemented in Photos. For this reason,
the Coulomb correction $\Omega_{\rm{C}}$ can be introduced in Photos as a
global factor to the uncorrected differential decay rate. The changes in the
shape of the $q^{2}$, $m_{\rm miss}^{2}$ and $E_{\mu}$ distributions are shown
in Fig. 7. While for the $\mu^{-}$ mode $\Omega_{\rm{C}}$ is mostly constant,
for the $\tau^{-}$ mode there is a dependence on each of the three variables
due to the smaller relative velocity. To quantify the effect of the Coulomb
correction, the above analysis is repeated, without including any
$E_{\rm{max}}$ cuts. The Coulomb correction is applied to the pseudo-
experiments, but not to the fit templates, resulting in a relative shift of
about -1.0% on $\mathcal{R}({{D}^{+}})$. This effect can even be amplified by
selecting certain regions of phase space.

Figure 7:  Coulomb correction as a function of $q^{2}$, $m_{\rm miss}^{2}$,
and $E_{\mu}$ for ${{\kern
1.79993pt\overline{\kern-1.79993ptB}{}}{}^{0}}\rightarrow{{D}^{+}}{\ell^{-}}{{\overline{\nu}}_{\ell}}$
decays, where ${\ell^{-}}={\mu^{-}},{\tau^{-}}$.

##  4 Conclusions and recommendations

The work in Ref. deBoer:2018ipi  describes QED corrections which are not fully
included in Photos. These corrections affect the semimuonic and semitauonic
modes differently at the level of a few percent. Ignoring the Coulomb
correction, there is more radiated energy in the calculation in Ref.
deBoer:2018ipi  than in Photos for the ${\kern
1.79993pt\overline{\kern-1.79993ptB}{}}{}^{0}$ decays, while this is the other
way around for the ${{B}^{-}}$ decays. In the ratio $\mathcal{R}({D})$, this
small discrepancy mostly cancels out. However, the main difference between the
QED corrections on $\mathcal{R}({{D}^{+}})$ and $\mathcal{R}({{D}^{0}})$,
which is up to 1%, is due to the Coulomb correction that only affects
$\mathcal{R}({{D}^{+}})$.

Coulomb interactions are not simulated by Photos and mainly affect the
kinematics of semitauonic decays, which in turn influence the shape of
distributions used to determine the signal yields in LHCb, BaBar, and Belle
analyses. These effects can alter values of $\mathcal{R}({D})$ up to 1% in an
LHCb-like analysis, and should be evaluated precisely for each measurement.

Using a simplified LHCb-like analysis, it is shown that over- or
underestimating radiative corrections could bias measurements of
$\mathcal{R}({D})$ up to 7% in an extreme case. This results in a bias of 0.02
on the value of $\mathcal{R}({D})$, and should be studied further when
performing these types of measurements, including a realistic evaluation of
cuts on $E_{\gamma}$. These effects could potentially be enhanced in
measurements from Belle II Kou:2018nap  where the resolution on the kinematic
variables is better than at LHCb.

When measuring values of $\mathcal{R}({D})$ with higher precision, additional
calculations of QED corrections for ${B}\rightarrow{D}\ell{{\nu}_{\ell}}$
decays are necessary. Especially calculations involving high-energy and
structure-dependent photons are currently mostly missing.

###### Acknowledgements.

We are grateful to S. de Boer, T. Kitahara, and I. Nisandzic for the fruitful
collaboration, and to U. Egede for his thoughtful comments. In addition, we
thank the Semileptonics $B$ decays working group of the LHCb collaboration,
and in particular M. De Cian and L. Grillo, for their useful feedback
throughout the development of this paper. Finally, we thank Z. Was for his
helpful insights into the Photos package.

##  Appendix A Additional Plots

The fits on the pseudo-experiments are performed on the three variables
$E_{\mu}$, $m_{\rm miss}^{2}$ and $q^{2}$, as described in Sect. 3. The effect
of cutting on $E_{\rm{max}}$ on the shape of the $E_{\mu}$ and $m_{\rm
miss}^{2}$ templates is shown in Fig. 8 for the ${{\kern
1.79993pt\overline{\kern-1.79993ptB}{}}{}^{0}}\rightarrow{{D}^{+}}{\mu^{-}}{{\overline{\nu}}_{\mu}}$
decay. Analogous plots for the semitauonic mode show a negligible dependence
on the $E_{\rm{max}}$ cut.

Figure 8: Ratios of the cut $q^{2}$ (left) and $E_{\mu}$ (right) distributions
over the corresponding default distributions for ${{\kern
1.79993pt\overline{\kern-1.79993ptB}{}}{}^{0}}\rightarrow{{D}^{+}}{\mu^{-}}{{\overline{\nu}}_{\mu}}$
decays, for various $E_{\rm{max}}$ cuts.

The results of performing the simplified LHCb-like analysis with the
alternative strategy are shown in Fig. 9. These results are obtained using
templates with an $E_{\gamma}$ distribution in agreement with Photos
predictions, and pseudo-experiments with cuts on $E_{\gamma}$ applied.

Figure 9: Values of $\mathcal{R}({D})$ as a function of $E_{\rm{max}}$ from a
simplified LHCb-like analysis of $\mathcal{R}({{D}^{+}})$ (left) and
$\mathcal{R}({{D}^{0}})$ (right). In this alternative approach, templates are
generated according to Photos predictions and $E_{\rm{max}}$ cuts are applied
to the pseudo-experiments. The error bars reflect the statistical
uncertainties from the generated MC samples. The blue bands correspond to fit
results obtained with no $E_{\rm{max}}$ cuts.

## References

  * (1) D. Bigi and P. Gambino, _Revisiting $B\rightarrow D\ell\nu$_, Phys. Rev. D94 (2016) 094008, arXiv:1606.08030
  * (2) F. U. Bernlochner, Z. Ligeti, M. Papucci, and D. J. Robinson, _Combined analysis of semileptonic $B$ decays to $D$ and $D^{*}$: $R(D^{(*)})$, $|V_{cb}|$, and new physics_, Phys. Rev. D95 (2017) 115008, arXiv:1703.05330
  * (3) S. Jaiswal, S. Nandi, and S. K. Patra, _Extraction of $|V_{cb}|$ from $B\rightarrow D^{(*)}\ell\nu_{\ell}$ and the Standard Model predictions of $R(D^{(*)})$_, JHEP 12 (2017) 060, arXiv:1707.09977
  * (4) Flavour Lattice Averaging Group, S. Aoki et al., _FLAG Review 2019_ , arXiv:1902.08191
  * (5) BaBar, J. P. Lees et al., _Evidence for an excess of $\bar{B}\rightarrow D^{(*)}\tau^{-}\bar{\nu}_{\tau}$ decays_, Phys. Rev. Lett. 109 (2012) 101802, arXiv:1205.5442
  * (6) Belle, M. Huschle et al., _Measurement of the branching ratio of $\bar{B}\rightarrow D^{(\ast)}\tau^{-}\bar{\nu}_{\tau}$ relative to $\bar{B}\rightarrow D^{(\ast)}\ell^{-}\bar{\nu}_{\ell}$ decays with hadronic tagging at Belle_, Phys. Rev. D92 (2015) 072014, arXiv:1507.03233
  * (7) Belle, A. Abdesselam et al., _Measurement of $\mathcal{R}(D)$ and $\mathcal{R}(D^{\ast})$ with a semileptonic tagging method_, arXiv:1904.08794
  * (8) Heavy Flavor Averaging Group, Y. Amhis et al., _Averages of $b$-hadron, $c$-hadron, and $\tau$-lepton properties as of summer 2016_, Eur. Phys. J. C77 (2017) 895, arXiv:1612.07233, updated results and plots available at https://hflav.web.cern.ch
  * (9) S. de Boer, T. Kitahara, and I. Nisandzic, _Soft-Photon Corrections to $\bar{B}\rightarrow D\tau^{-}\bar{\nu}_{\tau}$ Relative to $\bar{B}\rightarrow D\mu^{-}\bar{\nu}_{\mu}$_, Phys. Rev. Lett. 120 (2018) 261804, arXiv:1803.05881
  * (10) E. Barberio and Z. Was, _PHOTOS: A Universal Monte Carlo for QED radiative corrections. Version 2.0_ , Comput. Phys. Commun. 79 (1994) 291
  * (11) P. Golonka and Z. Was, _PHOTOS Monte Carlo: A precision tool for QED corrections in $Z$ and $W$ decays_, Eur. Phys. J. C45 (2006) 97, arXiv:hep-ph/0506026
  * (12) P. Golonka, Computer simulations in high energy physics: a case for PHOTOS, MC-TESTER, TAUOLA and at2sim, PhD thesis, Cracow, INP, 2006, CERN-THESIS-2006-098
  * (13) J. Bijnens, G. Colangelo, G. Ecker, and J. Gasser, _Semileptonic kaon decays_ , in 2nd DAPHNE Physics Handbook:315-389, pp. 315–389, 1994.  arXiv:hep-ph/9411311
  * (14) G. Nanava and Z. Was, _Scalar QED, NLO and PHOTOS Monte Carlo_ , Eur. Phys. J. C51 (2007) 569, arXiv:hep-ph/0607019
  * (15) T. Sjöstrand, S. Mrenna, and P. Skands, _PYTHIA 6.4 physics and manual_ , JHEP 05 (2006) 026, arXiv:hep-ph/0603175
  * (16) T. Sjöstrand, S. Mrenna, and P. Skands, _A brief introduction to PYTHIA 8.1_ , Comput. Phys. Commun. 178 (2008) 852, arXiv:0710.3820
  * (17) D. J. Lange, _The EvtGen particle decay simulation package_ , Nucl. Instrum. Meth. A462 (2001) 152
  * (18) D. Atwood and W. J. Marciano, _Radiative Corrections and Semileptonic $B$ Decays_, Phys. Rev. D41 (1990) 1736
  * (19) F. U. Bernlochner and H. Lacker, _A phenomenological model for radiative corrections in exclusive semileptonic B-meson decays to (pseudo)scalar […]_ , arXiv:1003.1620
  * (20) LHCb, R. Aaij et al., _Measurement of the ratio of branching fractions $\mathcal{B}(\bar{B}^{0}\rightarrow D^{*+}\tau^{-}\bar{\nu}_{\tau})/\mathcal{B}(\bar{B}^{0}\rightarrow D^{*+}\mu^{-}\bar{\nu}_{\mu})$_, Phys. Rev. Lett. 115 (2015) 111803, arXiv:1506.08614, [Erratum: Phys. Rev. Lett.115,no.15,159901(2015)]
  * (21) G. Ciezarek, A. Lupato, M. Rotondo, and M. Vesterinen, _Reconstruction of semileptonically decaying beauty hadrons produced in high energy pp collisions_ , JHEP 02 (2017) 021, arXiv:1611.08522
  * (22) Belle-II, W. Altmannshofer et al., _The Belle II Physics Book_ , arXiv:1808.10567


MIDLMedical Imaging with Deep Learning 2019 MIDL 2019 – Extended Abstract
Track Azam Hamidinekoo1Corresponding Author<EMAIL_ADDRESS>
Erika Denton2<EMAIL_ADDRESS>
Reyer Zwiggelaar1<EMAIL_ADDRESS>
1 Department of Computer Science, Aberystwyth University, Aberystwyth, United
Kingdom
2 Department of Radiology, Norfolk & Norwich University Hospital, Norwich,
United Kingdom

# Automated Mammogram Analysis with a Deep Learning Pipeline

###### Abstract

Current deep learning based detection models tackle detection and segmentation
tasks by casting them to pixel or patch-wise classification. To automate the
initial mass lesion detection and segmentation on the whole mammographic
images and avoid the computational redundancy of patch-based and sliding
window approaches, the conditional generative adversarial network (cGAN) was
used in this study. Subsequently, feeding the detected regions to the trained
densely connected network (DenseNet), the binary classification of benign
versus malignant was predicted. We used a combination of publicly available
mammographic data repositories to train the pipeline, while evaluating the
model’s robustness toward our clinically collected repository, which was
unseen to the pipeline.

###### keywords:

Mammography, Generative Adversarial Network, DenseNet, Detection,
Segmentation, Classification

##  1 Introduction

Breast mass lesions are mostly dense and appear in grey to white pixel
intensity values on mammograms, with various size, distribution, shape and
density. Computer-Aided Detection (CADe) and Diagnosis (CADi) systems based on
deep learning (DL) methods have been improved with regards to their
performances but still cannot identify all cancerous cases [Hamidinekoo et
al.(2018a)Hamidinekoo, Denton, Rampun, Honnor, and Zwiggelaar]. DL can be
applied directly on the whole image or in a patch-based approach. Majority of
the approaches proposed for mass detection [Hamidinekoo et
al.(2018a)Hamidinekoo, Denton, Rampun, Honnor, and Zwiggelaar] use a patch-
based method for training and a sliding window approach for testing. These
methods basically involve patch classification followed by localisation, which
are time-consuming with substantial redundancy for overlapping patches. To
avoid the computational redundancy of these approaches, we proposed using the
cGAN [Isola et al.(2017)Isola, Zhu, Zhou, and Efros] to increase the
efficiency by training on whole mammographic images.

The contribution of this work was to develop an automatic detection and
segmentation model (CADe model) using images from multi-centres and detect
mass abnormalities on our clinically collected unseen mammograms. To complete
our pipeline, the detected regions were fed into the DenseNet [Huang et
al.(2017)Huang, Liu, Van Der Maaten, and Weinberger] (CADi model) to predict
the benign or malignant nature of the detected regions.

##  2 Methodology

###  2.1 Datasets

Four multi-centred, publicly available mammographic databases were used to
build the CADe and CADi models. The combination of them is referred to as
Set-1. After building the CADe and CADi models, their performances were
evaluated using an unseen dataset, called Set-2. This dataset was clinically
collected from the Norfolk & Norwich University hospital, UK. The information
about these datasets are summarised in Table 1. In the pre-processing step,
all images were segmented into background and tissue and the intensity values
of the segmented regions were normalised [Hamidinekoo et
al.(2018b)Hamidinekoo, Suhail, Denton, and Zwiggelaar].

For building the CADe, an additional dataset was created from the raw images
in Set-1, where Gaussian noise with $\sigma=2$ was added to increase the
number of samples.

For building the CADi, Regions of Interest (RoIs) were extracted with size
equal to double the square bounding box of the lesion to include the
neighbourhood information. These extractions were scaled to 256$\times$256
followed by random crops of 224$\times$224, which were used for training the
CADi [Hamidinekoo et al.(2017)Hamidinekoo, Suhail, Qaiser, and Zwiggelaar].

Table 1: Utilised databases containing mass lesions.

|  Set-1: used for training the CADe and CADi models |  Set-2: used for testing
---|---|---
|  BCDR-F03 [Lopez et al.(2012)Lopez, Posada, Moura, Pollán, Valiente, Ortega, Solar, Diaz-Herrero, Ramos, Loureiro, Fernandes, and Ferreira de Araújo] |  BCDR-D01 [Lopez et al.(2012)Lopez, Posada, Moura, Pollán, Valiente, Ortega, Solar, Diaz-Herrero, Ramos, Loureiro, Fernandes, and Ferreira de Araújo] |  DDSM [Heath et al.(2001)Heath, Bowyer, Kopans, Moore, and Kegelmeyer] |  Inbreast [Moreira et al.(2012)Moreira, Amaral, Domingues, Cardoso, Cardoso, and Cardoso] | Private-Dataset
Number of cases |  341 |  51 |  975 |  102 | 103
Number of images |  664 |  105 |  1930 |  102 | 210
Benign images |  369 |  69 |  1023 |  34 | 46
Malignant images |  295 |  36 |  907 |  68 | 164
Resolution (bits/pixel) |  8 |  14 |  12, 16 |  14 | 12, 14
Image mode |  digitised |  digital |  digitised |  digital | digital
View |  MLO, CC |  MLO, CC |  MLO, CC |  MLO, CC | MLO, CC
Age distribution |  58.4$\pm$15.3 |  57.7$\pm$13.5 |  58.9$\pm$11.5 |  - | 60.9$\pm$17.7

###  2.2 Model Architecture

#### CADe development

Mass detection and segmentation on mammograms can be defined as translating
each image into the corresponding semantic label map representing mass lesion
in the breast tissue. Therefore, we propose using the conditional GAN (cGAN)
[Isola et al.(2017)Isola, Zhu, Zhou, and Efros] in order to apply a specific
condition on the input images to train a conditional generative model for the
CADe. In the architecture of the utilised cGAN, a U-Net-based structure and a
convolutional PatchGAN classifier were used as the generator and the
discriminator, respectively. Both the generator and the discriminator used
modules of the form convolution-BatchNorm-ReLu. The generator was trained to
generate mask images from the mammographic images, which were expected to be
similar to the mask images of the real observed images from Set-1 samples. The
cGAN-based CADe model was trained on the prepared training set via the
stochastic gradient descent (SGD) solver with learning rate=0.00004 and batch-
size=16 (using parameter tuning) for epochs=50. When several regions were
detected, connected component analysis was used in the post-processing section
and the three largest detected regions were extracted as RoIs.

#### CADi development

Based on a comparative study [Hamidinekoo et al.(2018b)Hamidinekoo, Suhail,
Denton, and Zwiggelaar], among the well-known deep CNNs, the DenseNet was
found as an appropriate model for mass classification due to its key
characteristic to bypass signals from the preceding layers to the subsequent
layers. In our implementations, the DenseNet’s growth rate was set to 4 to
construct 4 dense-blocks and 3 transition layers in the architecture. In this
model, the final Softmax classifier made a binary decision based on the
created features. The rest of the model’s parameters (kernel, stride and
padding sizes) were kept as default [Huang et al.(2017)Huang, Liu, Van Der
Maaten, and Weinberger]. The objective of training for the CADi was to
minimise the difference error between the network prediction and the expected
output (benign vs. malignant). DenseNet was trained via the SGD solver with
Gamma=0.1, momentum=0.9 and weight-decay=10-5 along with batch-size=8 (based
on our hardware specifications), a dynamic learning rate with initial
value=0.001 for 30 epochs. We used transfer learning with the ImageNet
dataset, whilst the network was fine-tuned using Set-1. The trained CADi was
able to classify images in the validation set with the accuracy of 76% and 90%
and the AUC of 0.78 and 0.87 on the combination of digitised+digital images
and only digital images, respectively.

##  3 Results & Discussions

To evaluate CADe performance, each mammographic image (from Set-2) was fed to
the model and the probability map corresponding to the probable detected
lesion was computed and compared with their annotated images. Using the CADe,
mass abnormalities on mammographic images were detected with 34% accuracy, 34%
precision, 32% recall and an F.Score of 0.33. The detected RoIs were segmented
with Dice Similarity Coefficient=0.33$\pm$0.30 and Hausdorff
Distance=8.33$\pm$2.49. Considering these results, the detection and
segmentation values did not look good. This was caused due to several
predictable reasons: (1) Mass segmentation is not explicitly undertaken in
regular breast screening. So, the available annotations for the Set-1 data
repositories and Set-2 samples were prone to subjectivity, which was not
possible to consider due to the lack of such information; (2) All the lesion
boundaries that were provided by the Set-1, were roughly made. But the
annotated contours that were provided for the Set-2 samples were delicate and
very fine in texture and structure. The model was learned to find a rough
boundary not a very fine contour, which led to low DCS and large HD values;
(3) All the training images used in the Set-1, had single mass lesion on each
individual image but there were a significant number of testing samples (from
the Set-2) with multiple lesions annotated on them; (4) There were several new
appearances in the testing samples that were not seen in the training set,
like micro-calcification deposits or prosthesis as shown in Figure 1; (5) The
Set-2 images were originally in DICOM format and using a thresholding approach
were converted to .png format. In our experiments we discovered that the value
of the threshold applied to this conversion was very important, which could
affect the detection performance; (6) Comparing the lesion size distributions
of the annotated vs predicted lesions illustrated that the model was able to
detect smaller normalised lesion sizes in the range of (0,.1]. As shown in the
examples provided in Figure 1, some lesions were partially detected but were
not considered as an accepted performance in the evaluation section because of
the low dice score (DCS $<$ 0.5). According to these reasons, many of the
segmented bounding boxes (in the testing set) were not accurately aligned with
the annotated masses, but the preliminary results on 103 subjects showed the
promise of this cGAN model for initial lesion localisation to be combined in
the traditional image processing techniques (i.e. initial seed point for the
region growing approach).

Table 2 compares confusion matrices of image based classification results on
double the square bounding box of the detected lesions (through
$CADe_{(\sigma=2)}$) with the actual pre-detected mass lesions by our
radiologist for Set-2. The quantitative results illustrated that the trained
CADi could classify pre-detected mass patches with the accuracy of 73.3%,
which demonstrated (1) its representational capacity to learn different
abnormality features from various digital and digitised training samples; (2)
its robustness and generalisability for an unseen data repository for the task
of classification. However, comparing the classification performance for the
pre-detected lesions vs $CADe_{(\sigma=2)}$ detected regions, comparative
accuracy values were not obtained (73.33% on 210 images vs 41.90% on 139
images, respectively). This suggested that the model was better to be used for
the classification of the lesions that were initially localised and segmented
by the radiologist. Considering that during the reading of screening
mammograms, radiologists use multiple views, we combined the results of
different views for each patient. Giving the priority to the malignant
predictions, Table 3 states the diagnosis output. These results justified the
use of radiologists’ annotated regions for gaining better output from this
pipeline.

Table 2: Comparing image-based classification performance using the model
detected regions vs. annotated lesions by the radiologist for 210 mammograms
in Set-2.

|  | Predicted Label |  |  |
---|---|---|---|---|---
Patches From | Actual Label | Benign | Malignant | Per-class Accuracy | Not Detected | Prediction Accuracy
$CADe_{(\sigma=2)}$ | Benign | 15 | 10 | 60.00% | 21 | 41.90%
| Malignant | 41 | 73 | 64.03% | 50
Annotated lesions | Benign | 36 | 10 | 78.26% | - | 73.33%
Malignant | 46 | 118 | 71.95% | -

Table 3: Comparing subject-based classification performance using the model
detected regions vs. annotated lesions by the radiologist for 103 subjects
with biopsy proven diagnosis.

|  | Predicted Label |  |  |
---|---|---|---|---|---
Patches From | Actual Label | Benign | Malignant | Per-class Accuracy | Wrong Detection & Diagnosis | Diagnosis Accuracy
$CADe_{(\sigma=2)}$ | Benign | 11 | 6 | 64.70% | 5 | 63.10%
| Malignant | 18 | 54 | 75.00% | 9
Annotated lesions | Benign | 9 | 11 | 45.00% | - | 78.64%
Malignant | 11 | 72 | 86.74% | -

Figure 1: Examples of various detection, segmentation and classification
performances on Set-2. From left to right: 3 failed detections due to low DSC
and 2 accepted detections and classifications from class malignant and benign.

##  4 Conclusion

Our contribution in this paper was to implement an automated pipeline based on
cGAN and DenseNet models along with the required specifications to initially
analyse whole mammographic images. This model avoided the computational
expense of currently used patch-based or sliding window approaches, commonly
used for large size images (i.e. mammograms). This exploratory research work
could be further extended to using DICOM data instead of converting them into
different formats in order to keep the wealth of original captured data and
decrease the model sensitivity to intensity variations.

## References

  * [Hamidinekoo et al.(2017)Hamidinekoo, Suhail, Qaiser, and Zwiggelaar] Azam Hamidinekoo, Zobia Suhail, Talha Qaiser, and Reyer Zwiggelaar.  Investigating the effect of various augmentations on the input data fed to a convolutional neural network for the task of mammographic mass classification.  In _MIUA_ , pages 398–409, 2017.
  * [Hamidinekoo et al.(2018a)Hamidinekoo, Denton, Rampun, Honnor, and Zwiggelaar] Azam Hamidinekoo, Erika Denton, Andrik Rampun, Kate Honnor, and Reyer Zwiggelaar.  Deep learning in mammography and breast histology, an overview and future trends.  _Medical Image Analysis_ , 47:45–67, 2018a.
  * [Hamidinekoo et al.(2018b)Hamidinekoo, Suhail, Denton, and Zwiggelaar] Azam Hamidinekoo, Zobia Suhail, Erika Denton, and Reyer Zwiggelaar.  Comparing the performance of various deep networks for binary classification of breast tumours.  In _IWBI_ , volume 10718, pages 10710–10718, 2018b.
  * [Heath et al.(2001)Heath, Bowyer, Kopans, Moore, and Kegelmeyer] Michael Heath, Kevin Bowyer, Daniel Kopans, Richard Moore, and W Philip Kegelmeyer.  Digital database for screening mammography.  In _IWDM_ , pages 212–218. Medical Physics, 2001.
  * [Huang et al.(2017)Huang, Liu, Van Der Maaten, and Weinberger] Gao Huang, Zhuang Liu, Laurens Van Der Maaten, and Kilian Q Weinberger.  Densely connected convolutional networks.  In _CVPR_ , pages 4700–4708, 2017.
  * [Isola et al.(2017)Isola, Zhu, Zhou, and Efros] Phillip Isola, Jun-Yan Zhu, Tinghui Zhou, and Alexei A Efros.  Image-to-image translation with conditional adversarial networks.  In _CVPR_ , pages 1125–1134, 2017.
  * [Lopez et al.(2012)Lopez, Posada, Moura, Pollán, Valiente, Ortega, Solar, Diaz-Herrero, Ramos, Loureiro, Fernandes, and Ferreira de Araújo] MA Guevara Lopez, N Posada, Daniel C Moura, Raúl Ramos Pollán, José M Franco Valiente, César Suárez Ortega, M Solar, G Diaz-Herrero, IMAP Ramos, J Loureiro, Teresa Cardoso Fernandes, and Bruno M. Ferreira de Araújo.  BCDR: a breast cancer digital repository.  In _15th International Conference on Experimental Mechanics_ , 2012.
  * [Moreira et al.(2012)Moreira, Amaral, Domingues, Cardoso, Cardoso, and Cardoso] Inês C Moreira, Igor Amaral, Inês Domingues, António Cardoso, Maria João Cardoso, and Jaime S Cardoso.  Inbreast: toward a full-field digital mammographic database.  _Academic Radiology_ , 19(2):236–248, 2012\.