# Searching for Lepton Flavour (Universality) Violation and Collider Signals
from a Singly-Charged Scalar Singlet

Andreas Crivellin<EMAIL_ADDRESS>CERN Theory Division, CH–1211
Geneva 23, Switzerland  Physik-Institut, Universität Zürich,
Winterthurerstrasse 190, CH–8057 Zürich, Switzerland  Paul Scherrer Institut,
CH–5232 Villigen PSI, Switzerland  Fiona Kirk<EMAIL_ADDRESS>Physik-
Institut, Universität Zürich, Winterthurerstrasse 190, CH–8057 Zürich,
Switzerland  Paul Scherrer Institut, CH–5232 Villigen PSI, Switzerland
Claudio Andrea Manzari<EMAIL_ADDRESS>Physik-Institut,
Universität Zürich, Winterthurerstrasse 190, CH–8057 Zürich, Switzerland  Paul
Scherrer Institut, CH–5232 Villigen PSI, Switzerland  Luca Panizzi
<EMAIL_ADDRESS>Department of Physics and Astronomy, Uppsala
University, Box 516, SE-751 20 Uppsala, Sweden

###### Abstract

A singly charged $SU(2)_{L}$ singlet scalar ($\phi^{\pm}$) can only have
flavour off-diagonal couplings to neutrinos and charged leptons, therefore
necessarily violating lepton flavour (universality). In fact, it gives a
(necessarily constructive) tree-level effect in $\ell\to\ell^{\prime}\nu\nu$
processes, while contributing to charged lepton flavour violating only at the
loop-level. Therefore, it can provide a common explanation of lepton flavour
universality violation in $\tau\to\mu\nu\nu/\tau(\mu)\to e\nu\nu$ and of the
Cabibbo Angle Anomaly. For the best fit region we predict ${\rm Br}[\tau\to
e\gamma]$ to be of the order of a few times $10^{-11}$ while ${\rm Br}[\tau\to
e\mu\mu]$ can be of the order of $10^{-9}$ for order one couplings.
Furthermore, we derive a coupling-independent lower limit on the scalar mass
of $\approx 200\,$GeV by recasting LHC slepton searches. In the scenario
preferred by low energy precision data, the lower limit is even strengthened
to $\approx 300\,$GeV. Finally, we point out that our scenario can be tested
by reinterpreting DM mono-photon searches at future $e^{+}e^{-}$ colliders.

††preprint: CERN-TH-2020-213, PSI-PR-20-24, UZ-TH 55/20

##  I Introduction

While the LHC, in its quest for discovering beyond the Standard Model (SM)
physics, has not discovered any new particles directly Butler (2017); Masetti
(2018), intriguing indirect hints for the violation of lepton flavour
universality (LFU) were accumulated. In particular, global fits to $b\to
s\ell^{+}\ell^{-}$ Aaij  _et al._ (2014a, b, 2015a, 2016); Khachatryan  _et
al._ (2016); ATLAS-CONF-2017- (023); CMS-PAS-BPH-15- (008); Aaij  _et al._
(2017) and $b\to c\tau\nu$ Lees  _et al._ (2012, 2013); Aaij  _et al._ (2015b,
2018a, 2018b); Abdesselam  _et al._ (2019) data point convincingly towards new
physics (NP) with a significance of $>\\!5\,\sigma$ Capdevila  _et al._
(2018); Altmannshofer  _et al._ (2017); D’Amico  _et al._ (2017); Ciuchini
_et al._ (2017); Hiller and Nisandzic (2017); Geng  _et al._ (2017); Hurth
_et al._ (2017); Alok  _et al._ (2017); Algueró  _et al._ (2019); Aebischer
_et al._ (2020); Ciuchini  _et al._ (2019) and $>\\!3\,\sigma$ Amhis  _et al._
(2017); Murgui  _et al._ (2019); Shi  _et al._ (2019); Blanke  _et al._
(2019); Kumbhakar  _et al._ (2020), respectively. In addition, also the long
standing tension in the anomalous magnetic moment of the muon Bennett  _et
al._ (2006); Aoyama  _et al._ (2020) and the deficit in first row CKM
unitarity Belfatto  _et al._ (2020); Grossman  _et al._ (2020); Coutinho  _et
al._ (2020); Crivellin and Hoferichter (2020); Crivellin  _et al._ (2020a);
Kirk (2020); Alok  _et al._ (2020); Crivellin  _et al._ (2020b); Shiells  _et
al._ (2020); Seng  _et al._ (2020), known as the Cabibbo Angle Anomaly (CAA),
can be interpreted as signs of LFU violation.

Interestingly, not only the CAA can be explained by a constructive NP
contribution to the SM $\mu\to e\nu_{\mu}\bar{\nu}_{e}$ amplitude, but also
the analogous tau decays $\tau\to\mu\nu_{\tau}\bar{\nu}_{\mu}$ prefer a
constructive NP effect at the $2\,\sigma$ level Amhis  _et al._ (2019). Such
an effect can be most naturally generated at tree-level, as loop effects are
strongly constrained by LEP and LHC data. Furthermore, as data require NP to
interfere constructively with the SM, there are only four possible NP
candidates111Also a $SU(2)_{L}$ triplet scalar gives rise to a SM-like
amplitude but interferes destructively.: vector-like leptons Crivellin  _et
al._ (2020a), a left-handed vector $SU(2)_{L}$ triplet Capdevila  _et al._
(2020), a left-handed $Z^{\prime}$ with flavour violating couplings Buras  _et
al._ and a singly charged $SU(2)_{L}$ singlet scalar. Interestingly, the last
option even gives a necessarily constructive effect and due to hermiticity of
the Lagrangian automatically violates lepton flavour (universality).
Furthermore, as a singly charged scalar cannot couple to quarks and only
generates charged lepton flavour violation at the loop level, it is weakly
constrained experimentally by other processes and can therefore potentially
explain the CAA and the hints for LFU violation in $\tau$ decays. This letter
is thus dedicated to the study of the phenomenology of the singly charged
$SU(2)_{L}$ singlet scalar in the light of the hints for LFU violation.

Singly charged scalars have been proposed within the Babu-Zee model Zee
(1986); Babu (1988) and studied in Refs. Krauss  _et al._ (2003); Nebot  _et
al._ (2008); Cai  _et al._ (2015); Cheung and Seto (2004); Ahriche  _et al._
(2014); Chen  _et al._ (2014); Ahriche  _et al._ (2016); Herrero-Garcia  _et
al._ (2014); Centelles Chuliá  _et al._ (2018); Babu  _et al._ (2020) as part
of a larger NP spectrum, mostly with the aim of generating neutrino masses at
loop-level. Here, we focus on the SM supplemented only by the singly charged
scalar (which constitutes a UV complete model) and perform a comprehensive
analysis of flavour and collider constraints in the context of the existing
hints for LFU violation.

##  II Model and Observables

As motivated in the introduction, we supplement the SM by a $SU(2)_{L}\times
SU(3)_{C}$ singlet $\phi^{+}$ with hypercharge +1. Interestingly, this allows
only for Yukawa-type interaction with leptons

$\displaystyle\mathcal{L}=\mathcal{L}_{\rm
SM}-\left({\xi_{ij}}/{2}\,\bar{L}^{c}_{a,i}\,\varepsilon_{ab}\,L_{b,j}\,\Phi^{+}+{\rm
h.c.}\right)\,,$ (1)

but not with quarks. Here $L$ is the left-handed $SU(2)_{L}$ lepton doublet,
$c$ stands for charge conjugation, $a$ and $b$ are $SU(2)_{L}$ indices, $i$
and $j$ are flavour indices and $\varepsilon_{ab}$ is the two-dimensional
anti-symmetric tensor. Note that without loss of generality, $\xi_{ij}$ can be
chosen to be anti-symmetric in flavour space, $\xi_{ji}=-\xi_{ij}$, such that
$\xi_{ii}=0$ and our free parameters are $\xi_{12}$, $\xi_{13}$ and
$\xi_{23}$. In addition, there can be a coupling to the SM Higgs doublet
$\lambda H^{\dagger}H\phi^{+}\phi^{-}$ which contributes to the mass
$m_{\phi}$ but otherwise only has a significant impact on $h\to\gamma\gamma$.

Figure 1: Feynman diagrams showing the contribution of $\phi^{\pm}$ to $\mu\to
e\nu_{\mu}\bar{\nu}_{e}$, $\mu\to\mathrm{e}\gamma$ and $\tau\to\mu ee$. The
corresponding diagrams for analogous processes with different flavours are not
depicted but can be deduced by straightforward substitutions.

###  II.1 $\ell\to\ell^{\prime}\nu\nu$

The SM decay of a charged lepton into a lighter one and a pair of neutrinos is
modified at tree-level in our model. Applying Fierz identities (see e.g. Ref.
Nieves and Pal (2004)) one can remove the charge conjugation and transform the
amplitude to the $V-A$ structure of the corresponding SM amplitude. Taking
only into account interfering effects with the SM we have

$\displaystyle\delta(\ell_{i}\to\ell_{j}\nu\nu)=\frac{\mathcal{A}_{NP}(\ell_{i}\to\ell_{j}\nu_{i}\bar{\nu}_{j})}{\mathcal{A}_{SM}(\ell_{i}\to\ell_{j}\nu_{i}\bar{\nu}_{j})}=\frac{\left|\xi_{ij}^{2}\right|}{g_{2}^{2}}\frac{m_{W}^{2}}{m_{\phi}^{2}}\,.$
(2)

This has to be compared to Amhis  _et al._ (2019)

$\displaystyle\begin{split}{\left.{\frac{{{\cal
A}(\tau\to\mu\nu\bar{\nu})}}{{{\cal A}(\mu\to
e\nu\bar{\nu})}}}\right|_{{\rm{EXP}}}}&=1.0029(14){\mkern 1.0mu},\\\
{\left.{\frac{{{\cal A}(\tau\to\mu\nu\bar{\nu})}}{{{\cal A}(\tau\to
e\nu\bar{\nu})}}}\right|_{{\rm{EXP}}}}&=1.0018(14){\mkern 1.0mu},\\\
{\left.{\frac{{{\cal A}(\tau\to e\nu\bar{\nu})}}{{{\cal A}(\mu\to
e\nu\bar{\nu})}}}\right|_{{\rm{EXP}}}}&=1.0010(14){\mkern 1.0mu},\end{split}$
(3)

with the correlations also given in Ref. Amhis  _et al._ (2019).

Furthermore, the effect in $\mathcal{A}(\mu\to
e\,\nu_{\mu}\overline{\nu}_{e})$ leads to a modification of the Fermi constant
which enters not only the electroweak (EW) precision observables but also the
determination of $V_{ud}$ from beta decays. Super-allowed beta decays provide
the most precise determination of $V_{ud}$ leading to Seng  _et al._
(2020)222Alternative determinations can be found in Refs. Czarnecki  _et al._
(2019); Shiells  _et al._ (2020). In addition, there is the possibility of
“new nuclear corrections” (NNCs) Seng  _et al._ (2019); Gorchtein (2019).
However, as this issue is debated, we will not consider them here for the sake
of argument (i.e. pointing out the potential NP implications).

$\displaystyle V_{us}^{\beta}=0.2280(6)\,.$ (4)

This value of $V_{us}^{\beta}$ can now be compared to $V_{us}$ from kaon Aoki
_et al._ (2020) and tau decays Amhis  _et al._ (2019)

$\displaystyle\begin{split}V_{us}^{K_{\mu
3}}&=0.22345(67)\,,\;\;\;V_{us}^{K_{e3}}=0.22320(61)\,,\\\ V_{us}^{K_{\mu
2}}&=0.22534(42)\,,\;\;\;V_{us}^{\tau}=0.2195(19)\,,\end{split}$ (5)

