\section{Introduction}

Hadronic light-by-light (HLbL) scattering currently gives the second-largest contribution to
the uncertainty in the Standard Model prediction for the anomalous magnetic moment of the muon~\cite{Aoyama:2020ynm,Aoyama:2012wk,Aoyama:2019ryr,Czarnecki:2002nt,Gnendiger:2013pva,Davier:2017zfy,Keshavarzi:2018mgv,Colangelo:2018mtw,Hoferichter:2019gzf,Davier:2019can,Keshavarzi:2019abf,Hoid:2020xjs,Kurz:2014wya,Melnikov:2003xd,Colangelo:2014dfa,Colangelo:2014pva,Colangelo:2015ama,Masjuan:2017tvw,Colangelo:2017qdm,Colangelo:2017fiz,Hoferichter:2018dmo,Hoferichter:2018kwz,Gerardin:2019vio,Bijnens:2019ghy,Colangelo:2019lpu,Colangelo:2019uex,Blum:2019ugy,Colangelo:2014qya}
\begin{equation}
\label{amuSM}
a_\mu^\text{SM}=116\,591\,810(43)\times 10^{-11}.
\end{equation}
While at present the uncertainty of hadronic vacuum polarization dominates in the comparison to experiment~\cite{bennett:2006fi,Abi:2021gix,Albahri:2021ixb,Albahri:2021kmg,Albahri:2021mtf}
\begin{equation}
\label{exp}
a_\mu^\text{exp}=116\,592\,061(41)\times 10^{-11},
\end{equation}
see Refs.~\cite{Borsanyi:2020mff,Lehner:2020crt,Crivellin:2020zul,Keshavarzi:2020bfy,Malaescu:2020zuc,Colangelo:2020lcg} for recent developments in the comparison to lattice QCD, at the level of the final Fermilab precision, $\Delta a_\mu^\text{exp}=16\times 10^{-11}$, also the HLbL contribution needs to be improved.
The phenomenological estimate from Ref.~\cite{Aoyama:2020ynm} (based on Refs.~\cite{Melnikov:2003xd,Colangelo:2014dfa,Colangelo:2014pva,Colangelo:2015ama,Masjuan:2017tvw,Colangelo:2017qdm,Colangelo:2017fiz,Hoferichter:2018dmo,Hoferichter:2018kwz,Gerardin:2019vio,Bijnens:2019ghy,Colangelo:2019lpu,Colangelo:2019uex,Pauk:2014rta,Danilkin:2016hnh,Jegerlehner:2017gek,Knecht:2018sci,Eichmann:2019bqf,Roig:2019reh})
\begin{equation}
\label{HLbL}
 a_\mu^\text{HLbL}=92(19)\times
10^{-11}
\end{equation}
agrees with $a_\mu^\text{HLbL}=82(35)\times 10^{-11}$ from lattice
QCD~\cite{Blum:2019ugy} (including the phenomenological estimate for the
charm contribution), and the average of both enters Eq.~\eqref{amuSM}. A
second, very recent lattice calculation~\cite{Chao:2021tvp} obtained
$a_\mu^\text{HLbL}=109.8(14.7)\times 10^{-11}$ (again after adding the
charm contribution), which agrees with both.

The modern phenomenological approach to HLbL scattering is based on dispersion relations~\cite{Hoferichter:2013ama,Colangelo:2014dfa,Colangelo:2014pva,Pauk:2014rfa,Colangelo:2015ama}, to identify in a model-independent way the contributions from hadronic intermediate states. So far, the light pseudoscalar states $\pi^0$, $\eta$, $\eta'$ have been addressed in the dispersive approach, in which case then the uncertainty simply propagates from the transition form factors (TFFs) used as input~\cite{Masjuan:2017tvw,Hoferichter:2018dmo,Hoferichter:2018kwz,Gerardin:2019vio}. Similarly, the contributions from two-pion intermediate states have been evaluated in Refs.~\cite{Colangelo:2017qdm,Colangelo:2017fiz}, including rescattering effects in the $S$-wave, which arise as unitarization of the Born-term contributions and can be interpreted as a model-independent implementation of scalar resonances, in the case of two-pion intermediate states the $f_0(500)$. Here, we extend the analysis to a coupled-channel description including $\bar K K$ intermediate states, which allows us to study in more detail the energy region of the $f_0(980)$. In particular, we can then compare the full description in terms of partial-wave helicity amplitudes for $\gamma^*\gamma^*\to\pi\pi/\bar K K$~\cite{GarciaMartin:2010cw,Hoferichter:2011wk,Moussallam:2013una,Danilkin:2018qfn,Hoferichter:2019nlq,Danilkin:2019opj} with a narrow-width approximation (NWA).

While the $\pi\pi/\bar K K$ system can be treated explicitly in terms of two-meson intermediate states,
the viability of a description in terms of narrow resonances is important for estimates of higher-multiplicity contributions, making the $f_0(980)$ an interesting test case. Such estimates will be particularly important for axial-vector intermediate states, which play a special role~\cite{Melnikov:2003xd,Jegerlehner:2017gek,Roig:2019reh,Leutgeb:2019gbz,Cappiello:2019hwh,Masjuan:2020jsf} in the transition to short-distance constraints~\cite{Melnikov:2003xd,Bijnens:2019ghy,Colangelo:2019lpu,Colangelo:2019uex,Knecht:2020xyr,Ludtke:2020moa,Bijnens:2020xnl,Bijnens:2021jqo,Colangelo:2021nkr}.
The required TFFs are subject to a set of short-distance constraints themselves~\cite{Hoferichter:2020lap}, which combined with the available experimental input~\cite{Zanke:2021wiq} should provide enough information to reduce the part of the uncertainty in Eq.~\eqref{HLbL} attributed to axial-vector intermediate states and their interplay with short-distance constraints.

However, existing estimates all rely on a simple Lagrangian definition of
such narrow-width (NW) contributions, which, in general, does not coincide
with a dispersive definition, and therefore cannot be combined with the
dispersive estimates for the one- and two-meson states. We will demonstrate
explicitly for the scalar case where the differences occur. In addition,
consistency of the dispersive approach requires a set of sum rules to be
fulfilled. For all single resonances other than pseudoscalar states these sum rules are in general not satisfied, which induces an ambiguity and renders
individual contributions dependent on the choice of the HLbL basis. Only the full result needs to fulfill the sum rules, restoring basis independence in the sum over all intermediate states.
For $S$-wave rescattering effects this potential ambiguity turns out to be
small, with the corresponding sum rule violated only at the level of $5\%$,
which allows us to provide an improved estimate of the scalar contributions
to HLbL scattering that by itself is essentially basis independent (and to be added to the pion- and kaon-box contributions). We also
comment on the role of even heavier scalar resonances and argue that those
should be included in the matching to short-distance constraints.

\section{Formalism}

\subsection{Hadronic light-by-light scattering}

We use the HLbL formalism established in Refs.~\cite{Colangelo:2017qdm,Colangelo:2017fiz}, and repeat here some of the salient features. First, following the general recipe by Bardeen, Tung~\cite{Bardeen:1969aw}, and Tarrach~\cite{Tarrach:1975tu}
(BTT), the HLbL tensor can be decomposed into $54$ Lorentz structures $T_i^{\mu\nu\lambda\sigma}$
\begin{equation}
\Pi^{\mu\nu\lambda\sigma} = \sum_{i=1}^{54} T_i^{\mu\nu\lambda\sigma} \Pi_i,
\end{equation}
with scalar functions $\Pi_i$ that encode the dynamical content of the HLbL
amplitude. These $54$ $\Pi_i$, however, form a redundant set: the number of independent structures has to match the number of helicity amplitudes, in general $41$, whereof $27$ are of relevance in the $g-2$ case of one on-shell photon. A large portion of the formalism in Refs.~\cite{Colangelo:2017qdm,Colangelo:2017fiz} is thus devoted to defining a singly-virtual basis with $27$ elements $\check\Pi_i$, in terms of which the contribution from partial-wave helicity amplitudes can be analyzed. In order to formulate the result for HLbL scattering, it is useful to return to a linear combination of the original $\Pi_i$, denoted by $\hat \Pi_i$ in Refs.~\cite{Colangelo:2017qdm,Colangelo:2017fiz}, a subset $\bar \Pi_i$ of which contribute to the master formula
\begin{align}
\label{master}
a_\mu^\text{HLbL} &= \frac{\alpha^3}{432\pi^2} \int_0^\infty d\Sigma\,\Sigma^3 \int_0^1 dr\, r\sqrt{1-r^2} \int_0^{2\pi} d\phi \notag\\
&\times\sum_{i=1}^{12} T_i(Q_1,Q_2,Q_3) \bar\Pi_i(Q_1,Q_2,Q_3),
\end{align}
where the $T_i$ are known kernel functions and the Euclidean momenta squared are given
by~\cite{Eichmann:2015nra}
\begin{align}
Q_{1}^2 &= \frac{\Sigma}{3} \left( 1 - \frac{r}{2} \cos\phi - \frac{r}{2}\sqrt{3} \sin\phi \right), \notag\\
Q_{1}^2 &= \frac{\Sigma}{3} \left( 1 - \frac{r}{2} \cos\phi + \frac{r}{2}\sqrt{3} \sin\phi \right), \notag\\
Q_3^2 &= \frac{\Sigma}{3} \left( 1 + r \cos\phi \right).
\end{align}
While the general result for higher partial waves becomes rather involved, the $S$-wave contribution can be written as~\cite{Colangelo:2017fiz}
\begin{align}
		\label{eq:SWaveDR}
		\hat\Pi_4^{J=0} &= \frac{1}{\pi} \int_{s_\text{thr}}^\infty ds' \frac{-2}{\lambda_{12}(s')(s'-q_3^2)^2} \Big( 4s' \Im h^0_{++,++}(s')\notag\\
		&- (s'+q_1^2-q_2^2)(s'-q_1^2+q_2^2) \Im h^0_{00,++}(s') \Big) , \notag\\
		\hat\Pi_{17}^{J=0} &= \frac{1}{\pi} \int_{s_\text{thr}}^\infty ds' \frac{4}{\lambda_{12}(s')(s'-q_3^2)^2} \Big( 2 \Im h^0_{++,++}(s')\notag\\
		&- (s'-q_1^2-q_2^2) \Im h^0_{00,++}(s') \Big) ,
\end{align}
where $\lambda_{12}(s)=\lambda(s,q_1^2,q_2^2)$, $\lambda(a,b,c)=a^2+b^2+c^2-2(ab+ac+bc)$ (plus crossed versions). In the special case of two-pion intermediate states, the threshold becomes $s_\text{thr}=4M_\pi^2$ and the imaginary parts for the two-pion rescattering contribution for given isospin $I=0,2$ are
\begin{align}
\label{unitarity_pipi}
	\Im h^{J,I}_{\lambda_1\lambda_2,\lambda_3\lambda_4}(s) &= \frac{\sigma_\pi(s)\theta\big(s-4M_\pi^2\big)}{32\pi} \Big[ h_{J,\lambda_1\lambda_2}^I(s) h_{J,\lambda_3\lambda_4}^{I*}(s)\notag\\
	&\qquad-c_I N_{J,\lambda_1\lambda_2}(s) N_{J,\lambda_3\lambda_4}(s)\Big],
\end{align}
subtracting the Born-term contributions
 involving $N_{J,\lambda_1\lambda_2}(s)$ and isospin factors $c_0=4/3$, $c_2=2/3$ (all amplitudes on the right-hand side refer to the $\gamma^*\gamma^*\to\pi\pi$ partial waves). The phase-space factor is $\sigma_\pi(s)=\sqrt{1-4M_\pi^2/s}$. The subtraction of the Born terms in Eq.~\eqref{unitarity_pipi} is required to avoid a double counting of the pion-box contribution.

 The formulation as in Eq.~\eqref{eq:SWaveDR} corresponds to a particular
 choice of the singly-virtual basis $\check\Pi_i$, and the requirement that
 different bases be equivalent leads to a set of sum rules that need to be
 fulfilled (ultimately, as a consequence of the $T_i^{\mu\nu\lambda\sigma}$ not all having the same mass dimension). For the $S$-waves, there is only a single combination that
 contributes to these sum rules, amounting to a relation between integrals over the two possible $S$-wave helicity projections
\begin{align}
		\label{eq:SR}
		0 &= \frac{1}{\pi} \int_{s_\text{thr}}^\infty ds' \frac{1}{\lambda_{12}(s')(s'-q_3^2)} \Big( 2 \Im h^0_{++,++}(s')\notag\\
		&- (s'-q_1^2-q_2^2) \Im h^0_{00,++}(s') \Big) +\text{higher waves}.