which are significantly lower. This is what constitutes the CAA. The tension
can be alleviated by the NP effect given by

$\displaystyle V_{us}^{\beta}$
$\displaystyle\equiv\sqrt{1\\!-\\!{{(V_{ud}^{\beta})}^{2}}\\!-\\!|{V_{ub}}{|^{2}}}$
$\displaystyle\simeq
V_{us}^{\rm{L}}\\!\left[{1-{{\left({\frac{{V_{ud}^{\rm{L}}}}{{V_{us}^{\rm{L}}}}}\right)}^{2}}\,\delta(\mu\to
e\nu\nu)}\right].$ (6)

where $V_{us(ud)}^{\rm L}$ is the value appearing in the CKM matrix. As
$G_{F}$ enters also the calculation of the EW gauge boson masses and $Z$ pole
observables, a global fit is necessary. Adding the determinations of the CKM
elements to the these standard EW observables (see e.g. Ref. Crivellin  _et
al._ (2020c) for details on our input and implementation) calculated by HEPfit
De Blas  _et al._ (2020) we find

$\delta(\mu\to e\nu\nu)=0.00065(15)\,.$ (7)

###  II.2 $\ell\to\ell^{\prime}\gamma$

The singly charged scalar generates $\ell\to\ell^{\prime}\gamma$ (see Fig. 1).
Using the results of Ref. Crivellin  _et al._ (2018) we obtain

$\displaystyle\mathrm{Br}[\mu\to
e\gamma]=\frac{m_{\mu}^{3}}{4\pi\Gamma_{\mu}}\left(|c_{L}^{{e}\mu}|^{2}+|c_{R}^{{e}\mu}|^{2}\right)\,,$
(8)

with $\Gamma_{\mu}$ being the total width of the muon and

$\displaystyle\begin{split}c_{L}^{{e}\mu}&=\frac{e\,\xi_{13}^{*}\,\xi_{23}}{384\,\pi^{2}}\frac{m_{e}}{m_{\phi}^{2}}\,,\;\;c_{R}^{{e}\mu}=\frac{e\,\xi_{13}^{*}\,\xi_{23}}{384\,\pi^{2}}\frac{m_{\mu}}{m_{\phi}^{2}}\,.\end{split}$
(9)

In what follows we will neglect the mass of the electron and thus
$c_{L}^{{e}\mu}$. Similarly, the expressions for $\tau\to\mu(e)\gamma$ can be
obtained by a straightforward exchange of indices. The current experimental
limits at $90\%$ C.L. are Bertl  _et al._ (2006); Aubert  _et al._ (2010);
Baldini  _et al._ (2016):

$\displaystyle\begin{split}{\rm Br}\\!\left[\mu\rightarrow e\gamma\right]&\leq
4.2\times 10^{-13}\,,\\\ {\rm Br}\\!\left[\tau\rightarrow\mu\gamma\right]&\leq
4.4\times 10^{-8}\,,\\\ {\rm Br}\\!\left[\tau\rightarrow e\gamma\right]&\leq
3.3\times 10^{-8}\,.\end{split}$

Note that in principle also contributions to anomalous magnetic moments of
charged leptons are generated. However, since the effect in our model is not
chirally enhanced, the effect is numerically small and can be safely
neglected. Interestingly, note that the $\phi^{\pm}$ interactions do not
generate electric dipole moments (EDMs) (disregarding very small quark and
neutrino effects already present in the SM) and therefore automatically agree
with the latest very stringent bound on the electron EDM from measurements of
Rb atoms Andreev  _et al._ (2018).

Figure 2: Diagram showing the Drell-Yan pair production of singly charged
scalars. Their decay necessarily give rise to a signal with an oppositely
charged lepton pair and missing transverse energy.

###  II.3 $\ell\to\ell^{\prime}\ell^{\prime(\prime)}\ell^{\prime(\prime)}$

The singly charged scalar contributes to three-body decays to charged leptons
at the loop level. Here the dominant contribution for sizable couplings $\xi$
is the box diagram shown in Fig. 1. For concreteness, we give the results for
$\tau\to 3e$ and $\tau\to\mu ee$, while the other decays can be obtained by an
appropriate exchange of the flavour indices:

$\displaystyle\mathrm{Br}[\tau\to e\mu\mu)=$
$\displaystyle\;\frac{m_{\tau}^{5}}{1536\,\pi^{3}\,\Gamma_{\tau}}\left|\frac{\xi_{12}^{*}\,\xi_{23}\left(\left|\xi_{12}^{2}\right|+\left|\xi_{23}^{2}\right|-\left|\xi_{13}^{2}\right|\right)}{64\,\pi^{2}\,m_{\phi}^{2}}\right|^{2}\,,$
$\displaystyle\mathrm{Br}[\tau\to eee]=$
$\displaystyle\;\frac{m_{\tau}^{5}}{768\,\pi^{3}\,\Gamma_{\tau}}\left|\frac{\xi_{12}^{*}\,\xi_{23}\left(\left|\xi_{12}^{2}\right|+\left|\xi_{13}^{2}\right|\right)}{64\,\pi^{2}\,m_{\phi}^{2}}\right|^{2}\,,$
(10)

where $\Gamma_{\tau}$ is the total decay width of the tau. Here we did not
include the small on- and off-shell photon contributions (they are given in
appendix A, together with our results for $\mu\to e$ conversion) and we did
not give the branching ratios for the decays involving more than one flavour
change (such as $\tau\to e\mu e$) which must be tiny in our model due to the
measured smallness of $\mu\to e\gamma$. The corresponding experimental bounds
(95% CL) are Bellgardt  _et al._ (1988); Hayasaka  _et al._ (2010); Lees  _et
al._ (2010); Aaij  _et al._ (2015c); Amhis  _et al._ (2019)

$\displaystyle\begin{aligned} \operatorname{Br}[\mu^{-}\rightarrow
e^{-}e^{+}e^{-}]&\leq 1.0\times 10^{-12}\,,\\\
\operatorname{Br}[\tau^{-}\rightarrow e^{-}e^{+}e^{-}]&\leq 1.4\times
10^{-8}\,,\\\ \operatorname{Br}[\tau^{-}\rightarrow e^{-}\mu^{+}\mu]&\leq
1.6\times 10^{-8}\,,\\\
\operatorname{Br}[\tau^{-}\rightarrow\mu^{-}e^{+}e^{-})&\leq 1.1\times
10^{-8}\,,\\\ \operatorname{Br}[\tau^{-}\rightarrow\mu^{-}\mu^{+}\mu^{-}]&\leq
1.1\times 10^{-8}\,.\end{aligned}$ (11)

###  II.4 LHC searches

The singly charged $SU(2)_{L}$ singlet scalar has the same quantum numbers as
the right-handed slepton in supersymmetry. Therefore, bounds from direct
searches for smuons and selectrons can be recast to set bounds on our model
Cao  _et al._ (2018); Alcaide  _et al._ (2018); Alcaide and Mileo (2020). The
dominant contribution is given by Drell-Yan pair production of $\phi^{\pm}$,
represented by the Feynman diagram in Fig. 2. We assume that interference with
the SM background (mostly $W^{+}W^{-}$ production in this case) can be
neglected in the limit of a large enough $m_{\phi}$ and a narrow $\phi^{\pm}$
width.

Figure 3: Recast ATLAS bounds on $m_{\phi}$ and $\rm
BR(\phi^{+}\to\ell^{+}\nu)$. The red (green) region is excluded by
$e^{+}e^{-}$ ($\mu^{+}\mu^{-}$) searches (see main text for details). The
dashed lines represent the projected exclusion reach for an integrated
luminosity of 3/ab at the HL-LHC.  Figure 4: Preferred regions at the
$1\sigma$ level in the $\delta(\tau\to\mu\nu\nu)$–$\delta(\mu\to e\nu\nu)$
plane together with the predictions for $\tau\to e\gamma$ (magenta), $\tau\to
e\mu\mu$ (black) and $|\xi_{12}^{2}|/m_{\phi}^{2}$ (blue) which can be
constrained from mono photon searches at future $e^{+}e^{-}$ colliders.

For the reinterpretation of bounds, we consider the most recent ATLAS analysis
Aad  _et al._ (2020) with 139/fb of data, searching for final states with an
oppositely charged lepton pair ($e^{+}e^{-}$ or $\mu^{+}\mu^{-}$) and missing
transverse energy. The search targets sleptons decaying into leptons and
neutralinos, which corresponds to our setup in the case of a vanishing
neutralino mass. The ATLAS bounds on the right-handed slepton mass in this
limit is $\approx 425$ GeV for both the $e^{+}e^{-}$ and $\mu^{+}\mu^{-}$
channels and for a 100% branching ratio of the slepton into the given channel.
To reinterpret this result, we have simulated the pair production cross-
section at leading order with MG5_aMC Alwall  _et al._ (2014) and rescaled it
with a constant K-factor, obtained by matching our values with the production
cross section to the one given by ATLAS (for a right-handed slepton mass of
500 GeV). A conservative error of $10\%$ has been added on the cross-section
to account for the differences in the simulation procedures.

Fig. 3 shows the bounds in the $\rm m_{\phi}$–$\rm Br(\phi^{\pm}\to
e^{\pm}(\mu^{\pm})\nu)$ plane extracted from the analysis of the $e^{+}e^{-}$
and $\mu^{+}\mu^{-}$ channels of ATLAS. The red (green) hatched region is
excluded by the $e^{+}e^{-}$ ($\mu^{+}\mu^{-}$) channel. The colored bands
indicate the change in the limit obtained by linearly varying the efficiency
calculated on the value of the ATLAS bound by $\pm 40\%$ between 200 GeV and
425 GeV for $m_{\phi}$. The solid line corresponds to the estimated limit
without taking into account the additional uncertainties discussed above. As
due to the anti-symmetry of the couplings the sum of the branching ratio to
muons and electrons can never be smaller than $1/2$, we can set a coupling
independent limit $\approx 200$ GeV on $m_{\phi}$.

###  II.5 Mono Photon Searches

We can also set a lower limit of $|\xi_{12,13}^{2}|/m_{\phi}^{2}$ from dark
matter (DM) searches at LEP with mono-photon signatures. Using the DELPHI
analysis of Refs. Abdallah  _et al._ (2005, 2009), Ref. Fox  _et al._ (2011)
was able to exploit the kinematic distributions to obtain a bound of $\approx
480$ GeV for zero DM mass on the DM mediator mass for unit coupling strength
and vectorial interactions (in the effective theory). Taking into account that
we have neutrinos and therefore interference with the SM, this translates into
a bound of $\approx 1\,$TeV.

Assuming that $m_{\phi}$ is sufficiently above the LEP production threshold,
as suggested by LHC searches discussed above, we can recast these results.
Taking into account that we have a left-handed vector current, we find
$(|\xi_{12,13}^{2}|)/m_{\phi}^{2}\lessapprox 1/(175\,{\rm GeV})^{2}$. This
bound would be strengthened if $\xi_{12}$ and $\xi_{13}$ are simultaneously
non zero, but further weakens as $m_{\phi}$ approaches the LEP beam energy.
Therefore, it is not yet competitive with flavour bounds but could be
significantly improved at future $e^{+}e^{-}$ colliders.

##  III Phenomenology

Let us start our phenomenology by considering the NP effect in
$\tau\to\mu\nu\nu$ and $\mu\to e\nu\nu$. The currently preferred regions (at
the $1\,\sigma$ level) for $\delta(\tau\to\mu\nu\nu)$ and $\delta(\mu\to
e\nu\nu)$ is shown in Fig. 4 as the orange and red regions, respectively,
while the combined region is shown in green. Note that for any point within
the combined region, $\xi_{13}$ must be vanishingly small in order not to
violate the bounds from $\mu\to e\gamma$ or $\mu\to e$ conversion. Therefore,
we can neglect its effect in the following.