\end{align}
Basis independence requires that this sum rule be satisfied by the sum over all intermediate hadronic states. It is automatically fulfilled by the scalar QED amplitudes and thus the pion- and kaon-box contributions, but needs to be monitored when calculating, e.g., rescattering corrections~\cite{Colangelo:2017qdm,Colangelo:2017fiz}. For contributions that do not individually satisfy the sum rules, a basis change amounts to a reshuffling between contributions of different partial waves and/or hadronic intermediate states.

\subsection{Coupled-channel amplitudes for $\gamma^*\gamma^*\to\pi\pi/\bar K K$}

\begin{figure*}[t]
	\centering
	\includegraphics[width=0.4\linewidth,clip]{hpp}\qquad
	\includegraphics[width=0.4\linewidth,clip]{h00}\\
	\includegraphics[width=0.4\linewidth,clip]{kpp}\qquad
	\includegraphics[width=0.4\linewidth,clip]{k00}
	\caption{Modulus of the Born-term-subtracted partial waves $\bar h^0_{0,\lambda_1\lambda_2}$ and $\bar k^0_{0,\lambda_1\lambda_2}$, for a representative set of photon virtualities.}
	\label{fig:hk}
\end{figure*}

Including $\bar K K$ intermediate states, the unitarity relation~\eqref{unitarity_pipi} receives new contributions involving the full and Born-term amplitudes $k_{J,\lambda_1\lambda_2}$ and $M_{J,\lambda_1,\lambda_2}$, respectively. Note that these partial-wave amplitudes are normalized in such a way as to ensure the same unitarity condition for identical and non-identical particles~\cite{GarciaMartin:2010cw,Danilkin:2019opj}. Concentrating on the rescattering contribution for $I=J=0$, we have
\begin{align}
\label{unitarity_pipiKK}
	\Im h^{0,0}_{\lambda_1\lambda_2,\lambda_3\lambda_4}(s) &= \frac{\sigma_\pi(s)\theta\big(s-4M_\pi^2\big)}{32\pi} \Big[ h_{0,\lambda_1\lambda_2}^0(s) h_{0,\lambda_3\lambda_4}^{0*}(s)\notag\\
	&\qquad-\frac{4}{3} N_{0,\lambda_1\lambda_2}(s) N_{0,\lambda_3\lambda_4}(s)\Big]\notag\\
	&+\frac{\sigma_K(s)\theta\big(s-4M_K^2\big)}{32\pi} \Big[ k_{0,\lambda_1\lambda_2}^0(s) k_{0,\lambda_3\lambda_4}^{0*}(s)\notag\\
	&\qquad-\frac{1}{2} M_{0,\lambda_1\lambda_2}(s) M_{0,\lambda_3\lambda_4}(s)\Big],