This means that in this setup ($\xi_{13}\\!\simeq\\!0$) and we have ${\rm
Br}(\phi^{+}\to\mu^{+}\nu)=0.5$, which leads to a bound of $\approx 300$ GeV
from the $\mu^{+}\mu^{-}$ channel. This bound could be further improved at the
HL-LHC Apo (2017) (around 30% more as shown in Fig. 3, rescaling the ATLAS
bound for an integrated luminosity of 3/ab) or at the FCC-hh Abada  _et al._
(2019a) where, considering the projections for other scenarios in absence of a
dedicated analysis, we estimate a potential improvement up to a factor of few
Baumholzer  _et al._ (2020). Furthermore, we can correlate
$\delta(\tau\to\mu\nu\nu)$ and $\delta(\mu\to e\nu\nu)$ directly to $\tau\to
e\gamma$, as indicated by the magenta lines in Fig. 4, while ${\rm Br}[\tau\to
e\gamma]\approx 0$. The predicted branching ratio for $\tau\to e\gamma$ is of
the order of a few times $10^{-11}$. Furthermore, we can also obtain
correlations with $\tau\to 3e$ and $\tau\to e\mu\mu$. Since the branching
ratio of the latter is predicted to be larger (for the region preferred by
data), we depict it in Fig. 4 as black lines. However, here the correlations
is not direct since it depends on $m_{\phi}$ and we find ${\rm Br}[\tau\to
e\mu\mu]\approx 10^{-10}m_{\phi}^{4}/(5\,{\rm TeV})^{4}$. Interestingly, this
lies within the reach of BELLE II Inami (2016) or the FCC-ee Pich (2020). We
also depict constant values of $|\xi_{12}^{2}|/m_{\phi}^{2}$ as dashed blue
lines. Even though their values are significantly below the LEP bounds
discussed above, future $e^{+}e^{-}$ colliders like ILC Baer  _et al._ (2013),
CLIC Aicheler  _et al._ (2012), CEPC An  _et al._ (2019) or FCC-ee Abada  _et
al._ (2019b) could test the predicted mono-photon signature. In particular,
the ILC can improve the bound on the Wilson coefficient by a factor 50
Habermehl  _et al._ (2020), CEPC by a factor 40 Liu  _et al._ (2019) and even
bigger improvements could be expected CLIC, and FCC-ee for which a dedicated
study is stongly motivated.

##  IV Conclusions

In this article we studied the phenomenology of the singly charged $SU(2)_{L}$
singlet scalar in the light of the hints for LFU violation encoded in the
measured values for $\tau\to\mu\nu\nu/\tau(\mu)\to e\nu\nu$ and $V_{us}$ from
super allowed beta decays (using CKM unitarity) compared to $V_{us}$ from Kaon
and tau decays. Interestingly, the singly charged scalar has only three free
couplings (due to hemiticity of the Lagrangain) and is therefore very
predictive: it necessarily violates LF(U) and leads to a positive definite
effect in $\ell\to\ell^{\prime}\nu\nu$ as preferred by data while there is no
(pure) NP contribution to the otherwise so stringently constraining electron
EDM.

Recasting ATLAS searches for right-handed sleptons we derive a coupling
independent limit of $m_{\phi}\approx 200\,$GeV. In the region preferred by
LFU violation in tau decays and the CAA, $\xi_{13}\approx 0$ is required by
$\mu\to e\gamma$, leading to an LHC bound of $m_{\phi}\approx 300\,$GeV.
Concerning LFV, we predicted ${\rm Br}[\tau\to e\gamma]$ to be of order of a
few times $10^{-11}$ and ${\rm Br}[\tau\to e\mu\mu]\approx
10^{-10}m_{\phi}^{2}/(5\,{\rm TeV})^{2}$. Furthermore, we pointed out that our
model can be tested not only by future experiments searching for these LFV
decays, but also via direct searches at the HL (HE) LHC and FCC-hh and by mono
photon searches at future $e^{+}e^{-}$ colliders. In particular, FCC-hh could
improve the bound on $m_{\phi}$ and push the predicted value for ${\rm
Br}{[\tau\to e\mu\mu]}$ towards the region observable by BELLE II and FCC-ee,
providing a prime example of complementarity between low energy precision
experiments and direct searches for NP.

###### Acknowledgements.

We thank Joachim Kopp and Marc Montull for useful discussions. The work of
A.C., C.A.M and F.K. is supported by a Professorship Grant (PP00P2_176884) of
the Swiss National Science Foundation. A.C. thanks CERN for the support via
the scientific associate program. The work of LP is supported by the Knut and
Alice Wallenberg foundation under the SHIFT project, grant KAW 2017.0100.

##  Appendix A $\mu\to e$ conversion

The SM-contributions to $\mu\to e$ conversion can safely be neglected. We
parametrise the NP-contributions by the effective Lagrangian

$\displaystyle\mathcal{L}_{\text{eff}}=\sum_{q=u,d}\left(C_{qq}^{V,LL}O_{qq}^{V,LL}+C_{qq}^{V,LR}O_{qq}^{V,LR}\right)+\mathrm{h.c.}\,,$

with

$\displaystyle O_{qq}^{V,LL(R)}$
$\displaystyle=(\bar{e}\gamma^{\mu}P_{L}\mu)(\bar{q}\gamma_{\mu}P_{L(R)}q),$
(12)

The singly charged scalar contributes to this process via the off-shell photon
penguin. In this case the vectorial nature of the photon coupling leads to
$C_{qq}^{V,LL}=C_{qq}^{V,LR}$. At leading order we have:

$\displaystyle\begin{split}C_{qq}^{V,LL}&=\frac{e^{2}Q_{q}}{288\,\pi^{2}\,m_{\phi}^{2}}\xi_{13}^{*}\,\xi_{23}\,,\end{split}$
(13)

where $Q_{q}$ is the electric charge of the quarks
($Q_{u}=+\frac{2}{3},\;Q_{d}=-\frac{1}{3}$).
The transition rate $\Gamma_{\mu\to e}^{N}\equiv\Gamma(\mu N\to eN)$ is given
by

$\displaystyle\Gamma_{\mu\to e}^{N}=$
$\displaystyle\,4\,m_{\mu}^{5}\;\Bigg{|}\sum_{q=u,d}\left(C_{qq}^{V,RL}+C_{qq}^{V,RR}\right)$
$\displaystyle\times\left(f_{Vp}^{(q)}V_{N}^{p}\,+\,f_{Vn}^{(q)}V_{N}^{n}\right)\Bigg{|}^{2}+\left(L\leftrightarrow
R\right).$ (14)

The nucleon vector form factors are the same as the ones measured in elastic
electron-hadron scattering, i.e.

$\displaystyle
f_{Vp}^{(u)}=2,\;f_{Vn}^{(u)}=1,\;f_{Vp}^{(d)}=1,\;f_{Vn}^{(d)}=2,$ (15)

The overlap integrals $V_{p/n}^{N}$ depend on the nature of the target $N$. We
use the numerical values Kitano  _et al._ (2002)

$\displaystyle\begin{split}V_{\text{Au}}^{p}&=0.0974\,,\quad
V_{\text{Au}}^{n}=0.146\,.\end{split}$ (16)

The branching ratio of $\mu\to e$ conversion is defined as the transition rate
divided by the $\mu$ capture rate:

$\displaystyle\mathrm{Br}(\mu\to
e)={\Gamma^{\text{conv}}}/{\Gamma^{\text{capt}}}\,,$ (17)

and for the latter we use Suzuki  _et al._ (1987)

$\displaystyle\Gamma_{\rm Au}^{\rm capt}=8.7\times
10^{-18}\;\text{GeV}\,,\quad\Gamma_{\rm Al}^{\rm capt}=4.6\times
10^{-19}\;\text{GeV}\,.$ (18)

The experimental limit on $\mu\to e$ conversion isBertl  _et al._ (2006)

$\displaystyle\frac{\Gamma_{\rm Au}^{\text{conv}}}{\Gamma_{\rm
Au}^{\text{capt}}}$ $\displaystyle<7.0\times 10^{-13}\quad\text{SINDRUM
II}\,.$ (19)

Adding the on-shell (see Eq. 9) and off-shell (see Eq. 13) photon
contributions to the $\tau$-decays of Eq. 10, we obtain

$\displaystyle{\rm Br}(\tau\to 3e)=$
$\displaystyle\,\frac{e^{2}\,m_{\tau}^{3}}{192\,\pi^{3}\,\Gamma_{\tau}}\left|c_{R}^{e\tau}\right|^{2}\left(4\log\left[\frac{m_{\tau}^{2}}{m_{e}^{2}}\right]-11\right)$
$\displaystyle+\frac{m_{\tau}^{5}}{3072\,\pi^{3}\,\Gamma_{\tau}}\left(\left|\left(\left|\xi_{12}\right|^{2}+\left|\xi_{13}\right|^{2}\right)\frac{\xi_{12}^{*}\,\xi_{23}}{32\,\pi^{2}\,m_{\phi}^{2}}+\,\frac{e^{2}}{288\,\pi^{2}}\frac{\xi_{12}^{*}\,\xi_{23}}{m_{\phi}^{2}}\right|^{2}+\left|\frac{e^{2}}{288\,\pi^{2}}\frac{\xi_{12}^{*}\,\xi_{23}}{m_{\phi}^{2}}\right|^{2}\right)$
$\displaystyle+\frac{e\,m_{\tau}^{4}}{384\,\pi^{3}\,\Gamma_{\tau}}\left(2\left(\left|\xi_{12}\right|^{2}+\left|\xi_{13}\right|^{2}\right)\frac{{\rm
Re}\left(c_{R}^{e\tau}\,\xi_{12}\,\xi_{23}^{*}\right)}{32\,\pi^{2}\,m_{\phi}^{2}}+\frac{3e^{2}}{288\,\pi^{2}}\frac{{\rm
Re}\left(c_{R}^{e\tau}\,\xi_{12}\,\xi_{23}^{*}\right)}{m_{\phi}^{2}}\right)\,,$
(20) $\displaystyle{\rm Br}(\tau\to e^{\mp}\mu^{\pm}\mu^{\mp})=$
$\displaystyle\,\frac{e^{2}\,m_{\tau}^{3}}{48\,\pi^{3}\,\Gamma_{\tau}}\left|c_{R}^{e\tau}\right|^{2}\left(\log\left[\frac{m_{\tau}^{2}}{m_{\mu}^{2}}\right]-3\right)$
$\displaystyle+\frac{m_{\tau}^{5}}{3072\,\pi^{3}\Gamma_{\tau}}\left(\left|\left(\left|\xi_{12}\right|^{2}-\left|\xi_{13}\right|^{2}+\left|\xi_{23}\right|^{2}\right)\frac{\xi_{12}^{*}\,\xi_{23}\,}{64\,\pi^{2}\,m_{\phi}^{2}}+\frac{e^{2}}{288\,\pi^{2}}\frac{\xi_{12}^{*}\,\xi_{23}}{m_{\phi}^{2}}\right|^{2}+\left|\frac{e^{2}}{288\,\pi^{2}}\frac{\xi_{12}^{*}\,\xi_{23}}{m_{\phi}^{2}}\right|^{2}\right)$
$\displaystyle+\frac{e\,m_{\tau}^{4}}{384\,\pi^{5}\,\Gamma_{\tau}}\left(2\left(\left|\xi_{12}\right|^{2}-\left|\xi_{13}\right|^{2}+\left|\xi_{23}\right|^{2}\right)\frac{{\rm
Re}\left(c_{R}^{e\tau}\,\xi_{12}\,\xi_{23}^{*}\right)}{64\,\pi^{2}\,m_{\phi}^{2}}+\frac{3e^{2}}{288\,\pi^{2}}\frac{{\rm
Re}\left(c_{R}^{e\tau}\,\xi_{12}\,\xi_{23}^{*}\right)}{m_{\phi}^{2}}\right)\,.$

## References

  * Butler (2017) J. N. Butler (CMS), in  _5th Large Hadron Collider Physics Conference_ (2017) arXiv:1709.03006 [hep-ex] .
  * Masetti (2018) L. Masetti (ATLAS), Nucl. Part. Phys. Proc. 303-305, 43 (2018).
  * Aaij  _et al._ (2014a) R. Aaij _et al._ (LHCb), JHEP 06, 133 (2014a), arXiv:1403.8044 [hep-ex] .
  * Aaij  _et al._ (2014b) R. Aaij _et al._ (LHCb), Phys. Rev. Lett. 113, 151601 (2014b), arXiv:1406.6482 [hep-ex] .
  * Aaij  _et al._ (2015a) R. Aaij _et al._ (LHCb), JHEP 09, 179 (2015a), arXiv:1506.08777 [hep-ex] .
  * Aaij  _et al._ (2016) R. Aaij _et al._ (LHCb), JHEP 02, 104 (2016), arXiv:1512.04442 [hep-ex] .
  * Khachatryan  _et al._ (2016) V. Khachatryan _et al._ (CMS), Phys. Lett. B 753, 424 (2016), arXiv:1507.08126 [hep-ex] .
  * ATLAS-CONF-2017- (023) ATLAS-CONF-2017-023 (ATLAS), .
  * CMS-PAS-BPH-15- (008) CMS-PAS-BPH-15-008 (CMS), .
  * Aaij  _et al._ (2017) R. Aaij _et al._ (LHCb), JHEP 08, 055 (2017), arXiv:1705.05802 [hep-ex] .
  * Lees  _et al._ (2012) J. Lees _et al._ (BaBar), Phys. Rev. Lett. 109, 101802 (2012), arXiv:1205.5442 [hep-ex] .
  * Lees  _et al._ (2013) J. Lees _et al._ (BaBar), Phys. Rev. D 88, 072012 (2013), arXiv:1303.0571 [hep-ex] .
  * Aaij  _et al._ (2015b) R. Aaij _et al._ (LHCb), Phys. Rev. Lett. 115, 111803 (2015b), [Erratum: Phys.Rev.Lett. 115, 159901 (2015)], arXiv:1506.08614 [hep-ex] .
  * Aaij  _et al._ (2018a) R. Aaij _et al._ (LHCb), Phys. Rev. D 97, 072013 (2018a), arXiv:1711.02505 [hep-ex] .
  * Aaij  _et al._ (2018b) R. Aaij _et al._ (LHCb), Phys. Rev. Lett. 120, 171802 (2018b), arXiv:1708.08856 [hep-ex] .
  * Abdesselam  _et al._ (2019) A. Abdesselam _et al._ (Belle), (2019), arXiv:1904.08794 [hep-ex] .
  * Capdevila  _et al._ (2018) B. Capdevila, A. Crivellin, S. Descotes-Genon, J. Matias, and J. Virto, JHEP 01, 093 (2018), arXiv:1704.05340 [hep-ph] .
  * Altmannshofer  _et al._ (2017) W. Altmannshofer, P. Stangl, and D. M. Straub, Phys. Rev. D 96, 055008 (2017), arXiv:1704.05435 [hep-ph] .
  * D’Amico  _et al._ (2017) G. D’Amico, M. Nardecchia, P. Panci, F. Sannino, A. Strumia, R. Torre, and A. Urbano, JHEP 09, 010 (2017), arXiv:1704.05438 [hep-ph] .
  * Ciuchini  _et al._ (2017) M. Ciuchini, A. M. Coutinho, M. Fedele, E. Franco, A. Paul, L. Silvestrini, and M. Valli, Eur. Phys. J. C 77, 688 (2017), arXiv:1704.05447 [hep-ph] .
  * Hiller and Nisandzic (2017) G. Hiller and I. Nisandzic, Phys. Rev. D 96, 035003 (2017), arXiv:1704.05444 [hep-ph] .
  * Geng  _et al._ (2017) L.-S. Geng, B. Grinstein, S. Jäger, J. Martin Camalich, X.-L. Ren, and R.-X. Shi, Phys. Rev. D 96, 093006 (2017), arXiv:1704.05446 [hep-ph] .
  * Hurth  _et al._ (2017) T. Hurth, F. Mahmoudi, D. Martinez Santos, and S. Neshatpour, Phys. Rev. D 96, 095034 (2017), arXiv:1705.06274 [hep-ph] .
  * Alok  _et al._ (2017) A. K. Alok, B. Bhattacharya, A. Datta, D. Kumar, J. Kumar, and D. London, Phys. Rev. D 96, 095009 (2017), arXiv:1704.07397 [hep-ph] .
  * Algueró  _et al._ (2019) M. Algueró, B. Capdevila, A. Crivellin, S. Descotes-Genon, P. Masjuan, J. Matias, M. Novoa Brunet, and J. Virto, Eur. Phys. J. C 79, 714 (2019), [Addendum: Eur.Phys.J.C 80, 511 (2020)], arXiv:1903.09578 [hep-ph] .
  * Aebischer  _et al._ (2020) J. Aebischer, W. Altmannshofer, D. Guadagnoli, M. Reboud, P. Stangl, and D. M. Straub, Eur. Phys. J. C 80, 252 (2020), arXiv:1903.10434 [hep-ph] .
  * Ciuchini  _et al._ (2019) M. Ciuchini, A. M. Coutinho, M. Fedele, E. Franco, A. Paul, L. Silvestrini, and M. Valli, Eur. Phys. J. C 79, 719 (2019), arXiv:1903.09632 [hep-ph] .
  * Amhis  _et al._ (2017) Y. Amhis _et al._ (HFLAV), Eur. Phys. J. C 77, 895 (2017), arXiv:1612.07233 [hep-ex] .
  * Murgui  _et al._ (2019) C. Murgui, A. Peñuelas, M. Jung, and A. Pich, JHEP 09, 103 (2019), arXiv:1904.09311 [hep-ph] .
  * Shi  _et al._ (2019) R.-X. Shi, L.-S. Geng, B. Grinstein, S. Jäger, and J. Martin Camalich, JHEP 12, 065 (2019), arXiv:1905.08498 [hep-ph] .
  * Blanke  _et al._ (2019) M. Blanke, A. Crivellin, T. Kitahara, M. Moscati, U. Nierste, and I. Nišandžić,  (2019), 10.1103/PhysRevD.100.035035, [Addendum: Phys.Rev.D 100, 035035 (2019)], arXiv:1905.08253 [hep-ph] .
  * Kumbhakar  _et al._ (2020) S. Kumbhakar, A. K. Alok, D. Kumar, and S. U. Sankar, PoS EPS-HEP2019, 272 (2020), arXiv:1909.02840 [hep-ph] .
  * Bennett  _et al._ (2006) G. Bennett _et al._ (Muon g-2), Phys. Rev. D 73, 072003 (2006), arXiv:hep-ex/0602035 .
  * Aoyama  _et al._ (2020) T. Aoyama _et al._ , Phys. Rept. 887, 1 (2020), arXiv:2006.04822 [hep-ph] .
  * Belfatto  _et al._ (2020) B. Belfatto, R. Beradze, and Z. Berezhiani, Eur. Phys. J. C 80, 149 (2020), arXiv:1906.02714 [hep-ph] .
  * Grossman  _et al._ (2020) Y. Grossman, E. Passemar, and S. Schacht, JHEP 07, 068 (2020), arXiv:1911.07821 [hep-ph] .
  * Coutinho  _et al._ (2020) A. M. Coutinho, A. Crivellin, and C. A. Manzari, Phys. Rev. Lett. 125, 071802 (2020), arXiv:1912.08823 [hep-ph] .
  * Crivellin and Hoferichter (2020) A. Crivellin and M. Hoferichter, Phys. Rev. Lett. 125, 111801 (2020), arXiv:2002.07184 [hep-ph] .
  * Crivellin  _et al._ (2020a) A. Crivellin, F. Kirk, C. A. Manzari, and M. Montull, (2020a), arXiv:2008.01113 [hep-ph] .
  * Kirk (2020) M. Kirk, (2020), arXiv:2008.03261 [hep-ph] .
  * Alok  _et al._ (2020) A. K. Alok, A. Dighe, S. Gangal, and J. Kumar, (2020), arXiv:2010.12009 [hep-ph] .
  * Crivellin  _et al._ (2020b) A. Crivellin, C. A. Manzari, M. Alguero, and J. Matias, (2020b), arXiv:2010.14504 [hep-ph] .
  * Shiells  _et al._ (2020) K. Shiells, P. Blunden, and W. Melnitchouk, (2020), arXiv:2012.01580 [hep-ph] .
  * Seng  _et al._ (2020) C.-Y. Seng, X. Feng, M. Gorchtein, and L.-C. Jin, Phys. Rev. D 101, 111301 (2020), arXiv:2003.11264 [hep-ph] .
  * Amhis  _et al._ (2019) Y. S. Amhis _et al._ (HFLAV), (2019), arXiv:1909.12524 [hep-ex] .
  * Capdevila  _et al._ (2020) B. Capdevila, A. Crivellin, C. A. Manzari, and M. Montull, (2020), arXiv:2005.13542 [hep-ph] .
  * (47) A. Buras, A. Crivellin, F. Kirk, C. A. Manzari, and M. Montull, In preparation .
  * Zee (1986) A. Zee, Nucl. Phys. B 264, 99 (1986).
  * Babu (1988) K. Babu, Phys. Lett. B 203, 132 (1988).
  * Krauss  _et al._ (2003) L. M. Krauss, S. Nasri, and M. Trodden, Phys. Rev. D 67, 085002 (2003), arXiv:hep-ph/0210389 .
  * Nebot  _et al._ (2008) M. Nebot, J. F. Oliver, D. Palao, and A. Santamaria, Phys. Rev. D 77, 093013 (2008), arXiv:0711.0483 [hep-ph] .
  * Cai  _et al._ (2015) Y. Cai, J. D. Clarke, M. A. Schmidt, and R. R. Volkas, JHEP 02, 161 (2015), arXiv:1410.0689 [hep-ph] .
  * Cheung and Seto (2004) K. Cheung and O. Seto, Phys. Rev. D 69, 113009 (2004), arXiv:hep-ph/0403003 .
  * Ahriche  _et al._ (2014) A. Ahriche, S. Nasri, and R. Soualah, Phys. Rev. D 89, 095010 (2014), arXiv:1403.5694 [hep-ph] .
  * Chen  _et al._ (2014) C.-S. Chen, K. L. McDonald, and S. Nasri, Phys. Lett. B 734, 388 (2014), arXiv:1404.6033 [hep-ph] .
  * Ahriche  _et al._ (2016) A. Ahriche, K. L. McDonald, and S. Nasri, JHEP 02, 038 (2016), arXiv:1508.02607 [hep-ph] .
  * Herrero-Garcia  _et al._ (2014) J. Herrero-Garcia, M. Nebot, N. Rius, and A. Santamaria, Nucl. Phys. B 885, 542 (2014), arXiv:1402.4491 [hep-ph] .
  * Centelles Chuliá  _et al._ (2018) S. Centelles Chuliá, R. Srivastava, and J. W. Valle, Phys. Lett. B 781, 122 (2018), arXiv:1802.05722 [hep-ph] .
  * Babu  _et al._ (2020) K. Babu, P. B. Dev, S. Jana, and A. Thapa, JHEP 03, 006 (2020), arXiv:1907.09498 [hep-ph] .
  * Nieves and Pal (2004) J. F. Nieves and P. B. Pal, Am. J. Phys. 72, 1100 (2004), arXiv:hep-ph/0306087 .
  * Czarnecki  _et al._ (2019) A. Czarnecki, W. J. Marciano, and A. Sirlin, Phys. Rev. D 100, 073008 (2019), arXiv:1907.06737 [hep-ph] .
  * Seng  _et al._ (2019) C. Y. Seng, M. Gorchtein, and M. J. Ramsey-Musolf, Phys. Rev. D 100, 013001 (2019), arXiv:1812.03352 [nucl-th] .
  * Gorchtein (2019) M. Gorchtein, Phys. Rev. Lett. 123, 042503 (2019), arXiv:1812.04229 [nucl-th] .
  * Aoki  _et al._ (2020) S. Aoki _et al._ (Flavour Lattice Averaging Group), Eur. Phys. J. C 80, 113 (2020), arXiv:1902.08191 [hep-lat] .
  * Crivellin  _et al._ (2020c) A. Crivellin, M. Hoferichter, C. A. Manzari, and M. Montull, Phys. Rev. Lett. 125, 091801 (2020c), arXiv:2003.04886 [hep-ph] .
  * De Blas  _et al._ (2020) J. De Blas _et al._ , Eur. Phys. J. C 80, 456 (2020), arXiv:1910.14012 [hep-ph] .
  * Crivellin  _et al._ (2018) A. Crivellin, M. Hoferichter, and P. Schmidt-Wellenburg, Phys. Rev. D 98, 113002 (2018), arXiv:1807.11484 [hep-ph] .
  * Bertl  _et al._ (2006) W. H. Bertl _et al._ (SINDRUM II), Eur. Phys. J. C 47, 337 (2006).
  * Aubert  _et al._ (2010) B. Aubert _et al._ (BaBar), Phys. Rev. Lett. 104, 021802 (2010), arXiv:0908.2381 [hep-ex] .
  * Baldini  _et al._ (2016) A. Baldini _et al._ (MEG), Eur. Phys. J. C 76, 434 (2016), arXiv:1605.05081 [hep-ex] .
  * Andreev  _et al._ (2018) V. Andreev _et al._ (ACME), Nature 562, 355 (2018).
  * Bellgardt  _et al._ (1988) U. Bellgardt _et al._ (SINDRUM), Nucl. Phys. B 299, 1 (1988).
  * Hayasaka  _et al._ (2010) K. Hayasaka _et al._ , Phys. Lett. B 687, 139 (2010), arXiv:1001.3221 [hep-ex] .
  * Lees  _et al._ (2010) J. Lees _et al._ (BaBar), Phys. Rev. D 81, 111101 (2010), arXiv:1002.4550 [hep-ex] .
  * Aaij  _et al._ (2015c) R. Aaij _et al._ (LHCb), JHEP 02, 121 (2015c), arXiv:1409.8548 [hep-ex] .
  * Cao  _et al._ (2018) Q.-H. Cao, G. Li, K.-P. Xie, and J. Zhang, Phys. Rev. D 97, 115036 (2018), arXiv:1711.02113 [hep-ph] .
  * Alcaide  _et al._ (2018) J. Alcaide, M. Chala, and A. Santamaria, Phys. Lett. B 779, 107 (2018), arXiv:1710.05885 [hep-ph] .
  * Alcaide and Mileo (2020) J. Alcaide and N. I. Mileo, Phys. Rev. D 102, 075030 (2020), arXiv:1906.08685 [hep-ph] .
  * Aad  _et al._ (2020) G. Aad _et al._ (ATLAS), Eur. Phys. J. C 80, 123 (2020), arXiv:1908.08215 [hep-ex] .
  * Alwall  _et al._ (2014) J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli, and M. Zaro, JHEP 07, 079 (2014), arXiv:1405.0301 [hep-ph] .
  * Abdallah  _et al._ (2005) J. Abdallah _et al._ (DELPHI), Eur. Phys. J. C 38, 395 (2005), arXiv:hep-ex/0406019 .
  * Abdallah  _et al._ (2009) J. Abdallah _et al._ (DELPHI), Eur. Phys. J. C 60, 17 (2009), arXiv:0901.4486 [hep-ex] .
  * Fox  _et al._ (2011) P. J. Fox, R. Harnik, J. Kopp, and Y. Tsai, Phys. Rev. D 84, 014028 (2011), arXiv:1103.0240 [hep-ph] .
  * Apo (2017) 4/2017 (2017), 10.23731/CYRM-2017-004.
  * Abada  _et al._ (2019a) A. Abada _et al._ (FCC), Eur. Phys. J. ST 228, 755 (2019a).
  * Baumholzer  _et al._ (2020) S. Baumholzer, V. Brdar, P. Schwaller, and A. Segner, JHEP 09, 136 (2020), arXiv:1912.08215 [hep-ph] .
  * Inami (2016) K. Inami (Belle-II), PoS ICHEP2016, 574 (2016).
  * Pich (2020) A. Pich, (2020), arXiv:2012.07099 [hep-ph] .
  * Baer  _et al._ (2013) H. Baer _et al._ , (2013), arXiv:1306.6352 [hep-ph] .
  * Aicheler  _et al._ (2012) M. Aicheler _et al._ ,  (2012), 10.5170/CERN-2012-007.
  * An  _et al._ (2019) F. An _et al._ , Chin. Phys. C 43, 043002 (2019), arXiv:1810.09037 [hep-ex] .
  * Abada  _et al._ (2019b) A. Abada _et al._ (FCC), Eur. Phys. J. ST 228, 261 (2019b).
  * Habermehl  _et al._ (2020) M. Habermehl, M. Berggren, and J. List, Phys. Rev. D 101, 075053 (2020), arXiv:2001.03011 [hep-ex] .
  * Liu  _et al._ (2019) Z. Liu, Y.-H. Xu, and Y. Zhang, JHEP 06, 009 (2019), arXiv:1903.12114 [hep-ph] .
  * Kitano  _et al._ (2002) R. Kitano, M. Koike, and Y. Okada, Phys. Rev. D 66, 096002 (2002), [Erratum: Phys.Rev.D 76, 059902 (2007)], arXiv:hep-ph/0203110 .
  * Suzuki  _et al._ (1987) T. Suzuki, D. F. Measday, and J. Roalsvig, Phys. Rev. C 35, 2212 (1987).