\end{align}
where the $\gamma^*\gamma^* \to \pi\pi$ and $\gamma^*\gamma^* \to K\bar{K}$ amplitudes are obtained using a modified coupled-channel Muskhelishvili--Omn\`es (MO) formalism~\cite{Danilkin:2019opj} (whose notation we follow). Besides eliminating all kinematic constraints in the partial-wave helicity amplitudes~\cite{Colangelo:2017fiz}, the MO formalism requires as input the knowledge of the left-hand cuts and the  hadronic Omn\`es functions~\cite{Omnes:1958hv}. The latter we take from a data-driven $N/D$ analysis~\cite{Danilkin:2020pak}, in which the fit is performed to the most recent Roy and Roy--Steiner results on $\pi\pi \to \pi\pi$~\cite{GarciaMartin:2011cn} and $\pi\pi \to \bar K K$~\cite{Pelaez:2020gnd}, respectively. By analytic continuation to the complex plane, this solution produces $f_0(500)$ and $f_0(980)$ poles at $\sqrt{s_{f_0(500)}}=458(10)^{+7}_{-15} - i\, 256(9)^{+5}_{-8}\,\text{MeV}$ and $\sqrt{s_{f_0(980)}}=993(2)^{+2}_{-1} - i\,21(3)^{+2}_{-4}\,\text{MeV}$, in good agreement with Refs.~\cite{Caprini:2005zr,GarciaMartin:2011jx,Moussallam:2011zg} (see Ref.~\cite{Danilkin:2020pak} for more details on the uncertainty estimates). Since the constructed $N/D$ solution is based on a once-subtracted dispersion relation, the obtained Omn\`es matrix is bounded  asymptotically. As for the left-hand cuts, it has been verified by comparison to the on-shell data from Refs.~\cite{Marsiske:1990hx,Boyer:1990vu,Behrend:1992hy,Mori:2007bu,Uehara:2008ep,Uehara:2009cka} that the MO formalism based on the pion- and kaon-pole left-hand cuts alone provides a good description of the $f_0(500)$ and $f_0(980)$ regions; see, e.g., Ref.~\cite{Aoyama:2020ynm}.\footnote{This statement no longer holds true for the $D$-waves, see Refs.~\cite{GarciaMartin:2010cw,Danilkin:2018qfn,Hoferichter:2019nlq,Danilkin:2019opj}, for which vector-meson left-hand cuts need to be included, as determined by the respective TFFs~\cite{Niecknig:2012sj,
Hoferichter:2012pm,Schneider:2012ez,Danilkin:2014cra,Hoferichter:2014vra,Hoferichter:2017ftn,Albaladejo:2020smb}.} In particular, the two-photon widths $\Gamma_{\gamma\gamma}[f_0(500)]=1.37(13)^{+0.09}_{-0.06}\,\text{keV}$ and $\Gamma_{\gamma\gamma}[f_0(980)]=0.33(16)^{+0.04}_{-0.16}\,\text{keV}$~\cite{Danilkin:2020pak}, come out consistent with other dispersive  extractions~\cite{ Hoferichter:2011wk, Moussallam:2011zg, Dai:2014zta}. For the pion and kaon electromagnetic form factors that enter in the pion- and kaon-pole contributions for virtual photons, we use the vector-meson-dominance (VMD)  expressions
\begin{align}
 F_\pi^V(s)&=\frac{M_\rho^2}{M_\rho^2-s},\notag\\
 F_K^V(s)&=\frac{1}{2}\frac{M_\rho^2}{M_\rho^2-s}
 +\frac{1}{6} \frac{M_\omega^2}{M_\omega^2-s}
 +\frac{1}{3} \frac{M_\phi^2}{M_\phi^2-s},
\end{align}
given that the difference to the full description is negligible~\cite{Colangelo:2017fiz}. The resulting rescattering contributions to the pion and kaon partial waves are shown in Fig.~\ref{fig:hk} for a representative set of photon virtualities. Both the $f_0(500)$ and the $f_0(980)$ are clearly visible in $\bar h^0_{0,\lambda_1\lambda_2}$, while the impact of the $f_0(980)$ is also reflected by the threshold enhancement in $\bar k^0_{0,\lambda_1\lambda_2}$.

\subsection{Narrow-width approximation}
\label{sec:NWA}

For a narrow scalar resonance with mass $m_S$ we decompose the matrix element with two electromagnetic currents $j_\mathrm{em}^\mu$ according to~\cite{Hoferichter:2020lap}
\begin{align}
\mathcal{M}^{\mu\nu}(p \rightarrow q_1,q_2) &= i \int d^4x \, e^{i q_1 \cdot x} \< 0 | T \{ j_\mathrm{em}^\mu(x) j_\mathrm{em}^\nu(0) \} | S(p) \> \notag\\
&=\frac{\mathcal{F}_1^S(q_1^2,q_2^2)}{m_S} T_1^{\mu\nu}  + \frac{\mathcal{F}_2^S(q_1^2,q_2^2)}{m_S^3} T_2^{\mu\nu} ,
\end{align}
with Lorentz structures
\begin{align}
\label{BTT_scalar}
		T_1^{\mu\nu} &= q_1 \cdot q_2 g^{\mu\nu} - q_2^\mu q_1^\nu , \notag\\
		T_2^{\mu\nu} &= q_1^2 q_2^2 g^{\mu\nu} + q_1 \cdot q_2 q_1^\mu q_2^\nu - q_1^2 q_2^\mu q_2^\nu - q_2^2 q_1^\mu q_1^\nu,
\end{align}
and TFFs $\mathcal{F}_1^S$, $\mathcal{F}_2^S$.
This BTT decomposition is again free of kinematic singularities, and the normalization of $\mathcal{F}_1^S$ can be related to the $S\to\gamma\gamma$ partial width $\Gamma_{\gamma\gamma}$:
\begin{equation}
	| \mathcal{F}_1^S(0,0) |^2 = \frac{4}{\pi \alpha^2 m_S} \Gamma_{\gamma\gamma} .
\end{equation}
The contribution to HLbL scattering follows most easily by the replacement
\begin{align}
 \Im h^0_{++,++}(s)&=\bigg(\frac{m_S^2-q_1^2-q_2^2}{2m_S}\mathcal{F}_1^S(q_1^2,q_2^2)+\frac{q_1^2q_2^2}{m_S^3}\mathcal{F}_2^S(q_1^2,q_2^2)\bigg)\notag\\
 &\times\frac{m_S^2-q_3^2}{2m_S}\mathcal{F}_1^S(q_3^2,0)\times \pi\delta(s-m_S^2),\notag\\
 \Im h^0_{00,++}(s)&=\bigg(\frac{1}{m_S}\mathcal{F}_1^S(q_1^2,q_2^2)+\frac{m_S^2-q_1^2-q_2^2}{2m_S^3}\mathcal{F}_2^S(q_1^2,q_2^2)\bigg)\notag\\
 &\times\frac{m_S^2-q_3^2}{2m_S}\mathcal{F}_1^S(q_3^2,0)\times \pi\delta(s-m_S^2),
\end{align}
leading to
\begin{align}
	\label{eq:hatFunctionsNarrowScalar}
 \hat\Pi_4&=\frac{\mathcal{F}_1^S(q_3^2,0)}{q_3^2-m_S^2}\bigg(\frac{\mathcal{F}_1^S(q_1^2,q_2^2)}{m_S^2}-\frac{m_S^2+q_1^2+q_2^2}{2m_S^4}\mathcal{F}_2^S(q_1^2,q_2^2)\bigg),\notag\\
 \hat\Pi_{17}&=\frac{\mathcal{F}_1^S(q_3^2,0)}{q_3^2-m_S^2} \frac{\mathcal{F}_2^S(q_1^2,q_2^2)}{m_S^4},
\end{align}
together with crossed versions for $t$- and $u$-channel exchange. We stress that in a Lagrangian model as formulated in Refs.~\cite{Pauk:2014rta,Knecht:2018sci}, the numerator of the prefactor of $\mathcal{F}_2^S$ in the bracket is altered to $m_S^2 + q_1^2 + q_2^2 \mapsto q_3^2 + q_1^2 + q_2^2$, i.e., while the residues of the scalar pole agree, the dispersive and the model description differ by non-pole terms. This difference could be removed by a further change in the HLbL basis, but then of course all other contributions to HLbL scattering would also need to be evaluated in this new basis, including the comparison to the rescattering corrections.

In the NWA, the sum rule~\eqref{eq:SR} evaluates to
\begin{equation}
0=-\frac{1}{4m_S^4}\mathcal{F}_2^S(q_1^2,q_2^2) \mathcal{F}_1^S(q_3^2,0) + \text{other states},
\end{equation}
so that, unless $\mathcal{F}_2^S=0$, the contribution of a single resonance is not unique because the sum rules that ensure basis independence cannot be fulfilled by a narrow scalar alone. This statement applies to all resonances apart from pseudoscalar states, which are not affected by the sum rules. While a complete evaluation of HLbL must not depend on the choice of basis and fulfill the sum rules, this is not the case for individual intermediate states. In consequence, NW estimates for HLbL contributions necessitate the specification of the chosen HLbL basis, as a basis change amounts to a reshuffling of contributions from different intermediate states. Here, we take the basis from Refs.~\cite{Colangelo:2017qdm,Colangelo:2017fiz}, to contrast a NW description with one based on the $\gamma^*\gamma^*\to\pi\pi/\bar K K$ partial waves.

As follows from Eq.~\eqref{BTT_scalar}, the second TFF $\mathcal{F}_2^S$ only contributes to doubly-virtual processes, so that no direct information from experiment exists. However,
for large virtualities the TFFs can be analyzed in a light-cone expansion, whose leading result gives~\cite{Hoferichter:2020lap}
\begin{align}
 \mathcal{F}_1^S(q_1^2,q_2^2)&=F_S^\text{eff}m_S\int_0^1 d u\,\frac{3(2u-1)^2\phi(u)}{u q_1^2+(1-u)q_2^2},\notag\\
 \mathcal{F}_2^S(q_1^2,q_2^2)
 	&= F_S^\text{eff} m_S^3 \int_0^1 d u\, \frac{3u(1-u)\phi(u)}{( u q_1^2+(1-u)q_2^2)^2},
 \label{result_scalar}
\end{align}
with distribution amplitude $\phi(u)$ and an effective decay constant $F_S^\text{eff}$.
Inserting the asymptotic form $\phi(u)=6u(1-u)$, the two integrals become related, leading to the expansion
\begin{equation}
\label{LC_scalar}
 \mathcal{F}_1^S(q_1^2,q_2^2)=\frac{F_S^\text{eff} m_S}{Q^2}f^S(w),\qquad
 \mathcal{F}_2^S(q_1^2,q_2^2)=\frac{F_S^\text{eff} m_S^3}{Q^4}f^S(w),
 \end{equation}
with
\begin{equation}
\label{defQ}
Q^2=\frac{q_1^2+q_2^2}{2},\qquad w=\frac{q_1^2-q_2^2}{q_1^2+q_2^2},
\end{equation}
and
\begin{equation}
f^S(w)=\frac{3}{2w^4}\bigg(3-2w^2+3\frac{1-w^2}{2w}\log\frac{1-w}{1+w}\bigg).
\end{equation}
These results can be contrasted with the quark model from Ref.~\cite{Schuler:1997yw}
\begin{align}
\label{TFF_Schuler}
  \frac{\mathcal{F}_1^S(q_1^2,q_2^2)}{\mathcal{F}_1^S(0,0)}\bigg|_{\text{\cite{Schuler:1997yw}}}&=\frac{m_S^2(3m_S^2-q_1^2-q_2^2)}{3(m_S^2-q_1^2-q_2^2)^2}= -\frac{m_S^2}{6Q^2}+\mathcal{O}\Big(Q^{-4}\Big),\notag\\
 \frac{\mathcal{F}_2^S(q_1^2,q_2^2)}{\mathcal{F}_1^S(0,0)}\bigg|_{\text{\cite{Schuler:1997yw}}}&=-\frac{2m_S^4}{3(m_S^2-q_1^2-q_2^2)^2}=-\frac{m_S^4}{6Q^4}+\mathcal{O}\Big(Q^{-6}\Big),
\end{align}
which thus reproduces the asymptotic relation $\mathcal{F}_2^S(q_1^2,q_2^2)/\mathcal{F}_1^S(q_1^2,q_2^2)=m_S^2/Q^2$ from Eq.~\eqref{LC_scalar}, while correctly interpolating to the normalization for $\mathcal{F}_1^S$ (to also implement the correct $w$ dependence more complicated TFF parameterizations would be necessary, see, e.g., Refs.~\cite{Hoferichter:2018kwz,Zanke:2021wiq}). We will use the prescription~\eqref{TFF_Schuler} to evaluate our NW estimates in Sec.~\ref{sec:NW}.

\section{Two-meson rescattering}

In Refs.~\cite{Colangelo:2017qdm,Colangelo:2017fiz} the effect of the $f_0(500)$ was estimated in terms of $S$-wave $\pi\pi$ rescattering to $a_\mu^\text{HLbL}[f_0(500)]=-9(1)\times 10^{-11}$, where the uncertainty was mainly attributed to the high-energy continuation of the amplitudes and, in practice, assessed by the violation of the sum rule that needs to be fulfilled to make the contribution basis independent. In particular, the $\pi\pi$ rescattering was implemented using the $\pi\pi$ phase shift from the inverse-amplitude method  in the Omn\`es solution for $\gamma^*\gamma^*\to\pi\pi$, to explicitly remove the effects of the $f_0(980)$ and ensure a smooth high-energy behavior.

\begin{table}[t]
	\centering
	\begin{tabular}{l c r  r}
	\toprule
	$\Lambda\ [\text{GeV}]$&  & $0.89$ & $2.0$\\\midrule
	\multicolumn{2}{l}{pion (+ kaon) Born terms ($S$-waves)} & $-11.4$ & $-11.8$\\
	\multicolumn{2}{l}{$S$-wave $I=0$ rescattering} & $-10.0$ & $-9.8$\\\midrule
	\multirow{3}{4cm}{sum rule pion (+ kaon) \\ Born terms ($S$-waves)} & $++,++$ & $8.0$ & $8.4$\\
	& $00,++$ & $-9.2$ & $-9.6$\\
	& total & $-1.2$ & $-1.2$\\\midrule
	\multirow{3}{4cm}{sum rule\\ $S$-wave  $I=0$ rescattering} & $++,++$ & $6.9$ & $6.8$\\
	 & $00,++$ & $-7.3$& $-7.2$\\
	& total & $-0.4$ & $-0.4$\\
\bottomrule
	\end{tabular}
	\caption{$S$-wave Born-term and $I=0$ rescattering contributions (upper panel), and helicity components of the sum rule (lower panels), all in units of $10^{-11}$. $\Lambda^2$ denotes the cutoff in the $s$ integration.}
	\label{tab:rescattering}
\end{table}

For the coupled-channel solution from Ref.~\cite{Danilkin:2019opj} we find the results shown in Table~\ref{tab:rescattering}, very close to the estimate
from Refs.~\cite{Colangelo:2017qdm,Colangelo:2017fiz}.\footnote{The Born terms correspond to the contributions that are subtracted in Eq.~\eqref{unitarity_pipiKK}. Resumming all Born-term partial waves would result in $\frac{3}{2}$ times the corresponding box contributions due to the different double-spectral regions of box and rescattering contributions~\cite{Colangelo:2017fiz}.} The numerical result is completely dominated by the two-pion contribution, but we verified that the imaginary part cancels between the $\pi\pi$ and $\bar K K$ channels, as required by the unitarity relation. We also checked the contribution to the sum rule~\eqref{eq:SR}, estimated by the difference when evaluated with the alternative basis from
Ref.~\cite{Colangelo:2014dfa}. As in Refs.~\cite{Colangelo:2017qdm,Colangelo:2017fiz} we observe a cancellation of the helicity components up to a remainder of $5\%$, indicating that sum-rule violations either due to the high-energy region in the $S$-waves or higher partial waves are indeed small. In fact, for the $S$-wave contributions of the pion and kaon Born terms the sum-rule violations are at least twice this relative size, but of course in this case the cancellation via the higher partial waves is exact. This implies that the $S$-waves can essentially be regarded as basis independent and the residual sum-rule violation be treated as an uncertainty estimate for the high-energy continuation. In combination with the $I=2$ estimates from Refs.~\cite{Colangelo:2017qdm,Colangelo:2017fiz}, $a_\mu^\text{HLbL}[\text{$S$-wave, $I=2$, $\Lambda=2\,\text{GeV}$}]=1.1\times 10^{-11}$, we quote for the complete $S$-wave rescattering
\begin{equation}
\label{rescattering_final}
a_\mu^\text{HLbL}[\text{$S$-wave rescattering}]=-8.7(1.0)\times 10^{-11},
\end{equation}
where the uncertainty covers the high-energy region, sum-rule violations, and input quantities.

One could then define an $f_0(980)$ contribution by integrating
over a window in $\sqrt{s}$, an obvious choice being
$\sqrt{s}\in[M_{f_0(980)}-\Gamma_{f_0(980)},M_{f_0(980)}+\Gamma_{f_0(980)}]$. The parameters quoted in Ref.~\cite{Zyla:2020zbs} are
\begin{equation}
M_{f_0(980)}=990(20)\,\text{MeV},\qquad \Gamma_{f_0(980)}=(10\text{--}100)\,\text{MeV},
\end{equation}
while dispersive analyses~\cite{GarciaMartin:2011jx,Moussallam:2011zg} favor a width around $50(20)\,\text{MeV}$, which already reflects the complicated line shape of the $f_0(980)$ due to the close proximity of the $\bar K K$ threshold~\cite{Flatte:1976xu,Baru:2003qq}.
We therefore suggest a different definition based on the decomposition
\begin{equation}
\label{amu_int}
a_\mu^\text{HLbL}=\int_{4M_\pi^2}^\infty ds' a_\mu^\text{HLbL}(s'),
\end{equation}
see Fig.~\ref{fig:amus} for the integrand produced by the rescattering corrections. This shows again that the $f_0(500)$ region by far dominates, while concentrating on the energy region around the $f_0(980)$ allows one to visualize the resulting line shape. Integrating the deficit below the baseline around $a_\mu^\text{HLbL}(s')=0.3\times 10^{-11}\,\text{GeV}^{-2}$ gives $-0.18\times 10^{-11}$, which suggests an $f_0(980)$
contribution around
\begin{equation}
\label{rescatt}
a_\mu^{\text{HLbL}}[f_0(980)]\big|_\text{rescattering}= -0.2(1)\times 10^{-11}.
\end{equation}
The strong suppression compared to Eq.~\eqref{rescattering_final} happens to be of similar size as the suppression of the kaon to the pion box, $a_\mu^\text{HLbL}[\text{$\pi$-box}]/a_\mu^\text{HLbL}[\text{$K$-box}]\sim 32$~\cite{Aoyama:2020ynm}, which in this case is partly driven by $M_K^2/M_\pi^2$, but given the complicated dependence of the kernel functions on the virtualities and the different line shapes the general mass scaling is difficult to anticipate.
We stress that a contribution such as Eq.~\eqref{rescatt}, defined via a particular window in $\sqrt{s}$, is in general not basis independent, so that the comparison to the NWA described in the following section should refer to the same HLbL basis.

\begin{figure*}[t]
	\centering
	\includegraphics[width=0.4\linewidth,clip]{amu_f0500}\qquad
	\includegraphics[width=0.4\linewidth,clip]{amu_f0980}
	\caption{Integrand from Eq.~\eqref{amu_int} for the rescattering contribution, for the $f_0(500)$ (left) and the $f_0(980)$ (right). Note the different scales in both cases.}
	\label{fig:amus}
\end{figure*}

\section{Narrow-width estimates}
\label{sec:NW}

The result from the rescattering approach can be contrasted with a description in terms of a narrow resonance. We use a dispersive definition in line with the HLbL basis chosen for the rescattering, and use the TFFs from the quark model of Ref.~\cite{Schuler:1997yw}, which provides a plausible interpolation between the normalization and the short-distance constraints derived in Ref.~\cite{Hoferichter:2020lap}. Using $M_{f_0(980)}=0.99\,\text{GeV}$ and $\Gamma_{\gamma\gamma}[f_0(980)]=0.31(5)\,\text{keV}$~\cite{Zyla:2020zbs}, we find
\begin{equation}
\label{NW}
a_\mu^{\text{HLbL}}[f_0(980)]\big|_\text{NWA}=-0.37(6)\times 10^{-11},
\end{equation}
not too far away from the rescattering definition Eq.~\eqref{rescatt} (and an uncertainty referring only to $\Gamma_{\gamma\gamma}$). If the scale in the TFF parameterization~\eqref{TFF_Schuler} were lowered to a VMD expectation, $m_S\to M_\rho$, the NW result would move to $-0.26(4)\times 10^{-11}$, even closer to Eq.~\eqref{rescatt}.
Taken together with Eq.~\eqref{rescatt}, this would suggest the estimate
\begin{equation}
a_\mu^{\text{HLbL}}[f_0(980)]=-0.2(2)\times 10^{-11},
\end{equation}
and with $M_{a_0(980)}=0.98\,\text{GeV}$ and $\Gamma_{\gamma\gamma}[a_0(980)]=0.3(1)\,\text{keV}$~\cite{Zyla:2020zbs} a similar range would be expected for the $a_0(980)$. An improved evaluation of the isospin $I=1$ channel could be obtained from a coupled-channel analysis of the doubly-virtual helicity amplitudes for $\gamma^*\gamma^*\to\pi\eta/\bar KK$~\cite{Danilkin:2017lyn,Deineka:2018nuh,Lu:2020qeo}, following the same strategy as for the $f_0(980)$.
In fact, the coupled-channel analysis from Ref.~\cite{Lu:2020qeo} (based on the data from Refs.~\cite{Albrecht:1989re,Uehara:2009cf,Uehara:2013mbo}) prefers a width $\Gamma_{\gamma\gamma}[a_0(980)]=0.5^{+0.2}_{-0.1}\,\text{keV}$, which translates to
\begin{equation}
\label{a0_estimate}
a_\mu^{\text{HLbL}}[a_0(980)]=-\Big(0.6_{-0.1}^{+0.2}\Big)\times 10^{-11},
\end{equation}
or a slightly lower range, $a_\mu^{\text{HLbL}}[a_0(980)]=-\big(0.4_{-0.1}^{+0.2}\big)\times 10^{-11}$, if instead a VMD scale were used in the TFFs.

We stress that when combining or comparing different contributions to HLbL, one should work within one unified framework, provided here by the dispersive framework and tensor basis of Refs.~\cite{Colangelo:2017qdm,Colangelo:2017fiz}. The basis dependence of the narrow scalar contribution can be illustrated by either setting $\mathcal{F}_2^S = 0$, or using the Lagrangian description by the replacement $m_S^2 + q_1^2 + q_2^2 \mapsto q_3^2 + q_1^2 + q_2^2$ as mentioned after Eq.~\eqref{eq:hatFunctionsNarrowScalar}, which with otherwise unchanged input for the TFFs (in particular, scale $m_S$) would result in
\begin{align}
	\label{eq:NoF2Lagrangian}
	a_\mu^{\text{HLbL}}[f_0(980)]\big|_\text{NWA}^{\mathcal{F}_2^S=0}=-0.47(8)\times 10^{-11} , \nonumber\\
	a_\mu^{\text{HLbL}}[f_0(980)]\big|_\text{NWA}^\text{Lagrangian}=-0.25(4)\times 10^{-11} .
\end{align}

In Ref.~\cite{Pauk:2014rta}, the $f_0(980)$ contribution is estimated in a Lagrangian model, keeping only the transverse helicity amplitude, which is then parameterized using a monopole form factor with scale varied between $(1\text{--}2)\,\text{GeV}$, leading to a range $-(0.19\text{--}0.61)\times
10^{-11}$. The representation in terms of a single helicity amplitude combined with the Lagrangian definition resulted in kinematic singularities, which were removed by hand through angular averages.
We emphasize that we cannot use the same input as Ref.~\cite{Pauk:2014rta} to reproduce these results using the BTT master formula for HLbL, in which a priori there are no kinematic singularities. A transverse $f_0(980)$ is obtained for $\mathcal{F}_2 = -2m_S^2/(m_S^2-q_1^2-q_2^2)\mathcal{F}_1$~\cite{Hoferichter:2020lap}, which with a monopole ansatz for $\mathcal{F}_1$ gives for the pole contribution a range $-(0.40\text{--}0.49)\times 10^{-11}$. Multiplying both form factors with an additional factor of $(m_S^2 - q_1^2 - q_2^2)/\lambda_{12}(m_S^2)$ would be closer in spirit to Ref.~\cite{Pauk:2014rta}, leading to a range $-(0.55\text{--}0.97)\times 10^{-11}$. Both variants are quite different from the range from Ref.~\cite{Pauk:2014rta} quoted above.

A NWA for the $f_0(980)$ is also considered in Ref.~\cite{Knecht:2018sci}, which uses the tensor decomposition~\eqref{BTT_scalar} without kinematic singularities, but again is based on a Lagrangian definition of the scalar contribution. The results are given as a function of a parameter $\kappa_S$, where $\kappa_S=0$ corresponds to switching off $\mathcal{F}_2^S$ and reduces $\mathcal{F}_1^S$ to a pure VMD form. The result without $\mathcal{F}_2^S$, $-0.42(9)\times 10^{-11}$, is close to Eq.~\eqref{eq:NoF2Lagrangian}, and quantifies the difference in the TFF input, where we believe that the quark model from Ref.~\cite{Schuler:1997yw} is more reliable because of the better implementation of the short-distance constraints~\cite{Hoferichter:2020lap} described in Sec.~\ref{sec:NWA} (see Refs.~\cite{Hoferichter:2020lap,Kroll:2016mbt} for the comparison to the singly-virtual data from Belle~\cite{Masuda:2015yoh}).
The difference to the results including $\mathcal{F}_2^S$, $-0.07(2)\times10^{-11}$, is mainly explained by their Lagrangian definition that includes non-pole pieces and to a lesser degree by the different TFF parameterizations. We checked that their spread for $\kappa_S \in [0,1]$ is much reduced when the dispersive basis of Refs.~\cite{Colangelo:2017qdm,Colangelo:2017fiz} is used instead.

References~\cite{Pauk:2014rta,Knecht:2018sci} also consider even heavier scalars, based on estimates of their two-photon coupling, e.g., $\Gamma_{\gamma\gamma}[f_0(1370)]=3.8(1.5)\,\text{keV}$~\cite{Boglione:1998rw}. This estimate, however, describes a combined effect of $f_0(500)$ and $f_0(1370)$, which could not be reliably separated at the time. In more recent partial-wave analyses the $f_0(500)$ couplings can be isolated, while the effect of the $f_0(1370)$ is barely seen in $\gamma\gamma$ reactions. In fact, the number $\Gamma_{\gamma\gamma}[f_0(1370)]=4.0(1.9)\,\text{keV}$ given in Ref.~\cite{Dai:2014zta} is accompanied by an explicit warning that even with its large error this number does not have the credibility of the other two-photon couplings (associating one star with the result).
The situation is slightly better for the $a_0(1450)$, for which Ref.~\cite{Lu:2020qeo} quotes $\Gamma_{\gamma\gamma}[a_0(1450)]=1.05^{+0.50}_{-0.30}\,\text{keV}$. Using U(3) assumptions, the decay widths of the excited scalars are related by
\begin{align}
\Gamma_{\gamma\gamma}[a_0(1450)]
&=\frac{\Gamma_{\gamma\gamma}[f_0(1370)]}{3\cos^2(\theta_A-\theta_0)}\frac{M_{a_0(1450)}}{M_{f_0(1370)}}\notag\\
&=\frac{\Gamma_{\gamma\gamma}[f_0(1500)]}{3\sin^2(\theta_A-\theta_0)}\frac{M_{a_0(1450)}}{M_{f_0(1500)}},
\end{align}
where $\theta_0=\arcsin (1/3)$ and $\theta_A$ is the mixing angle between $f_0(1370)$ and $f_0(1500)$~\cite{Zyla:2020zbs} (octet/singlet mixing is reproduced for $\theta_A=\pi/2$, see, e.g., Ref.~\cite{Zanke:2021wiq}). Since the $f_0(1500)$ has not been seen in $\gamma\gamma$ collisions~\cite{Zyla:2020zbs,Barate:1999ze,Acciarri:2000ex,Uehara:2008ep}, one could determine $\theta_A$ by the requirement that $\Gamma_{\gamma\gamma}[f_0(1500)]=0$. This choice, $\theta_A=\theta_0$, essentially defines an upper limit for $\Gamma_{\gamma\gamma}[f_0(1370)]\leq 3.0^{+1.4}_{-0.9}\,\text{keV}$, where the uncertainties are propagated from Ref.~\cite{Lu:2020qeo} but do not include an additional U(3) uncertainty. Translating these couplings into a contribution to HLbL scattering is further complicated by the absence of information on the corresponding TFFs, and for the heavy scalars it makes a bigger difference if the scale is set by $m_S$ or VMD. The corresponding results are
\begin{align}
 a_\mu^{\text{HLbL}}[f_0(1370)]&=-(1.5^{+0.7}_{-0.4})\times 10^{-11}\quad \big[-(0.6^{+0.3}_{-0.2})\times 10^{-11}\big],\notag\\
 a_\mu^{\text{HLbL}}[a_0(1450)]&=-(0.5^{+0.2}_{-0.1})\times 10^{-11}\quad \big[-(0.2^{+0.1}_{-0.05})\times 10^{-11}\big],
\end{align}
where for the $f_0(1370)$ we adopted the above U(3) estimate and the
numbers in brackets are obtained for $m_S\to M_\rho$ in
Eq.~\eqref{TFF_Schuler}. The comparison for the $f_0(980)$ to the
implementation in terms of partial waves suggests that the latter results
may be more reliable, pointing to a combined effect of
$f_0(1370)$ and $a_0(1450)$ of at most $-1\times 10^{-11}$.

However, given that even the two-photon couplings of the heavy scalar resonances are
highly uncertain, let alone their TFFs, we propose a different point of view here. In comparison to the tensor mesons $f_2(1270)$ and $a_2(1320)$ these states are not seen prominently in two-photon reactions, and in contrast to axial-vector states they are not expected to play a special role in the implementation of short-distance constraints. Moreover, with a mass around $1.5\,\text{GeV}$ and in view of the substantial uncertainties, it is not clear that a description in terms of hadronic degrees of freedom is useful, while it should be more promising to cover the respective physics in the asymptotic matching~\cite{Bijnens:2019ghy,Bijnens:2020xnl,Bijnens:2021jqo}.
Comparing to the additional scalar contribution given in Ref.~\cite{Aoyama:2020ynm}, we have now included the $f_0(980)$ in the rescattering part~\eqref{rescattering_final}, and combined with the $a_0(980)$ NW estimate from Eq.~\eqref{a0_estimate} we quote
\begin{equation}
\label{scalars_final}
a_\mu^\text{HLbL}[\text{scalars}]=-9(1)\times 10^{-11}
\end{equation}
as our final result for the $S$-wave contribution to HLbL scattering to be added to the pion and kaon boxes.
With the effects of even heavier states moved into the asymptotic matching, this eliminates the need for an additional scalar contribution.

\section{Summary and outlook}

In this paper we addressed scalar contributions to HLbL scattering in the framework of dispersion relations. First, we extended previous work on the $f_0(500)$ resonance, implemented via rescattering corrections to two-pion intermediate states, to the coupled-channel system of $\pi\pi/\bar K K$, which allowed us to identify a contribution from the $f_0(980)$ resonance. The result is reasonably close to an estimate using narrow resonances, given the uncertainties inherent in the scalar transition form factors. With a similar estimate for the $a_0(980)$ we find that the combined effect
$a_\mu^{\text{HLbL}}[f_0(980)+a_0(980)]$ is well below $1\times 10^{-11}$,
so that our final result for the scalar contributions~\eqref{scalars_final}
is by far dominated by the $f_0(500)$ region. As for the contribution of
even heavier resonances, we argued that given scant
experimental input for their two-photon couplings and transition form
factors a description in terms of hadronic degrees of freedom is not
particularly useful and that their contribution should be included in the
asymptotic matching.

Another important goal of this paper is to emphasize conceptual issues
that first arise for the scalar contributions, but will become more
critical for axial-vector and tensor resonances. First of all, it is only
in a dispersive framework that the narrow-width approximation corresponds
to including a pure pole term. Definitions based on phenomenological
Lagrangians usually include non-pole terms, which are model dependent,
modify the high-energy behavior, and have a significant numerical impact.
This observation is critical, because to ensure consistency of the entire
HLbL result, each contribution needs to be defined and evaluated
within the same framework. Even within the dispersive approach a basis
change for the HLbL tensor leads to an ambiguity in the definition of
each individual contribution, because a set of sum rules that guarantees basis independence
in general only needs to be satisfied by the sum over all hadronic intermediate states. We
demonstrated that for the $S$-wave rescattering the only relevant sum rule
is well fulfilled thanks to a cancellation between helicity components, so
that the result~\eqref{rescattering_final} is essentially basis
independent. For future estimates of axial-vector and tensor contributions
these subtleties will require a careful treatment, and together with the
matching to short-distance constraints will be critical to improve the
precision in the evaluation of HLbL scattering.

\section*{Acknowledgments}

We thank G.~Colangelo for very fruitful discussions and collaboration on parts of this work.
Financial support by the DOE (Grant Nos.\
DE-FG02-00ER41132 and DE-SC0009919) and the SNSF (Project No.\ PCEFP2\_181117) is gratefully acknowledged.  The work of I.D.\ was supported by the Deutsche Forschungsgemeinschaft (German Research Foundation) in part through the Collaborative Research Center (The Low-Energy Frontier of the Standard Model, Projektnummer 204404729 -- SFB 1044) and in part through the Cluster of Excellence (Precision Physics, Fundamental Interactions, and Structure of Matter, PRISMA$^+$ EXC 2118/1) within the German Excellence Strategy (Project No.\ 39083149).
We thank the Institute for Nuclear Theory at the University of Washington for its kind hospitality and  the DOE for partial support during INT Workshop INT-19-74W, where this project was initiated.

\bibliographystyle{apsrev4-1_mod}
\balance
\biboptions{sort&compress}


\subsection{Operation Procedure}
\begin{figure}[tbp]
    \centering
    \includegraphics[width=.75\columnwidth]{figs/flowchart.pdf}
    \caption{Real-time operation procedure}
    \label{fig:flowchart}
\end{figure}

This section introduces the real-time operation procedure for the DSO to manage the grid. PV capacity, EV power and energy capacity, HP capacity, water tank volume and temperature range are assumed known parameters to the DSO upon their installation. As shown in Fig. \ref{fig:flowchart}, before each operation step, the DSO collects power profiles from the customers via ICT infrastructures for the next time step, and runs a power flow program to evaluate the grid states. In case bus voltage limits or cable loading limits are violated, the DSO again collects data for model parameters and activates the coordination model to optimize the active and reactive power setpoints of various controllable DERs. It is assumed water tanks are used for domestic hot water heating and space heating. Their temperature is thus used as an important parameter to determine the urgentness of their power demand. Correspondingly, EVs' state of charge (SoC) and departure time are used to determine their urgentness. Finally, the optimized active and reactive power setpoints are re-fed into the power flow program for evaluation and grid state extraction considering that the proposed coordination model adopts relaxation for the power flow relations. Control actions are implemented afterwards.

\subsection{DER Coordination Model}
\begin{table}[tbp]
\begin{mdframed}
\input{nomenclature}
\printnomenclature
\end{mdframed}
\end{table}

This section presents the convex DER coordination model, which is implemented when grid limits are violated. Notations used are listed in the nomenclature. The model aims to restore grid limits by exploiting flexibility from various DERs, in which the decision variables include active and reactive power setpoints of EVs, HPs, and PVs connected to the grid. The objective function (\ref{obj}) minimizes the deviation from the DER users' original schedules collected from the customers, considering meanwhile losses in the grid. The cost terms are determined such that DER users with more urgent demand are assigned with higher values, which lowers the chance of their demand being curtailed. In this way, user satisfaction is more likely to be met.

The constraints in the model are presented in the follows. The active power and reactive power injections into the grid from each household are defined with (\ref{pbalance}) and (\ref{qbalance}), respectively. The PVs' apparent power ratings are enforced in (\ref{pv1}), while (\ref{pv2}) reflects PV's power factor limits. Similar constraints for EVs are defined in (\ref{ev1}) and (\ref{ev2}). Fixed power factors for HPs are enforced in (\ref{hp}). Constraints (\ref{evcurtail}) and (\ref{pvcurtail}) respectively define active power curtailment for EVs and PVs. Constraints (\ref{hpcurtail1})-(\ref{hpcurtail3}) combined define the downward and upward regulation of active power for HPs.

The power flow relations are defined in (\ref{pflow})-(\ref{ilimits}), which adopts a second-order cone formulation \cite{Farivar2013}. Specifically, (\ref{pflow}) and (\ref{qflow}) define the active power and reactive power balance for each node, respectively. Voltage relation is enforced in (\ref{vrelations}), which is an exact formulation. The second-order cone relaxation is implemented in (\ref{socp}), which is relaxed from its equality form, defining apparent power flow for each cable. Constraints (\ref{vlimits}) and (\ref{ilimits}) define bus voltage limits and cable loading limits. Although there exist other convex relaxation models of the power flow relations such as semidefinite programming formulation \cite{Lavaei2012}, the second-order cone relaxation generally has a high accuracy and is exact under mild conditions for radial networks \cite{Farivar2013}, thus is widely used.

\begin{subequations}
\begin{IEEEeqnarray}{rl} \label{obj}
&\underset{(p_j^{ev},p_j^{hp},p_j^{pv},q_j^{ev},q_j^{pv},\forall j \in \mathcal{H})}{\text{minimize}}\Big\{c^{loss}\sum_{\forall (i, j) \in \mathcal{E}} r_{i j} \ell_{i j} \nonumber\\&   +  \sum_{\forall j \in \mathcal{H}}(c_j^{ev} p_j^{ev\downarrow} + c_j^{pv} p_j^{pv\downarrow} + c_j^{hp\uparrow} p_j^{hp\uparrow}+ c_j^{hp\downarrow} p_j^{hp\downarrow} )\Big\}
\end{IEEEeqnarray}
subject to
\begin{equation} \label{pbalance}
   p_j = p_j^{pv}-p_j^{load}-p_j^{ev}-p_j^{hp}  , \forall j \in \mathcal{H}
\end{equation}
\begin{equation} \label{qbalance}
   q_j = q_j^{pv}+q_j^{ev}-q_j^{hp}-q_j^{load} , \forall j\in \mathcal{H}
\end{equation}
\begin{equation} \label{pv1}
    (p_j^{pv})^2 +(q_j^{pv})^2 \leq (S_j^{pv})^2, \forall j\in \mathcal{H}
\end{equation}
\begin{equation} \label{pv2}
    - p_j^{pv} \tan \Phi^{pv,max} \leq q_j^{pv} \leq p_j^{pv} \tan \Phi^{pv,max}, \forall j\in \mathcal{H}
\end{equation}
\begin{equation} \label{ev1}
    (p_j^{ev})^2 +(q_j^{ev})^2 \leq (S_j^{ev})^2, \forall j\in \mathcal{H}
\end{equation}
\begin{equation} \label{ev2}
    - p_j^{ev} \tan \Phi^{ev,max} \leq q_j^{ev} \leq p_j^{ev} \tan \Phi^{ev,max}, \forall j\in \mathcal{H}
\end{equation}
\begin{equation} \label{hp}
    q_j^{hp} = p_j^{hp} \tan \Phi^{hp}, \forall j\in \mathcal{H}
\end{equation}
\begin{equation} \label{evcurtail}
    0 \leq p_j^{ev} \leq p_j^{ev,max}, p_j^{ev\downarrow} = p_j^{ev,max}-p_j^{ev}, \forall j\in \mathcal{H}
\end{equation}
\begin{equation}\label{pvcurtail}
    0 \leq p_j^{pv} \leq p_j^{pv,fore}, p_j^{pv\downarrow} = p_j^{pv,fore}-p_j^{pv}, \forall j\in \mathcal{H}
\end{equation}
\begin{equation}\label{hpcurtail1}
    0 \leq p_j^{hp} \leq p_j^{hp,max},  \forall j\in \mathcal{H}
\end{equation}
\begin{equation}\label{hpcurtail2}
    p_j^{hp\downarrow} \geq \max(0, p_j^{hp,set}-p_j^{hp}), \forall j\in \mathcal{H}
\end{equation}
\begin{equation}\label{hpcurtail3}
    p_j^{hp\uparrow} \geq \max(0, p_j^{hp}-p_j^{hp,set}), \forall j \in \mathcal{H}
\end{equation}
\begin{equation}\label{pflow}
p_{j}=\sum_{k: j \rightarrow k} P_{j k}-\sum_{i: i \rightarrow j}\left(P_{i j}-r_{i j} \ell_{i j}\right), \forall j \in \mathcal{N}
\end{equation}
\begin{equation}\label{qflow}
q_{j}=\sum_{k: j \rightarrow k} Q_{j k}-\sum_{i: i \rightarrow j}\left(Q_{i j}-x_{i j} \ell_{i j}\right), \forall j\in \mathcal{N}
\end{equation}
\begin{IEEEeqnarray}{rl}\label{vrelations}
v_{j}=v_{i}-2\left(r_{i j} P_{i j}+x_{i j} Q_{i j}\right)+&\left(r_{i j}^{2}+x_{i j}^{2}\right) \ell_{i j},\\
&\forall(i, j) \in \mathcal{E} \nonumber
\end{IEEEeqnarray}
\begin{equation}\label{socp}
P_{i j}^{2}+Q_{i j}^{2} \leq \ell_{i j} v_i , \forall(i, j) \in \mathcal{E}
\end{equation}
\begin{equation}\label{vlimits}
   v_j^{min} \leq v_j \leq v_j^{max},  \forall j\in \mathcal{N}
\end{equation}
\begin{equation}\label{ilimits}
   0 \leq \ell_{i,j} \leq \ell_{i,j}^{max},
\forall(i, j) \in \mathcal{E}
\end{equation}
\end{subequations}

In summary, this model can restore grid limits by making full use of DER flexibility. It is noteworthy that this model requires only a one-step-ahead forecast, without the need for a multi-period setting and a multi-period-ahead forecast. The time-coupling DERs are coped with by deviating as little as possible from their original schedules while prioritizing DERs with the most urgent demand using dynamic cost terms.
\subsection{Case Description}
This section presents a case study to evaluate the proposed procedure and DER coordination model for real-time grid operations. The simulation is performed on a Dutch LV network, shown in Fig. \ref{fig:epse}. A 100\% penetration scenario for EVs, PVs, and HPs is assumed for 67 connected residential households seeing the energy transition.
\begin{figure}[tbp]
    \centering
    \includegraphics[width=\linewidth]{figs/epse.png}
    \caption{Dutch LV network topology}
    \label{fig:epse}
\end{figure}

\begin{figure}[tbp]
    \centering
    \includegraphics[width=\linewidth]{figs/dataShow.pdf}
    \caption{Base load, PV, EV, HP profiles for a household for the selected day, starting from 8:00 with a 15-minute resolution}
    \label{fig:datashow}
\end{figure}

At the first step, a power flow program is implemented for a year with a 15-minute resolution. From which, a winter day with the most significant under-voltage and overloading issues is selected to evaluate the proposed real-time operation procedure and coordination model. The base load profiles are based on anonymous smart meter data from a Dutch DSO, while the PV generation profile is calculated from an annual solar irradiance dataset from the Royal Netherlands Meteorological Institute assuming a 4.5kW PV installation for each household. The HP consumption profiles are based on anonymous smart meter gas usage data from a Dutch DSO, assuming 70\% of gas is used for domestic heating \cite{Ruhnau2019}. Literature \cite{Brunner2016, Sichilalu2017} can be referred to for HP and water tank models. For EV charging, the Monto Carlo simulation is performed by assuming normally distributed EV arrival time and a fixed electricity demand of 7.5kWh per day. A 230V, 16A charger with a 90\% efficiency is assumed for EV charging. For illustration purpose, base load, PV, EV, and HP profiles for a random household on the selected winter day is presented in Fig. \ref{fig:datashow}, where the coincidence of uncontrolled EV charging and HP consumption during nighttime is likely to result in under-voltage and overloading issues in the grid.

\subsection{Simulation Flowchart}
\begin{figure}[tbp]
    \centering
    \includegraphics[width=.75\linewidth]{figs/simFlow.pdf}
    \caption{Simulation flowchart}
    \label{fig:simFlow}
\end{figure}

\begin{figure}[tbp]
    \centering
    \includegraphics[width=\linewidth]{figs/CDF_voltage.pdf}
    \caption{Bus voltage CDF}
    \label{fig:voltageCDF}
\end{figure}

\begin{figure}[tbp]
    \centering
    \includegraphics[width=\linewidth]{figs/CDF_loading.pdf}
    \caption{Cable loading CDF}
    \label{fig:loadingCDF}
\end{figure}

This section introduces the setting for the simulation on the selected day. The programs are implemented on a laptop with an Intel Core i5-9300H CPU running at 2.4GHz with 8GB RAM. The power flow program is implemented with Pandapower \cite{pandapower.2018}, while the optimization model is solved with Gurobi \cite{gurobi} with Python interface. The power flow program takes on average 1 second to solve, while the optimization model takes on average 4 seconds. It is noteworthy that a much higher resolution simulation could be performed under such short execution time, a 15-minute resolution simulation is still used considering the granularity of the available data.

The simulation process is described in Fig. \ref{fig:simFlow}, which starts with collecting various datasets as specified. At each simulation time step, a power flow program is implemented, followed by the DER coordination model in case any grid limit is violated. The optimized schedule from the coordination model is implemented, following which the EV charging and HP consumption profiles are updated for subsequent time slots to catch up on the deviation from the original schedules. For EVs, the charging process is extended in case the charging is curtailed. For HPs, the consumption is adjusted within their full capacity ranges to cancel out the deviation with the least possible delay. Finally, the EV state of charge (SoC) and water tank temperature are updated, which are used to calculate their corresponding cost terms in the DER coordination model and to evaluate user comfort.

\begin{table}[tbp]
    \centering
        \caption{Cost terms for the DER coordination model}
    \begin{tabular}{|c|c|c|} \hline
  \textbf{Cost term} & \textbf{Value} & \textbf{Note}\\ \hline
       $c^{loss}$  & 32\texteuro/MWh & Average day-ahead price \\ \hline
       $c^{pv}$  & 200\texteuro/MWh & Behind-the-meter\\\hline
       $c^{ev}$  & $\frac{c_0(1-SoC)}{(T_{max}-t+1)}$\texteuro/MWh & $c_0=1440, T_{max}=96$\\\hline
       $c^{hp\uparrow}$  & $\max(10,150\Delta T)$\texteuro/MWh & $\Delta T$ normalized \\\hline
       $c^{hp\downarrow}$  & $\max(10,-150\Delta T)$\texteuro/MWh & $\Delta T$ normalized \\\hline
    \end{tabular}
    \label{tab:cost}
\end{table}

For EVs, the cost terms are assumed to be functions of SoC and the remaining charging time. A unified departure time is assumed for the EVs, at the end of the simulation period. User satisfaction is determined by whether the EVs are fully charged before departure. For HPs, water tanks are assumed to supply hot water and space heating. A water tank of 1000L is assumed for households with no more than 5kW heating demand, while an additional 200L is assumed for each kW above 5kW \cite{Williams2012}. The water tanks are allowed for ${5}^{\circ} \mathrm{C}$ deviation up and down from the reference temperature, while not affecting user comfort \cite{Brunner2016, Sichilalu2017}. The cost terms are assumed as linear functions of water tank temperature deviation with a lower bound. User satisfaction is evaluated with respect to the water tank temperature limits. The used cost terms for this case study are listed in Table \ref{tab:cost}. It is noteworthy that in real-world applications, these cost terms can be derived through extensive simulation to guarantee the model's performance, or directly from end-users from a market perspective. With these model settings, the simulation results are presented in the following sections.

\subsection{Network Issues}

\begin{figure}[tbp]
    \centering
    \includegraphics[width=\linewidth]{figs/comparison.pdf}
    \caption{Voltage magnitude, EV/HP power and EV reactive power generation for a time slot}
    \label{fig:compare}
\end{figure}

As described above, the major issues envisioned in this case study are under-voltage and overloading, which have been confirmed in Fig. \ref{fig:voltageCDF} and Fig. \ref{fig:loadingCDF} respectively. The reported grid operation states are both from the power flow programs, instead of the optimization model. Using the DER coordination model, both issues have been completely resolved, where the bus voltage stays higher than 0.9 pu and the cable loading percentage is no more than 70\%. Fig. \ref{fig:compare} presents the voltage magnitude, EV and HP power demand, and EV reactive power generation for a time slot before and after implementation of the DER coordination model, which suggests that the restoration of voltage magnitude is attributed to curtailed EV and HP power demand, and EV reactive power compensation.

\subsection{User Satisfaction}

\begin{figure}[tbp]
    \centering
    \includegraphics[width=\linewidth]{figs/evsatisfaction.pdf}
    \caption{EV charging profiles for various households}
    \label{fig:evsatisfy}
\end{figure}

\begin{figure}[tbp]
    \centering
    \includegraphics[width=\linewidth]{figs/hpsatisfaction.pdf}
    \caption{Water tank temperature deviation for various households}
    \label{fig:hpsatisfy}
\end{figure}

In this study, two indicators are used to evaluate user satisfaction: whether EVs have been fully charged, and whether the temperature deviation of water tanks has exceeded the upper and lower limits. Fig. \ref{fig:evsatisfy} shows 3 different charging profiles, where the first EV's charging process is not interrupted, while the rest two have been curtailed during the charging process. However, all EVs have been fully charged at the end of the simulation. Fig. \ref{fig:hpsatisfy} presents the water tank temperature changes for various households. The temperature has been quickly restored to the normal state after a HP curtailment occurs. It is shown that all water tanks have not violated the temperature limits (${5}^{\circ} \mathrm{C}$ up and down). Combining these results, it is demonstrated that user satisfaction has been kept using this DER coordination model.

\subsection{SOCP Relaxation Accuracy}
\begin{figure}[tbp]
    \centering
    \includegraphics[width=\linewidth]{figs/socp.pdf}
    \caption{Voltage magnitude and cable loading from both SOCP relaxation model and power flow program}
    \label{fig:SOCP}
\end{figure}

This section deals with the accuracy of the adopted SOCP relaxation of the power flow relations. Fig. \ref{fig:SOCP} shows the voltage magnitude and cable loading for one time slot from the DER coordination model and the power flow program. It is seen that the voltage has been precisely captured by the SOCP model, while the cable loading has been over-estimated, due to the relaxation of the apparent power flow relation. However, an important finding is that the SOCP relaxation model ensures the feasibility of the optimized active and reactive power setpoints, which implies that when these control actions from the DER coordination model are implemented in the grid, it would not result in loading or voltage issues that are not captured in the SOCP model.

\section{Introduction}
\input{1.intro}
\section{DER Coordination model}\label{sec:model}
\input{2.mathmodel}
\section{Case study}\label{sec:case}
\input{3.case}
\section{Conclusions}\label{sec:conclusion}
\input{4.conclusions}


\section{INTRODUCTION}
In recent years, much research has been focused on autonomous driving technology. LiDAR is one of the most important sensors for perception tasks such as drivable region segmentation, object detection and vehicle tracking. Different from images captured by cameras, point cloud generated by LiDARs supplies 3D spatial information of the objects in the form of (X, Y, Z) coordinates and intensity. This alleviates the barrier of distance estimation and makes 3D object detection or tracking much more accurate. However, the price of high-resolution LiDARs is much higher than their low-resolution counterparts. The specifications of the most popular Velodyne 64-line LiDAR HDL-64E and 16-line LiDAR VLP-16 are compared in Table~\ref{tab:lidar_cmp}. As we can see, the cost of a low-resolution LiDAR is only about 1/18 of the high-resolution ones. Therefore, it is more economical to consider low-resolution LiDARs in order to build low-cost autonomous driving systems. However, it is a major challenge to perform object detection from the point cloud produced by a low-resolution LiDAR since it is too sparse to even show the shapes of objects. As illustrated in Fig.~\ref{fig:dmap_comp}, we can barely find objects from the depth map captured from a 16-line LiDAR, while in the 64-line LiDAR objects are more visible.

\begin{table}[htbp]
    \centering
    \caption{A comparison of VLP-16 and HDL-64E LiDARs}
    \label{tab:lidar_cmp}
    \begin{tabular}{|l|l|l|}
    \hline
    \textbf{LiDAR type} & \textbf{VLP-16} & \textbf{HDL-64E} \\
    \hline
    \hline

    Channel number & 16 & 64 \\
    \hline

    Res.(vertical / horizontal) & 2\textdegree / 0.2\textdegree & 0.4\textdegree / 0.1728\textdegree\\
    \hline

    Power (Watts) & 8 & 60 \\
    \hline
    Price (USD) & \$4'000 & \$75'000 \\
    \hline
    \end{tabular}
\end{table}

\begin{figure*}[htbp]
    \captionsetup[subfigure]{aboveskip=-0.5pt,belowskip=-1pt}
    \centering
    \begin{subfigure}{0.32\textwidth}
        \includegraphics[width=0.98\textwidth]{fig/000943_img.jpg}
        \caption{image scenario 1}
    \end{subfigure}%
    \begin{subfigure}{0.32\textwidth}
        \includegraphics[width=0.98\textwidth]{fig/000948_img.jpg}
        \caption{image scenario 2}
    \end{subfigure}%
    \begin{subfigure}{0.32\textwidth}
        \includegraphics[width=0.98\textwidth]{fig/000965_img.jpg}
        \caption{image scenario 3}
    \end{subfigure}%
    \newline
    \begin{subfigure}{0.32\textwidth}
        \includegraphics[width=0.98\textwidth]{fig/000943_64_1.jpg}
        \caption{dense depth map scenario 1}
    \end{subfigure}%
    \begin{subfigure}{0.32\textwidth}
        \includegraphics[width=0.98\textwidth]{fig/000948_64_1.jpg}
        \caption{dense depth map scenario 2}
    \end{subfigure}%
    \begin{subfigure}{0.32\textwidth}
        \includegraphics[width=0.98\textwidth]{fig/000965_64_1.jpg}
        \caption{dense depth map scenario 3}
    \end{subfigure}%
    \newline
    \begin{subfigure}{0.32\textwidth}
        \includegraphics[width=0.98\linewidth]{fig/000943_16_1.jpg}
        \caption{sparse depth map scenario 1}
    \end{subfigure}%
    \begin{subfigure}{0.32\textwidth}
        \includegraphics[width=0.98\linewidth]{fig/000948_16_1.jpg}
        \caption{sparse depth map scenario 2}
    \end{subfigure}%
    \begin{subfigure}{0.32\textwidth}
        \includegraphics[width=0.98\linewidth]{fig/000965_16_1.jpg}
        \caption{sparse depth map scenario 3}
    \end{subfigure}%
    \caption{Comparison of depth map from 16-line LiDAR (bottom) and 64-line LiDAR (middle) to their RGB image (top), on which  red boxes represent the short-range vehicles, orange boxes show the medium range vehicles and the far vehicles are marked by blue boxes. Green boxes illustrate the occluded vehicles.}
    \label{fig:dmap_comp}
\end{figure*}

\section{RELATED WORKS}

\subsection{Low-Resolution LiDAR and Depth Completion}
Some research works focused on segmentation using low-resolution LiDARs. \cite{gigli2020road} introduced the local normal vector for the LiDAR’s spherical coordinates as an input channel. Based on the existing LoDNN architectures \cite{calt2017fast}, its road segmentation performance using low-resolution LiDAR was close to that from high-resolution LiDAR within a reasonable degradation. A supervised domain adaptation was utilized by \cite{elhadidy2020improved} to predict the low-resolution point cloud into high-resolution point cloud in spherical coordinate and further evaluated the results in 3D semantic segmentation task. Low-resolution LiDARs had been also employed for object tracking tasks. In \cite{pino2017low}, a LiDAR-based system was proposed for estimation of actual positions and velocities of the detected vehicles.
Some other works utilized depth completion for 2D object detection, such as \cite{seikavandi2020deep} and \cite{farahnakian2020fusing}. In \cite{seikavandi2020deep}, a weighted depth filling algorithm was proposed to make the high-resolution (HDL-64E) LiDAR depth map even denser. Subsequently, this dense depth map was concatenated with the corresponding RGB image as the input of YOLOv3 \cite{redmon2018yolov3} network for 2D object detection. Similarly, the authors of \cite{farahnakian2020fusing} introduced a self-supervised depth completion network to fill the high-resolution depth map before detection 2D objects.

\subsection{High Resolution LiDAR for BEV Object Detection} Nearly all state-of-the-art object detectors utilize high-resolution LiDAR. In \cite{beltran2018birdnet}, it first transformed the point cloud into \textcolor{black}{Bird-Eye View (BEV)} map, and then extracted the ground and proposed the objects in two branches separately. Finally, the objects were predicted by a post-processing block. \cite{barrera2020birdnet+} further refined the previous version into an end-to-end model and achieved better performance. Single-stage detector, PIXOR, was proposed in \cite{yang2018pixor} by using 2D convolution on the voxelized BEV map. Without any anchor, it achieved real-time processing speed.

As mentioned earlier, due to the extreme sparsity, low-resolution LiDAR depth map does not supply enough shape information of the objects, but some sub-samples of the precise depth information. Meanwhile, the RGB image supplies rich context information. Thus, we argue that when fusing sparse depth map and RGB image together, 3D object detection becomes possible.

\section{PROPOSED 3D OBJECT DETECTION FRAMEWORK}

In this paper, we investigate the possibility of low-resolution LiDAR usage in BEV object detection task. In Fig.~\ref{fig:dmap_comp}, red box, orange box and blue box represent the vehicle in short range, medium range and long range respectively. For short range vehicles, their shapes are clearly visible from dense depth maps. In sparse depth maps, the shapes are very blurry but still recognizable since the number of points hitting on the vehicles is still large enough. Concerning the medium and long range vehicles (in orange and blue boxes), we can only get a small number of points even using 64-line LiDAR. While in the sparse depth map from 16-line LiDAR, the number of hit points is few to none. Taking the medium range vehicles in orange boxes in Fig.~\ref{fig:dmap_comp} (h) for example, it is easy to recognize them as obstacles due to sharp distance distinction but difficult to recognize them as vehicles. This also applies to vehicles with occlusion (green boxes in Fig.~\ref{fig:dmap_comp} (c), (f) and (i)). The long range vehicles in blue box (in both Fig.~\ref{fig:dmap_comp} (e) and (h)) get too few points to be correctly localized and classified. According to the analysis above, we found that unlike the depth map from 64-line LiDAR, 16-line LiDAR depth map does not show reliable context information but accurate distance information. This implies that 16-line LiDAR depth map is more useful for depth estimation rather than context information extraction. Therefore, to better use the information from 16-line depth map, we put a depth completion network prior to the object detector to generate a dense depth map with context information. After the dense depth map is generated, it is sent to 3D object detector, as demonstrated in Fig.~\ref{fig:framework}.

\begin{figure}[htbp]
    \centering
    \includegraphics[width=\linewidth]{fig/framework.jpg}
    \caption{The proposed framework for 3D object detection using low-resolution point cloud and RGB image}
    \label{fig:framework}
\end{figure}

\subsection{Depth Completion Network}\label{sec:depth_com}
The depth completion network aims to fill the sparse depth map from 16-line LiDAR point cloud with the help of RGB image. The state-of-the-art depth completion network \cite{ma2019self} is adopted here with some modifications. It requires two inputs, RGB image and low-resolution sparse depth map. The RGB image supplies the context information in detail, while the sparse depth map supplies the precise depth information for some pixels on the image. The sensor fusion strategy adopted here is also referred as early fusion. To make the network more compact, we first replace the ResNet-34 backbone with ResNet-18. For performance improvement, global attention modules and an Atrous Spatial Pyramid Pooling (ASPP) \cite{chen2018deeplabv3+} module are placed to bridge the encoder and decoder.

As shown in Fig.~\ref{fig:gam}, the global attention module is used to extract global context information of the feature map by global pooling layer, and then fuse the global information back to guide the feature learning. Through adding this module, the global information is merged into features without up-sampling layer. This helps the decoder part to achieve better performance. Besides, an ASPP module is placed between encoder and decoder, with each convolution dilated rate 2, 4, 8 and 16. The ASPP module concatenates feature maps with different fields of perception, so that decoder has a better understanding of the context information.

\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.6\linewidth]{fig/gam.jpg}
    \caption{Structure of global attention module}
    \label{fig:gam}
\end{figure}

The loss function of depth completion network is the Mean Square Error (MSE) between the predicted depth map and the ground truth.

\subsection{Object Detection Network}
The object detection network used in this framework is PIXOR \cite{yang2018pixor}. Its main idea is to take advantage of 2D convolution and anchor-free network to realize super-fast point cloud object detection in BEV. PIXOR consists of two steps. The first step is to reform the representation of input point cloud. It reduces 3 degrees of freedom to 2 in BEV, and extracts the 3rd freedom (z or height) as another input feature map channel. So that 2D convolution instead of 3D convolution is necessary to greatly decrease the computation complexity. The second step is to feed the reformed input feature map into an anchor-free one-stage object detector network. For the highly efficient computation on dense predictions, a fully convolutional architecture is utilized to build the backbone and header of PIXOR. Without any pre-defined anchors and proposals, PIXOR outputs the predicted class and orientation from header in a single network.

Concerning the loss function, the total loss of object detection consists of the classification loss and the regression loss (Eq.~\ref{eq:detect}), where $\lambda_{cls}$ and $\lambda_{reg}$ are the corresponding coefficients. The classification loss $\mathcal{L}_{cls}$ targets to correctly predict the object (cars in our case) and the regression loss $\mathcal{L}_{reg}$ aims to refine the size, center and the orientation of the predicted bounding boxes.

\begin{equation}
    \mathcal{L}_{detect} = \lambda_{cls}\mathcal{L}_{cls} + \lambda_{reg}\mathcal{L}_{reg}
    \label{eq:detect}
\end{equation}

\subsection{Implementation Details}
The depth completion network is first trained on KITTI Depth Completion dataset. The depth completion network is trained with batch size of 4, and learning rate starts at 1e-4 which decreases every 5 epochs. The total number of training epoch is 10. After training the depth completion network and keeping as it is, we move on to train the object detector from scratch. The KITTI Object Detection dataset has been split into training and validation parts according to \cite{qian2020end}.  The optimizer is Adam, with batch size 8. The learning rate starts at 1e-3 and reduces by a factor of 2 when the validation loss does not decrease. Finally, we fine-tune the entire framework with both depth completion network and object detection network together, with 16-line point cloud and images as input and vehicles in BEV as output.

\section{EXPERIMENT RESULTS}

\subsection{Evaluation Dataset}
Training and evaluation of the whole framework both employ KITTI dataset (both Depth Completion and Object Detection). Before feeding into the framework mentioned above, the point clouds are down-sampled to emulate the VLP-16 low-resolution LiDAR. KITTI depth completion dataset contains 85,898 training data and 1,000 selected validation data. Its ground truth is produced by aggregating consecutive LiDAR scan frames into a semi-dense depth map, about 30\% annotated pixels. KITTI object detection dataset has 7,481 training data and 7,518 testing data. Evaluation is categorized into three regimes: easy, moderate and hard, representing objects at different occlusion and truncation levels.

\subsection{Depth Completion Performance Evaluation}
As described in Sec.~\ref{sec:depth_com}, in order to enhance the depth completion performance, multiple GAM modules have been added to bridge the encoder and the decoder of depth completion network. The performance comparison on validation dataset is illustrated in Tab.~\ref{tab:dc_perform}. Adding GAM modules results in the performance improvement of about 3.6\% and 7.0\% measured by Root Mean Square Root (RMSE) and Mean Average Error (MAE) respectively.

\begin{table}[htbp]
    \centering
    \caption{Depth completion performance comparison with and without GAM modules}
    \label{tab:dc_perform}
    \begin{tabular}{|c|c|c|}
    \hline
    with GAMs & RMSE (1/mm) & MAE (1/mm) \\
    \hline
    Yes & 1592.74 & 537.81\\
    \hline
    No & 1651.68 & 578.14\\
    \hline
    \end{tabular}
\end{table}

Fig.~\ref{fig:dc_cmp} (b) and (c) demonstrate the predicted depth maps of depth completion networks with and without GAM modules respectively. And the bottom figure shows the ground truth. In this example, the depth map from depth completion network with GAM module gives objects slightly better shape representation.

\begin{figure}[htbp]
    \centering
    \begin{subfigure}{0.5\linewidth}
        \includegraphics[width=0.98\linewidth]{fig/img.jpg}
        \vspace{-0.1cm}
        \caption{Image}
        \includegraphics[width=0.98\linewidth]{fig/gt.jpg}
        \vspace{-0.1cm}
        \caption{Depth ground truth}
        \captionsetup{labelformat=empty}
    \end{subfigure}%
    \begin{subfigure}{0.5\linewidth}
        \includegraphics[width=0.98\linewidth]{fig/noGam.jpg}
        \vspace{-0.1cm}
        \caption{Depth map without GAM}
        \includegraphics[width=0.98\linewidth]{fig/withGam.jpg}
        \vspace{-0.1cm}
        \caption{Depth map with GAM}
    \end{subfigure}%
    \caption{Comparison of depth map from 16-line LiDAR with and without GAM module to their RGB image and ground truth}
    \label{fig:dc_cmp}
    \vspace{-0.6cm}
\end{figure}

\begin{table*}[htbp]
    \centering
    \caption{\textbf{BEV performance comparison on KITTI object detection validation dataset.} This table shows $AP_{BEV}$ (in \%) of the car category, corresponding to average precision of the bird's-eye view.}

    \begin{tabular}{|c|c|c|c|c|c|c|c|}
    \hline
    \multirow{2}*{Detection Networks} & \multirow{2}*{Input} & \multicolumn{3}{|c|}{IoU=0.5} & \multicolumn{3}{|c|}{IoU=0.7} \\
    \cline{3-8}
     & & Easy & Moderate & Hard & Easy & Moderate & Hard \\
    \hline
    PIXOR\cite{yang2018pixor} & LiDAR only (64-line) & 94.2 & 86.7 & 86.1 & 85.2 & 81.2 & 76.1 \\
    \hline
    PIXOR & LiDAR only (16-line) & 60.7 & 51.2 & 46.8 & 53.8 & 47.1 & 39.1 \\
    \hline
    \textbf{Ours} & \textbf{LiDAR (16-line) + Camera} & \textbf{89.0} & \textbf{75.8} & \textbf{68.1} & \textbf{75.4} & \textbf{61.2} & \textbf{55.2} \\
    \hline
    \end{tabular}
    \label{tab:results}
\end{table*}

\begin{figure}[htbp]
    \centering
    \begin{subfigure}{0.5\linewidth}
        \includegraphics[width=0.98\linewidth]{fig/img_25.jpg}
        \includegraphics[width=0.98\linewidth]{fig/fov_25.jpg}
    \end{subfigure}%
    \begin{subfigure}{0.5\linewidth}
        \includegraphics[width=0.98\linewidth]{fig/img_117.jpg}
        \includegraphics[width=0.98\linewidth]{fig/fov_117.jpg}
    \end{subfigure}%
    \newline
    \begin{subfigure}{0.5\linewidth}
        \includegraphics[width=0.98\linewidth]{fig/img_52.jpg}
        \includegraphics[width=0.98\linewidth]{fig/fov_52.jpg}
    \end{subfigure}%
    \begin{subfigure}{0.5\linewidth}
        \includegraphics[width=0.98\linewidth]{fig/img_33.jpg}
        \includegraphics[width=0.98\linewidth]{fig/fov_33.jpg}
    \end{subfigure}%
    \caption{Visualization of object detection from the proposed framework, where the green boxes are ground truth and the blue boxes represent the predicted results}
    \label{fig:result}
\end{figure}

\subsection{Object Detection Performance Evaluation}
The performance of our framework on KITTI object detection validation dataset is illustrated in Tab.~\ref{tab:results} and Fig.~\ref{fig:result}. The results are shown in two circumstances IoU=0.5 and IoU=0.7 respectively. When IoU=0.5, our framework achieves 89.0\%, 75.8\% and 68.1\% detection accuracy for easy, moderate and hard cases respectively. While in case of IoU=0.7, the prediction accuracy is decreased to 75.4\%, 61.2\% and 55.2\% respectively. Compared to feeding 16-line point cloud directly into PIXOR, our framework pulls up the detection accuracy significantly in all cases. If compared to PIXOR with 64-line point cloud as input, the performance of our framework is relatively comparable in easy and moderate cases. But in hard case, the prediction accuracy drops around 20\% in both IoU criteria. The precision-recall curve is demonstrated in Fig.~\ref{fig:curve}.

\textcolor{black}{In regard to the computations, when running on RTX 2080Ti, the inference time of proposed network is 25.4ms or 39.8 frames per second (fps). Besides, as the network is aiming for embedded systems, we also tested it on DRIVE PX2 that contains two discrete Pascal GPUs. The inference latency for each point cloud frame is 581.7ms. If two GPUs run as two threads, the throughput increases to 3.4fps.}

\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.9\linewidth]{fig/curve.jpg}
    \caption{Precision-Recall curve of the proposed framework on KITTI val dataset.}
    \label{fig:curve}
\end{figure}

\begin{comment}
\subsection{Analysis and Future Work}
\subsubsection{Location prediction}
The depth completion task in KITTI adopts RMSE and MAE as benchmark. However, if training with RMSE or MAE as loss function, the network cannot distinguish object boarder pixels as foreground pixels or background pixels. Instead, to reach higher RMSE and MAE, depth completion network attempts to predict the depth of the car edge as a value between the distance of the vehicle and the distance of the ground in front of the car (Fig.~\ref{fig:dc_out}). This "long-tail" boundary effect might make the object localization more difficult. In our future work, a new benchmark will be proposed for LiDAR depth completion in order to remove the boundary effects mentioned above.

\begin{figure}[htbp]
    \centering
    \begin{subfigure}[b]{\linewidth}
        \centering
        \includegraphics[width=0.9\linewidth]{fig/0000000045.jpg}
        \caption{Image}
    \end{subfigure}
    \begin{subfigure}[b]{0.45\linewidth}
        \centering
        \includegraphics[width=0.9\linewidth]{fig/dc_3d_local.jpg}
        \caption{Predicted depth map}
    \end{subfigure}
    \begin{subfigure}[b]{0.45\linewidth}
        \centering
        \includegraphics[width=0.9\linewidth]{fig/gt_3d_local.jpg}
        \caption{Ground truth}
    \end{subfigure}
    \caption{3D visualization of long-tail effect in both the predicted depth map and ground truth of depth completion network}
    \label{fig:dc_out}
\end{figure}

\subsubsection{Vehicles in long range}
Due to the extreme sparsity of low-resolution LiDAR depth map, vehicles in long distances from the LiDAR are only visible in images. Thus, for these vehicles, the depth map assisted depth prediction is practically downgraded as depth estimation from image only. This explains a sharp performance drop on hard cases in KITTI dataset.
\end{comment}

\section{Conclusion}
This paper presents a framework that enables 3D object detection using a low-resolution LiDAR. By cascading a depth completion network with an object detector, it first converts the sparse point cloud into a denser depth map that is subsequently processed for 3D object detection. It demonstrates 3D object detection with only a 16-line LiDAR and a camera. When evaluated on KITTI dataset, the proposed solution achieves high accuracy in object detection for both easy and moderate cases, comparable to the benchmarks using 64-line LiDARs.