11institutetext: GEPI, Observatoire de Paris, PSL University, CNRS, 5 Place
Jules Janssen, 92190 Meudon, France
11email<EMAIL_ADDRESS>22institutetext: ACRI-ST, 260 route du Pin
Montard, 06904, Sophia Antipolis, France 33institutetext: Univ. Grenoble
Alpes, CNRS, IPAG, 38000 Grenoble, France

# Updated Gaia-2MASS 3D maps of Galactic interstellar dust

R. Lallement 11 J.L. Vergely 22 C. Babusiaux 33 N.L.J. Cox 22

(Received ; accepted )

###### Abstract

Aims. Three-dimensional (3D) maps of Galactic interstellar dust are a tool for
a wide range of uses. We aim to construct 3D maps of dust extinction in the
Local Arm and surrounding regions.

Methods. To do this, Gaia EDR3 photometric data were combined with 2MASS
measurements to derive extinction towards stars with accurate photometry and
relative uncertainties on EDR3 parallaxes of less than 20%. We applied our
hierarchical inversion algorithm adapted to inhomogeneous spatial
distributions of target stars to this catalogue of individual extinctions.

Results. We present the updated 3D dust extinction distribution and provide an
estimate of the error on integrated extinctions from the Sun to each area in
the 3D map. The full computational area is similar to the one of the previous
DR2 map, that is to say with a 6 $x$ 6 $x$ 0.8 kpc3 volume around the Sun. Due
to the addition of fainter target stars, the volume in which the clouds can be
reconstructed has increased. Due to the improved accuracy of the parallaxes
and photometric data in EDR3, extinctions among neighbouring targets are more
consistent, allowing one to reach an increased contrast in the dense areas,
while cavity contours are more regular. We show several comparisons with
recent results on dust and star distributions. The wavy pattern around the
Plane of the dust concentrations is better seen and exists over large regions.
Its mean vertical peak-to-peak amplitude is of the order of 300 pc;
interestingly, it is similar to the vertical period of the spectacular snail-
shaped stellar kinematical pattern discovered in Gaia data.

Conclusions. The Gaia EDR3 catalogue allows for a significant improvement of
the extinction maps to be made, both in extent and quality. The hierarchical
technique confirms its efficiency in the inversion of massive datasets. Future
comparisons between 3D maps of interstellar matter and stellar distributions
may help to understand which mergers or internal perturbations have shaped the
Galaxy within the first 3 kpc.

###### Key Words.:

Dust: extinction; ISM: lines and bands; ISM: structure ; ISM: solar
neighborhood ; ISM: Galaxy

##  1 Introduction

Three-dimensional (3D) maps of extinction by the Galactic interstellar (IS)
dust are a useful tool in many respects. Among the main applications, they
allow de-reddening of spectra to provide intrinsic properties of observed
objects, they provide information on the foreground, environment, or
background for specific Galactic targets, and they can be used in conjunction
with stellar populations and gas distributions to infer the past history of
the Milky Way. In the context of interstellar matter studies, they can be used
in conjunction with other absorption data and/or multi-wavelength interstellar
emission measurements, including polarised emission, and they provide a frame
for models of light or particle propagation. Three-dimensional dust maps are
constructed through the synthesis of estimated starlight extinctions and they
have started to develop rapidly during the last ten years with the advent of
massive stellar surveys. Their development is considerably boosted by the Gaia
mission data releases.

Three-dimensional dust extinction mapping is based on the tomographic
inversion of measured extinction-distance pairs or extinction-distance
probability distributions for large numbers of targets distributed in 3D
space. What differs among the various 3D reconstructions is the choice of the
photometric and(or) spectro-photometric catalogues used to estimate individual
extinctions, the source of distances, that is parallax, photometry, or
combined parallax and photometry, and finally the inversion technique. Such
choices result in various extents of the mapped volumes. For descriptions of
recently constructed maps, readers can refer to Sale & Magorrian (2018); Green
et al. (2019); Lallement et al. (2019); Chen et al. (2019); Rezaei Kh. et al.
(2020); Chen et al. (2020); Leike et al. (2020); Hottier et al. (2020) and Guo
et al. (2021).

One of the main difficulties of 3D IS dust extinction tomography is the multi-
scale structure of the interstellar medium, from small dense cores to extended
warm clouds. Reconstructing the cores evidently requires a high concentration
of targets in their areas. This is also connected to the second difficulty,
the lack of homogeneity of the target spatial distribution. In general, the
number density of targets decreases with increasing distance due to
observational limitations. To adapt to the multi-scale structure, several
inversions allowed for a large spectrum of cloud sizes. The corresponding,
detailed maps were computed based on Gaussian random processes or techniques
derived from information field theory (Sale & Magorrian, 2018; Rezaei Kh. et
al., 2020; Leike et al., 2020). The present work uses a different technique
based on iterations of inversions at increasing resolutions. This technique
has been previously used to compute a dust map based on Gaia Data Release
(DR2) and 2MASS photometric data combined with Gaia DR2 parallaxes. The basic
principles of the 3D inversion have been detailed in Vergely et al. (2001) and
the new hierarchical technique has been described in Lallement et al. (2019).
Along this iterative process, the inversion is restricted to regions of 3D
space with a target density high enough to achieve the current resolution,
that is their total volume is gradually decreasing and there is a simultaneous
increase in the volumes where the retained solution is the solution of the
previous inversion. This allowed us to adapt the spatial resolution to the
concentration of the targets at each location. The advantage of this
hierarchical technique is the limited cost in computational time and the
possibility to increase the physical dimensions of the maps while avoiding
creating unreal structures due to the scarcity of targets.

The Gaia Early Data Release 3 (EDR3, Gaia Collaboration et al., 2021) contains
a new catalogue of parallaxes and luminosity measurements, with significant
improvements with respect to Gaia DR2, in particular the accuracy of
parallaxes has significantly increased and the homogeneity of photometric
measurements has thoroughly improved. We made use of the new measurements to
produce a larger and higher quality catalogue of individual extinctions to
enter the inversion process.

Section 2 describes the selection of EDR3 data and the extinction estimate
method. Section 3 recalls the main characteristics of the hierarchical
tomographic inversion and how it will be adapted to the new input data. It
contains the description of error estimates and potential biases. Section 4
contains comparisons between the map and locations of molecular clouds based
on different data, as well as a comparison between the dust distribution and
young stars close to the mid-plane. It describes several features of the new
3D map. Section 5 contains a summary and discusses mapping perspectives.

Figure 1: A0 extinction distribution (mag) in low extinction selections. In
blue, the local bubble selection is provided ($\varpi>14$ mas, corresponding
to the first 70 pc); and in green, a high latitude selection is shown
($|b|>80^{\circ}$, $|Z|<400$pc). Figure 2: De-reddened Hess diagram. The
greyscale corresponds to the square root of the stellar density. The red clump
is the most prominent feature whose shape can be compared to Fig. 2 of Ruiz-
Dern et al. (2018).

Figure 3: Extinction densities in the mid-plane for four scales: 200pc, 100pc,
50pc, and 25pc Figure 4: Differential extinction (or extinction density)
distribution, here in the plane parallel to the Galactic Plane and containing
the Sun. Units are in parsecs. The Sun is at (0,0) and the Galactic Centre
direction is to the right. The green isocontours surround the densest areas at
the level of 0.01 $mag\cdot pc^{-1}$. Other contours correspond to 0.00018,
0.0005, and 0.001 $mag\cdot pc^{-1}$. Figure 5: Comparison between the new
eDR3 inverted extinction density and the previous result using Gaia DR2. Here
we focus on the fraction of the 3D map containing the densest regions of the
Local Arm, the so-called ’split’, and Sagittarius. Left: Gaia DR2 map
(Lallement et al. (2019). Right: Gaia EDR3 map. The green contours mark the
densest areas, with the same level of 0.01 $mag\cdot pc^{-1}$ in the two maps.
We note the increased fraction of areas above this value, and, on the
contrary, the disappearance of faint artefacts within the cavities, which are
both due to reduced uncertainties on distances and estimated extinctions. The
local cavity around the Sun is better defined thanks to the addition of
nearby, bright targets, and it appears narrower. Figure 6: Relative error of
the extinction density in the Galactic Plane (top) and in a perpendicular cut
(bottom). Figure 7: Relative error on the integrated extinction in the
Galactic Plane (top) and in a perpendicular cut (bottom).

##  2 Extinction estimates

To derive the extinction of individual stars, we combined the Gaia EDR3 (Gaia
Collaboration et al., 2021) and the 2MASS (Skrutskie et al., 2006) photometry.
We used the catalogue’s cross-match provided within the Gaia archive. We
selected sources with a parallax uncertainty lower than 20%, astrometric
renormalised weight error $\texttt{\small{ruwe}}<1.4$, photometric
uncertainties lower than 2% in $G$, 5% in $G_{BP}$ and $G_{RP}$, $G_{BP}<20.3$
to avoid biases at low $G_{BP}$ fluxes (Riello et al., 2021), and $G>5$ to
avoid saturation. The corrected flux excess was used to remove sources
affected by crowding following the recipe of Riello et al. (2021) at
10$\sigma$. The $G$ photometry for the six-parameter solutions was corrected
following Riello et al. (2021). The parallax was corrected from the zero point
variation with magnitude, colour, and ecliptic latitude following Lindegren et
al. (2021). For 2MASS, we selected stars with 2MASS photometric quality flag
AAA and photometric uncertainties lower than 0.05 mag. Only stars with an
absolute magnitude $M_{G}<7.5$ were selected to avoid the near-infrared
discontinuity of the colour-colour relation for dwarfs redder than $G-K>3$. A
first selection on this criteria was done at the query level assuming no
extinction.

Similarly to Lallement et al. (2019), we derived intrinsic colour-colour
relations using low extinction stars. Here, we selected stars having
$A_{0}<0.05$ mag using the 3D extinction map of Lallement et al. (2019) (here
and throughout the paper, the subscript 0 means the monochromatic value at
5500Å). The main difference between our study and that of Lallement et al.
(2019) is that we went to smaller absolute magnitudes so that we had to derive
to separate colour-colour relations: one for the top of the HR diagram
($M_{G}<4$) and one for the dwarf sequence, with the hot dwarfs being part of
both calibrations, ensuring the continuity. For the top of the HR diagram, we
selected only stars within 500 pc of the Galactic Plane to avoid being
dominated by metal-poor giants. A 7 degree monotonic polynomial111R package
MonoPoly was adjusted to each colour-colour relation $G-X$ as a function of
$G-K$, removing strong outliers (ten times the median absolute deviation (MAD)
of the residuals) one by one. The colour range of our calibration is
$-1.0<(G-K)_{0}<3.7$ for the top of the HR diagram and $-1.0<(G-K)_{0}<3.0$
for the dwarfs. It is important to note that the colour range is smaller than
in Lallement et al. (2019) due to the selection of stars near the Galactic
Plane.

The extinction coefficients used are provided in the Gaia auxiliary data
pages222https://www.cosmos.esa.int/web/gaia/edr3-extinction-law. They have
been derived similarly to those used in Lallement et al. (2019) with the
update of the Gaia passbands from Riello et al. (2021) and the extinction law
of Fitzpatrick et al. (2019).

Similarly to Lallement et al. (2019), for each star the extinction $A_{0}$ and
the intrinsic colour $(G-K_{s})_{0}$ and their associated uncertainties were
determined through a maximum likelihood estimation333R package bbmle,

$\mathcal{L}=\prod_{X}P(G-X|A_{0},(G-K)_{0})$ (1)

for $X=G_{BP},G_{RP},J,H$. To avoid local minima, the following three
different initial values were tested: $(G-K)_{0}=G-K$ (no extinction), as well
as $(G-K)_{0}=1.5$ and $(G-K)_{0}=(G-K)_{\mathrm{max}}$ (the maximum colour of
the intrinsic colour relation). Both the top of the HR diagram and the dwarf
sequence calibrations were tested. We used the intrinsic colour-colour
relations to derive $(G-X)_{0}$ from $(G-K)_{0}$. Then
$P(G-X|A_{0},(G-K)_{0})=P(G-X|(G-X)_{0}+k[A_{0},(G-K)_{0}]A_{0})$. We modelled
this probability by a Gaussian, quadratically adding the photometric error in
the $X$ band and the intrinsic scatter of the colour-colour relation. Negative
values of $A_{0}$ were allowed to ensure a Gaussian uncertainty model, which
is needed for the inversion method. However no extrapolation of the extinction
coefficients was carried out, the negative values of $A_{0}$ being replaced by
0 to derive the extinction coefficients. Typical resulting uncertainties are
similar to Lallement et al. (2019), that is 0.3 mag in $A_{0}$ and 0.2 mag in
$(G-K)_{0}$. A chi-square test was performed to check the validity of the
resulting parameters, removing stars with a p-value limit smaller than 0.05.
We removed stars with uncertainties on the derived $A_{0}$ and $(G-K)_{0}$
higher than 0.5 and 0.4 mag, respectively, as well as outliers for which
$A_{0}$ is not compatible with being positive at 3 sigma. We also removed
stars for which the absolute magnitude did not correspond to the colour-colour
calibration data used at 3$\sigma$ and outside the calibration range, that is
$(G-K)_{0}>3$.

In contrast to Lallement et al. (2019), the high latitude stars present a
small bias towards negative values (Fig. 1). This is due to the fact that we
restricted our calibration sample to the rather metal-rich population of the
Galactic Plane, leading to a bias in the metal-poor-dominated population of a
high latitude sample. In the de-reddened HRD (Fig. 2), we see that the issue
affecting the blue stars has been resolved with EDR3. However a population of
stars on the left of the red giant branch is present, which should not be
there. This is due to a degeneracy between two different solutions for the
maximum likelihood. We suspect that the extinction coefficients we used do not
accurately reproduce the behaviour for the red giants in the near-infrared
bands. In future work, other spectral energy distributions will be tested to
derive the extinction coefficients. We note that those stars were too faint in
the $G_{BP}$ band to be in the Lallement et al. (2019) sample. Also, it is
important to note that those potential biases introduced by these objects are
taken into account in the errors (see next section).

More than 35 million location-extinction pairs were obtained from this
procedure and may enter the inversion. This is $\simeq$30% more than for the
catalogue based on DR2, and it allowed us to obtain better constraints on the
sources of extinction.

Figure 8: Online tool for extinction and extinction density radial profile: an
example. Figure 9: Online tool: example of extinction-density image. Here the
normal to the plane of the image is defined by its Galactic coordinates
(0,-90) and the origin of the image is defined by its distance from the Sun
r=100pc and its Galactic coordinates (0,90). As a result, we show the image of
the dust in a plane parallel to the mid-plane and it is displaced by 100 pc to
the south. We note that in this example we attempted an inversion up to five
kpc.

##  3 Hierarchical inversion of individual extinctions and error estimates

The 3D tomography algorithm used for inverting the new Gaia EDR3-2MASS
extinction catalogue is similar to the one used for inverting Gaia DR2-2MASS
data (Vergely et al., 2001; Lallement et al., 2019). We estimated the
extinction per unit distance $\rho$ = A0 / pc (or equivalently the extinction
density, a quantity directly related to the dust volume density) in all points
in 3D space based on the series of distance-extinction pairs for stars
spatially distributed around the Sun. The multi-scale approach we used
consists in estimating the extinction density in an iterative way, starting
with large spatial scales and moving towards decreasing scales (or increasing
spatial resolution). This method makes it possible to take account of the
heterogeneity of the distribution of the lines of sight. Indeed, some areas of
the sky are very well covered by data, while others are not; the second
situation happens beyond a certain distance or behind dense obscuring clouds.
For each spatial scale, two consecutive types of computations were performed.
First, the space was divided into quasi-cubic volumes of a linear size equal
to the considered scale divided by two and an estimate of the mean extinction
was computed in each of those volumes. Inversions were performed in radial
directions and radial extinction density profiles were then discretised and
the resulting data points were accumulated into an extinction density cube
which was updated direction by direction. Second, the resulting set of average
extinctions, associated with the corresponding average distances entered a 3D
inversion, using a correlation length equal to the current spatial scale. If,
in some of the individual volumes entering the 3D inversion, the number of
contributing stars was below a scale-dependent threshold, then this volume was
ignored, that is, we did not estimate the extinction density in the vicinity
of the barycentre of this volume and kept the value obtained at this location
during the previous iteration. If we had not proceeded in this way, some
spatially spread structures of extinction which have no physical reality could
have appeared in the vicinity of strong extinction areas (see for example the
results obtained without the hierarchical process and shown in Capitanio et
al. (2017)). Moreover, in such a way, the estimated density structures are
spatially more isotropic and the ’fingers’ of God’ effects are reduced.

Similarly to the algorithm used in Lallement et al. (2019), outliers were
detected by means of comparisons with neighbouring stars. Here a 3$\sigma$
selection was applied by comparing individual extinctions with the median of
the extinction at each given scale. This means that the outlier detection
depends on the scale investigated. At large spatial scales, some high
extinction stars may be removed if the majority of the stars are affected by
low extinction. However, the same high extinction stars may not be considered
as outliers at smaller spatial scales. This strategy avoids the extinction
density at large scales to be overestimated as was the case in the past in the
context of the first 3D map computations (Vergely et al. 2001). Each 3D
inversion is a Bayesian inversion which, for a given scale, uses the 3D
solution estimated at the previous scale as a prior. The extinction density
model was thus updated at each scale by considering data averaged over smaller
and smaller volumes, until a resolution of 25pc was reached (i.e. a sampling
of 10pc). Due to the target density threshold mentioned above, the achieved
resolution is space-dependent and is better (resp. poorer) in regions of high
(resp. low) target density. In order to illustrate the process, we show the
results of iterative inversions in Fig. 3 for the 200pc, 100pc, 50pc, and 25pc
scales. What is represented is the computed extinction density distribution in
the mid-plane. We note that by mid-plane, we mean the plane containing the Sun
which is parallel to the true Galactic mid-plane, that is, we neglected the
small distance between the Sun and the true Galactic Plane. In all figures we
use Cartesian coordinates with the Sun at centre X,Y,Z=0,0,0. The X axis is
directed to the Galactic Centre, the Y axis is along the direction of
rotation, and the Z axis points to the Northern Galactic Pole. Since the data
coverage decreases with distance, for the consecutive scales of 200pc, 100pc,
50pc, and 25pc, the computations were done in cubes of $20\,{\rm
kpc}\,\times\,20\,{\rm kpc}\,\times\,0.8\,{\rm kpc}$, $20\,{\rm
kpc}\,\times\,20\,{\rm kpc}\,\times\,0.8\,{\rm kpc}$, $10\,{\rm
kpc}\,\times\,10\,{\rm kpc}\,\times\,0.8\,{\rm kpc}$, and $6\,{\rm
kpc}\,\times\,6\,{\rm kpc}\,\times\,0.8\,{\rm kpc}$, respectively. However, we
show the results in the smallest, final mid-plane, that is a $6\,{\rm
kpc}\,\times\,6\,{\rm kpc}$ plane. The computation time to achieve this 25pc
resolution map was about 5 days using a four cores processor.

A frequent request from users of 3D dust maps is the availability of
accompanying distributions of errors, most of the time on integrated
extinctions or more rarely on extinction densities. As we have warned in
previous presentations of maps, our use of a correlation kernel prevents us
from reaching precision on the latter quantity because the extinction is
forced to be distributed in volumes that can be wider than the actual clouds.
On the other hand, it is important to note that errors on integrated
extinctions, the measured quantities, do not suffer from the imposed finite
resolution as much and they are expected to be much smaller. Here the error
budget could not be computed by means of a classical approach, that is
Gaussian error propagation. The main contribution to the error comes from
uncertainties on both extinction and distance; however, the target
distribution is an important factor too, and the correlation length has an
impact on extinction density errors. We have preferred a statistical approach
and used errors on target extinctions and distances to build randomly
distributed datasets of extinctions and distances around their observed
values. Such noisy data then entered the hierarchical approach separately
which we previously described as providing extinction densities. In such a
way, the actual spatial distribution of the targets was taken into account in
addition to distance and extinction uncertainties. Here we proceeded in this
way five times and we therefore built five cubes of extinction densities. The
standard deviation for these five realisations gives an idea of the errors on
the extinction density. In parallel, errors on the integrated extinction can
also be extracted by calculating, for the five cubes, the extinction integrals
from the Sun to any point and by calculating the standard deviation of these
integrated extinctions. We note that this error estimate may be affected by
the presence of remaining outliers or by the absence of data in certain areas.
Moreover, beyond a few kiloparsecs, the computed error is underestimated
because the model is drawn towards the prior which is constant in the absence
of data. Nevertheless, the order of magnitude of the error may be inferred
from the results at a shorter distance. In order to be more conservative, we
computed the errors in the case of the 50 pc kernel.

##  4 Results

###  4.1 Mapping improvements

Fig. 4 displays the new EDR3 extinction density distributions in the X,Y
plane. The Galactic Centre direction (or X axis) is to the right and the
direction of rotation (Y axis) is towards the top of the figure. The size of
each pixel is 5pc x 5pc. Units are $mag\cdot pc^{-1}$. Several iso-contours
have been drawn to help visualise both dense and faint structures. Compared
with the Gaia DR2-2MASS inverted distribution, there are additional dust
structures at large distances that were not recovered previously. This
extension is more visible in the first and fourth quadrants due to the
existence of numerous clouds from the Sagittarius and Scutum arms, but it is
also effective in the two other quadrants, for example we can now detect the
prolongation in the second quadrant of the Perseus clouds from the third
quadrant. Although distant structures are mapped at lower resolution, as
illustrated in Fig. 3, the increased number of targets allows for the 3 kpc x
3 kpc x 0.8 kpc volume for the large structures to be covered entirely.

The second main evolution is the increased contrast, that is the larger number
of peak densities and, conversely, the increased vacuity of the large
cavities. To illustrate this, isocontours for the same relatively high value
of differential opacity, namely 10-2 $mag\cdot pc^{-1}$, are drawn on both the
previous and new map in Fig. 5 in a smaller region around the Sun, and they
help to visualise the increased number of opaque areas in the new map. Because
our minimum kernel is similar in the two maps (25pc), this trend is a sign of
better consistency within the dataset, and it reflects the improved accuracy
of parallax distances and extinction estimates. Conversely, a careful look at
the differences between the two images shows that most of the cavities appear
to be emptier, and elongated structures within them that were due to a lack of
consistency of the data or distances that were too uncertain have disappeared
in the new map. This is also consistent with the increase in contrast and a
sign of a more reliable inversion. In particular, the ’avenue’ devoid of dust
between the ’split’ and the inner border of Sagittarius looks even more
conspicuous and extends up to 2.5 kpc from the Sun in the second quadrant.
Finally, we must draw attention to the local cavity (or the 100-200 pc low
density, wide area around the Sun). Due to the addition of very nearby target
stars, the cavity is more accurately mapped and appears narrower, as can be
seen in Fig. 5.

###  4.2 Error estimates

Orders of magnitudes of errors on the extinction densities (or dust volume
densities) are shown in Fig. 6. The performed computations show that the
average uncertainties on their amplitudes are around 30 or 40% in the main
clouds for a 50 pc kernel (Fig. 6). In the case of the 25 pc kernel, errors
should in principle be smaller; however, we have kept these 50 pc estimates to
be more conservative. Such large values are not unexpected, as explained in
section 3. They confirm that the maps give an order of magnitude of the dust
concentration, but accurate values within the dense clouds are out of reach.
We note that, because part of the uncertainties on the extinction density
arise from the fact that there are uncertainties on the positions of the
clouds (due to parallax errors), the standard deviation of different error
realisations increases in regions with a high density, characterised by small
clouds.

Errors on the integrated extinctions are shown in Fig. 7. They are of the
order of a few percent in most regions, that is to say they are much smaller
than the errors on the extinction densities. It may appear surprising since an
integrated extinction is a sum of extinctions’ densities, but the classical
scheme of the square root of the quadratic sum of errors weighted by distances
does not apply here since the integrated extinctions are the constrained
quantities.

###  4.3 Online tools

The new EDR3 extinction density distribution and corresponding error can be
downloaded under the form of a matrix in Cartesian Galactic coordinates
centred on the Sun at the EXPLORE website444 https://explore-platform.eu.
Several online tools are additionally available. The first tool is the
computation of cumulative extinction radial profiles for any direction entered
in Galactic or equatorial coordinate (ICRS) systems. Both figures and
corresponding numerical files can be downloaded. A series of coordinates can
also be entered in a single file. For users interested in the distribution of
clouds as a function of distance along a line of sight, the computation of the
extinction density can be computed and downloaded too. For both tools, the
step is 5pc. Fig. 8 shows an example of both outputs.

A second tool is the computation of the image of the extinction density
distribution in any plane that does or does not contain the Sun for any
orientation. The user is asked for the Galactic (or ICRS) coordinates and
distance of the image centre, and the Galactic (or ICRS) coordinates of the
normal to the desired plane. To avoid any ambiguity in the orientation of the
image, the coordinates of the X and Y axes are indicated. Fig. 9 shows an
example of a plane parallel to the mid-plane and displaced by 100 pc to the
south (i.e. crossing the Southern Galactic Pole axis at Z=-100pc.).

##  5 Comparisons with recent results on dust and stars

Figure 10: Dust in the mid-plane (black and white image) and projections onto
the plane of molecular clouds from the Zucker et al. (2020) catalogue (pink
markers). The size of the marker decreases with the distance to the plane
(maximum visible distance abs(Z)=100 pc.) Only a few objects from the
catalogue fall outside dense dust structures: W3, W4, W5 (l=134 to 137 ∘.
d=1-2 kpc), and IC1396 (l=99∘, d=900 pc (but see text). Figure 11: Massive O
and B stars from the recent catalogue of Pantaleoni González et al. (2021),
superimposed on the dust extinction image in the horizontal plane containing
the Sun. Only massive stars with identified stellar types, possessing EDR3
parallaxes, and located within 50 pc of the plane are represented. The colour
scale and the size of the marker distinguish types from O3 (blue, large
circles) to B9 (pink, small circles).

It is important to compare dust cloud mapping based on distinct photometric
data in order to detect potential differences. To do so, we have used the
recent catalogue of molecular cloud (MC) distance assignments of Zucker et al.
(2020). In this catalogue, molecular cloud locations are based on the
PanSTARRS photometric data and Gaia DR2 parallaxes. For clouds within
$\simeq$2 kpc, DR2 and EDR3 parallaxes slightly differ individually; however,
the MC distance assignments are based on large numbers of targets, and errors
are expected to be smoothed out in a large extent. Fig. 10 displays the
projection onto the midplane of the $\geq$300 MCs from the catalogue,
superimposed on the dust extinction in this plane. The figure shows a
remarkable agreement, that is to say clouds are located within or are
extremely close to the main dust concentrations deduced by inversion. There
are only very few exceptions, namely W3 (IC 1795), W4(IC 805), W5 (IC 1848),
and IC1396 (see Fig. 10). We note, however, that those particular clouds are
much more diffuse than other objects of the catalogue, and that assigned
locations of subsets of the first three are distributed over large distances.
We also note some discrepancies between the assigned locations and distances
of the dust fronts found using the Green et al. (2019) online tool in the
quoted directions. We attribute this variability to the faintness and large
extent of these nebulae. Apart from these rare objects, the excellent
agreement demonstrates that using different photometric bands does not
introduce systematic trends in the distance assignments of the obscuring dust
clouds.

As said in the introduction, 3D maps are a general tool and comparisons with
many types of observed quantities have a potential interest. We have selected
one example, the comparison between the distribution of the clouds and
cavities and the locations of massive early-type stars assembled in the recent
Alma catalogue of O and B stars from Pantaleoni González et al. (2021). The
catalogue is ideally suited for comparisons with our maps since the majority
of the targets are located within 3 kpc. Before comparisons with the present
dust distribution, and for consistency with the maps, we have searched for the
Gaia EDR3 parallaxes of all objects from the catalogue and removed the very
small number of objects without EDR3 parallaxes from the list. Fig. 11 shows,
superimposed on the dust map in the mid-plane, the location of O, B stars
identified as massive and located within 50 pc from the plane. A colour scale
distinguishes the spectral types. It can be seen that, in general, O stars are
located close to dense cloud complexes, which is in agreement with the
classical scheme of star formation, according to which young massive star
winds start to carve cavities in their parent cloud complexes immediately
after their birth. There is, however, a difference in the second quadrant,
characterised by a relatively large number of objects between 1 and 4 kpc, and
the absence of dense dust structures. This behavior was already perceptible in
our previous map; however, here the spatial coverage is improved and the
mapping limitations can no longer be at the origin of this peculiarity. Such a
difference deserves further studies, in particular it requires for potential
effects of extinction law variability to be disentangled from actual different
star formation episodes.

Figure 12: Distance along the Z axis of the majority of the dust mass (see
text for a definition). Distances are in parsecs and positive towards the
northern pole. In the third quadrant, most of the dust is below the mid-plane
(but see text), while it is the opposite in the second and fourth quadrants.
Also shown are the projections onto the mid-plane of the boundaries of the
thick vertical slices from Fig. 13, numbered from 1 to 4. Figure 13: Examples
of images of the dust in four vertical slabs. The orientations are indicated
on the right, by means of the Galactic longitude of the corresponding
extremity of the slab. The projections of the slabs onto the mid-plane are
shown in Fig. 12. The wavy pattern of the dust distribution with respect to
the mid-plane is seen to be a general feature.

Finally, a very peculiar and interesting aspect of this mapping is the
vertical distribution of the dust. Fig. 12 displays, for each pair of X,Y
coordinates, the quantity $Z_{mean}$ which represents the average Z coordinate
of the dust (with the Z axis oriented from the southern to northern poles).
More specifically, it was computed as the average value of the Z abscissa of
each cell, weighted by the square of the extinction density in the cell,

$Z_{mean}=\frac{\int_{-400}^{+400}\rho(Z)^{2}ZdZ}{\int_{-400}^{+400}\rho(Z)^{2}dZ}.$
(2)

The figure shows a complex structure, dominated by a conspicuous negative
trend in the third quadrant (dominant yellow and red colour in the figure)
where most of the dust is 100-150 pc below the mid-plane. Conversely, a large
fraction of positive values are found in the other quadrants, with most of the
dust being 150-200 pc above the mid-plane in a number of regions. In some
specific areas, however, this average value has a limited significance due to
the existence of two dust layers, such as in the anti-centre direction. In
addition to this general trend, there are vertical wavy structures in many
regions, with an amplitude in Z particularly large within 1 kpc from the Sun.
This can be seen in Fig. 13 which shows examples of the dust density in four
vertical slabs whose orientations, thicknesses, and intersections with the
mid-plane are indicated in Fig. 12. The quantity represented in Fig. 13 is,
for each Z value, the integral of the extinction density along horizontal axes
perpendicular to the slab orientation in the mid-plane. We have restricted the
figure to four different slabs; however, the vertical oscillations of the dust
layer appear to be a general feature. One of the vertical slabs (number three)
contains the gas clouds forming the Radcliffe wave discussed by Alves et al.
(2020), and it shows a remarkable sinusoidal pattern. It is beyond the scope
of this work to discuss such structures in detail; however, we argue that they
are likely connected with the ripples in the 3D velocity distribution of
nearby stars recently discovered by Antoja et al. (2018). The spectacular
spiral (or snail-shaped) pattern in the Z-VZ plane found by the authors is
reproduced in state-of-the-art simulations of Milky-Way crossings by dwarf
galaxies, but perturbations by the buckling bar and by arms also succeed in
reproducing similar structures (see, e.g. Hunt et al., 2021; Khoperskov et
al., 2019). The period of the spiral along the Z direction can be derived from
Fig. 1 of Antoja et al. (2018) and is of the order of 300-350 pc. On the other
hand, from Fig. 13, one can derive a mean vertical amplitude of the dust
oscillations around the mid-plane of the same order visually, namely $\simeq$
300-350 pc. We believe that this similarity is not accidental and that it is
the signature of a common origin for both gas and stars perturbations. Future
studies and, in particular, adjustments of self-consistent models for both
stars and gas to stellar and interstellar matter distributions may help to
infer the origin of both types of quasi-periodic structures. Fig. 2 from
Antoja et al. (2018) reveals another organised perturbation of the stellar
velocities, namely ridges in the distribution of azimuthal velocities as a
function of the Galactocentric radius. Pursuing the search for correspondences
between the dust distribution and the local stellar kinematics, we have used
the figure to derive a mean spatial period of roughly 550 pc. However, we
could not find clear signs of such a period along the X axis in the extinction
density map. This is probably due to the complexity of the dust structures and
especially the existence of two very different orientations of the main chains
of dust clouds, a duality we already pointed out based on the DR2 maps: a
l$\simeq$80∘direction similar to the spiral arm orientation and, second, a
l$\simeq$45∘orientation illustrated by the ’split’ and a large fraction of
Sagittarius (see Fig. 10 and 11 from Lallement et al. (2019)). Such a complex
pattern suggests the existence of past perturbations, other than the one
having generated the snail-shaped pattern.

##  6 Conclusions and perspectives

Thanks to the new version of the Gaia data (EDR3), and in particular to the
more accurate parallaxes and more homogeneous photometry, we have updated our
calculations of the 3D extinction density distribution in the local
interstellar medium at $\lambda$=5500Å. The computational volume is a 6 kpc x
6 kpc x 800 pc slab with the Sun at the centre. We have also provided error
estimates on both the extinction densities (or differential extinction, in
$mag\cdot pc^{-1}$) and on the integrated extinctions on the lines of sight
(from the Sun to any location within the volume). The maps are downloadable in
their entirety (3D distribution of 5pcx5pcx5pc voxels) and/or available for
inspection with several tools on the EXPLORE website555https://explore-
platform.eu. In particular, it is possible to download integrated extinction
profiles for series of directions captured in a single file, as well as dust
distribution images in selected planes.

The main improvements to the maps are improved cloud reconstruction quality,
increased distances at which clouds are reconstructed, and significantly
improved local cavity reconstruction. The new map confirms the highly complex
3D distribution of dust at 3-4 kpc, an unexpected structure before Gaia. The
45∘ pitch angle cloud chain nicknamed ’split’ and the huge dust-free ’avenue’
between ’split’ and Sagittarius are clearly visible. The undulations of the
dust sheet around the plane are extensive. The average vertical amplitude of
these oscillations around the mid-plane is of the order of $\Delta$Z=300 pc,
which is similar to the vertical period found in the 3D velocity distribution
of nearby stars by Antoja et al. (2018), namely the spectacular spiral (or
¿snail) pattern of azimuthal (tangential) and radial velocities in the Z-VZ
plane. We suggest a common origin for these two phenomena. If this is the
case, the grain and gas distributions could be associated with the kinetic 3D
tomography of the stars in order to help clarify their origin.

With the low-resolution blue and red spectrographs on board the satellite,
future Gaia data releases will provide a much larger number of reliable
extinction measurements. The dual hierarchical method has proven to be well
suited to the inversion of large input data sets due to its low computational
time and will be well adapted to these data. In the future, one may expect
much improved dust distributions, ideally accompanied by co-spatial gas
kinematics. Such results may enter self-consistent models of stars and gas to
help decipher the origin of the present-day Galactic structure and the epochs
and impacts of disk traversals and bar and arm disruptions. In this regard,
the construction of fully self-consistent distributions of extinction and
stars (Babusiaux et al., 2020; Hottier et al., 2020) will be a crucial
improvement.

###### Acknowledgements.

We thank our anonymous referee for useful comments.
J.L.V. and N.L.J.C. acknowledge support from the EXPLORE project. EXPLORE has
received funding from the European Union’s Horizon 2020 research and
innovation programme under grant agreement No 101004214. J.L.V. acknowledges
support from the THETA Observatoire des Sciences de l’Univers in Besançon.
This work has made use of data from the European Space Agency (ESA) mission
Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing
and Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement. This work also makes use of
data products from the 2MASS, which is a joint project of the University of
Massachusetts and the Infrared Processing and Analysis Center/California
Institute of Technology, funded by the National Aeronautics and Space
Administration and the National Science Foundation.
This research has made use of the SIMBAD database, operated at CDS,
Strasbourg, France.

## References

  * Alves et al. (2020) Alves, J., Zucker, C., Goodman, A. A., et al. 2020, Nature, 578, 237
  * Antoja et al. (2018) Antoja, T., Helmi, A., Romero-Gómez, M., et al. 2018, Nature, 561, 360
  * Babusiaux et al. (2020) Babusiaux, C., Fourtune-Ravard, C., Hottier, C., Arenou, F., & Gómez, A. 2020, A&A, 641, A78
  * Capitanio et al. (2017) Capitanio, L., Lallement, R., Vergely, J. L., Elyajouri, M., & Monreal-Ibero, A. 2017, A&A, 606, A65
  * Chen et al. (2019) Chen, B. Q., Huang, Y., Yuan, H. B., et al. 2019, MNRAS, 483, 4277
  * Chen et al. (2020) Chen, B. Q., Li, G. X., Yuan, H. B., et al. 2020, MNRAS, 493, 351
  * Fitzpatrick et al. (2019) Fitzpatrick, E. L., Massa, D., Gordon, K. D., Bohlin, R., & Clayton, G. C. 2019, ApJ, 886, 108
  * Gaia Collaboration et al. (2021) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2021, A&A, 649, A1
  * Green et al. (2019) Green, G. M., Schlafly, E., Zucker, C., Speagle, J. S., & Finkbeiner, D. 2019, ApJ, 887, 93
  * Guo et al. (2021) Guo, H. L., Chen, B. Q., Yuan, H. B., et al. 2021, ApJ, 906, 47
  * Hottier et al. (2020) Hottier, C., Babusiaux, C., & Arenou, F. 2020, A&A, 641, A79
  * Hunt et al. (2021) Hunt, J. A. S., Stelea, I. A., Johnston, K. V., et al. 2021, MNRAS, 508, 1459
  * Khoperskov et al. (2019) Khoperskov, S., Di Matteo, P., Gerhard, O., et al. 2019, A&A, 622, L6
  * Lallement et al. (2019) Lallement, R., Babusiaux, C., Vergely, J. L., et al. 2019, A&A, 625, A135
  * Leike et al. (2020) Leike, R. H., Glatzle, M., & Enßlin, T. A. 2020, A&A, 639, A138
  * Lindegren et al. (2021) Lindegren, L., Bastian, U., Biermann, M., et al. 2021, A&A, 649, A4
  * Pantaleoni González et al. (2021) Pantaleoni González, M., Maíz Apellániz, J., Barbá, R. H., & Reed, B. C. 2021, MNRAS, 504, 2968
  * Rezaei Kh. et al. (2020) Rezaei Kh., S., Bailer-Jones, C. A. L., Soler, J. D., & Zari, E. 2020, A&A, 643, A151
  * Riello et al. (2021) Riello, M., De Angeli, F., Evans, D. W., et al. 2021, A&A, 649, A3
  * Ruiz-Dern et al. (2018) Ruiz-Dern, L., Babusiaux, C., Arenou, F., Turon, C., & Lallement, R. 2018, A&A, 609, A116
  * Sale & Magorrian (2018) Sale, S. E. & Magorrian, J. 2018, MNRAS, 481, 494
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163
  * Vergely et al. (2001) Vergely, J.-L., Freire Ferrero, R., Siebert, A., & Valette, B. 2001, A&A, 366, 1016
  * Zucker et al. (2020) Zucker, C., Speagle, J. S., Schlafly, E. F., et al. 2020, A&A, 633, A51