#  Where is the lightest charmed scalar meson?

Meng-Lin Du<EMAIL_ADDRESS>Helmholtz-Institut für Strahlen- und
Kernphysik and Bethe Center for Theoretical Physics,
Universität Bonn, D-53115 Bonn, Germany  Feng-Kun Guo<EMAIL_ADDRESS>CAS Key
Laboratory of Theoretical Physics, Institute of Theoretical Physics,
Chinese Academy of Sciences, Beijing 100190, China  School of Physical
Sciences, University of Chinese Academy of Sciences, Beijing 100049, China
Christoph Hanhart<EMAIL_ADDRESS>Institute for Advanced Simulation,
Institut für Kernphysik and Jülich Center for Hadron Physics,
Forschungszentrum Jülich, D-52425 Jülich, Germany  Bastian Kubis
<EMAIL_ADDRESS>Helmholtz-Institut für Strahlen- und Kernphysik and
Bethe Center for Theoretical Physics,
Universität Bonn, D-53115 Bonn, Germany  Ulf-G. Meißner<EMAIL_ADDRESS>bonn.de Helmholtz-Institut für Strahlen- und Kernphysik and Bethe Center for
Theoretical Physics,
Universität Bonn, D-53115 Bonn, Germany  Institute for Advanced Simulation,
Institut für Kernphysik and Jülich Center for Hadron Physics,
Forschungszentrum Jülich, D-52425 Jülich, Germany  Tbilisi State University,
0186 Tbilisi, Georgia

###### Abstract

The lightest charmed scalar meson is known as the $D_{0}^{*}(2300)$, which is
one of the earliest new hadron resonances observed at modern $B$ factories. We
show here that the parameters assigned to the lightest scalar $D$-meson are in
conflict with the precise LHCb data of the decay $B^{-}\to
D^{+}\pi^{-}\pi^{-}$. On the contrary, these data can be well described by an
unitarized chiral amplitude containing a much lighter charmed scalar meson,
the $D_{0}^{*}(2100)$. We also extract the low-energy $S$-wave $D\pi$ phase of
the decay $B^{-}\to D^{+}\pi^{-}\pi^{-}$ from the data in a model-independent
way, and show that its difference from the $D\pi$ scattering phase shift can
be traced back to an intermediate $\rho^{-}$ exchange. Our work highlights
that an analysis of data consistent with chiral symmetry, unitarity, and
analyticity is mandatory in order to extract the properties of the ground-
state scalar mesons in the singly heavy sector correctly, in analogy to the
light scalar mesons $f_{0}(500)$ and $K_{0}^{*}(700)$.

Introduction.—Since the discovery of the $D_{s0}^{*}(2317)$ Aubert et al.
(2003), many hadrons were observed beyond the quark model expectations, which
have seriously challenged the understanding of the hadron spectrum in terms of
the conventional quark model that identifies mesons as $\bar{q}q$ states. The
observation that the $D_{s0}^{*}(2317)$ Aubert et al. (2003) and
$D_{s1}(2460)$ Besson et al. (2003) are significantly lighter than expected by
the quark model, around $2.48$ and $2.55\,\mbox{GeV}$ Godfrey and Isgur
(1985); Godfrey and Moats (2016); Ebert et al. (2010), has driven the
development of various models, including $D^{(\ast)}K$ hadronic molecules
Barnes et al. (2003); van Beveren and Rupp (2003); Szczepaniak (2003);
Kolomeitsev and Lutz (2004); Chen and Li (2004); Guo et al. (2006, 2007);
Gamermann et al. (2007), tetraquark states Cheng and Hou (2003); Maiani et al.
(2005), and mixtures of $c\bar{q}$ with tetraquarks Browder et al. (2004). In
2004, two new charm-nonstrange structures, the $D_{0}^{*}(2300)$ Abe et al.
(2004); Link et al. (2004), called $D_{0}^{*}(2400)$ previously, and
$D_{1}(2430)$ Abe et al. (2004), were reported as the SU(3) partners of the
$D_{s0}^{*}(2317)$ and $D_{s1}(2460)$, respectively. The observations posed a
puzzle: why are the masses of the two nonstrange mesons, $D_{0}^{*}(2300)$ and
$D_{1}(2430)$, almost equal to their strange siblings, i.e., the
$D_{s0}^{*}(2317)$ and $D_{s1}(2460)$? Thanks to new data from both lattice
quantum chromodynamics (QCD) Liu et al. (2013); Mohler et al. (2013); Lang et
al. (2014); Moir et al. (2016); Bali et al. (2017); Cheung et al. (2020) and
the LHCb experiment Aaij et al. (2016), it was recently demonstrated that the
various puzzles in the charm meson spectrum can be solved naturally in the
framework of unitarized chiral perturbation theory (UChPT) Albaladejo et al.
(2017); Du et al. (2018); Guo et al. (2018a) that allows one to calculate the
nonperturbative dynamics of Goldstone bosons scattering off the
$D_{(s)}^{(*)}$ mesons in a controlled way. The combination of UChPT and
lattice QCD not only reproduced the correct $D_{s0}^{*}(2317)$ mass Liu et al.
(2013), but also predicted its pion mass dependence Du et al. (2017); Bali et
al. (2017). The solution provided for the SU(3) mass hierarchy puzzle
mentioned above is that instead of only one heavy state, $D_{0}^{*}(2300)$ in
the channel $(S,I)\equiv(\text{strangeness},\text{isospin})=(0,1/2)$, there
are two states, one lighter and one heavier Kolomeitsev and Lutz (2004); Guo
et al. (2006); Gamermann et al. (2007); Guo et al. (2009); Liu et al. (2013);
Guo et al. (2015); Albaladejo et al. (2017); Guo et al. (2019, 2018b, 2018c).
The most recent studies revealed their pole locations to be at
$\left(2105^{+6}_{-8}-i\,102^{+10}_{-11}\right)$ and
$\left(2451^{+35}_{-26}-i\,134^{+7}_{-8}\right)\text{MeV}$ Albaladejo et al.
(2017); Du et al. (2018), respectively. The SU(3) partner of the
$D_{s0}^{*}(2317)$ is the lighter one, denoted as $D_{0}^{*}(2100)$ in the
following, which restores the expected mass hierarchy. The heavier pole on the
other hand is a member of a different multiplet. Support for the presence of
two poles comes from an analysis of the high-quality LHCb data on the decays
$B^{-}\to D^{+}\pi^{-}\pi^{-}$ Aaij et al. (2016),
$B_{s}^{0}\to\bar{D}^{0}K^{-}\pi^{+}$ Aaij et al. (2014),
$B^{0}\to\bar{D}^{0}\pi^{-}\pi^{+}$ Aaij et al. (2015a), $B^{-}\to
D^{+}\pi^{-}K^{-}$ Aaij et al. (2015b), and $B^{0}\to\bar{D}^{0}\pi^{-}K^{+}$
Aaij et al. (2015c) performed in Refs. Du et al. (2018, 2019), as well as from
the fact that their existence is consistent with the lattice energy levels Liu
et al. (2013); Mohler et al. (2013); Lang et al. (2014); Moir et al. (2016);
Bali et al. (2017) for the relevant two-body scattering Albaladejo et al.
(2017, 2018); Guo et al. (2019). This two-pole structure indeed emerges as a
more general pattern in the hadron spectrum, see, e.g., Ref. Meißner (2020).

Despite the phenomenological success of this picture in describing the
available lattice and LHCb data, the observation that the lightest $D_{0}^{*}$
has a mass around $2.1\,\mbox{GeV}$ has not entered the Review of Particle
Physics (RPP) Zyla et al. (2020) yet, which still lists the $D_{0}^{*}(2300)$
as the lightest charmed scalar meson and the $D_{1}(2430)$ as the
corresponding axial-vector meson. In this Letter, we demonstrate that the
$D_{0}^{*}(2300)$ as in the RPP is not consistent with the most precise data
for $B^{-}\to D^{+}\pi^{-}\pi^{-}$, contrary to the $D_{0}^{*}(2100)$
predicted in UChPT, and conclude that the positive-parity charm-nonstrange
meson spectrum in the RPP needs to be revised.

$D\pi$ $S$-wave phase of $B^{-}\to D^{+}\pi^{-}\pi^{-}$.—The decay amplitude
in the low-energy region of the $D\pi$ system can be decomposed into $S$-,
$P$-, and $D$-waves,

$\displaystyle\mathcal{A}_{B^{-}\to
D^{+}\pi^{-}\pi^{-}}(s,z)=\sum_{\ell=0}^{2}\sqrt{2\ell+1}\mathcal{A}_{\ell}(s)P_{\ell}(z_{s})\,,$
(1)

where $\mathcal{A}_{\ell}(s)$ with $\ell=0,1,2$ correspond to the amplitudes
with $D^{+}\pi^{-}$ in the $S$-, $P$-, and $D$-waves, respectively, $s$ is the
c.m. energy squared of the $D^{+}\pi^{-}$ system, and $P_{\ell}(z_{s})$ are
the Legendre polynomials with $z_{s}$ the cosine of the helicity angle of the
$D^{+}\pi^{-}$ system, i.e., the angle between the moving directions of the
two pions in the $D^{+}\pi^{-}$ c.m. frame. The angular moments are determined
by weighting the data with the Legendre polynomials $P_{\ell}(z)$ Aaij et al.
(2016). They contain contributions from certain partial waves and their
interference terms, and thus the corresponding phase variations. The first few
moments are given by Aaij et al. (2016); Du et al. (2018, 2019)

$\displaystyle\langle P_{0}\rangle$
$\displaystyle\propto|\mathcal{A}_{0}|^{2}+|\mathcal{A}_{1}|^{2}+|\mathcal{A}_{2}|^{2}\,,$
$\displaystyle\langle P_{2}\rangle$
$\displaystyle\propto\frac{2}{5}|\mathcal{A}_{1}|^{2}+\frac{2}{7}|\mathcal{A}_{2}|^{2}+\frac{2}{\sqrt{5}}|\mathcal{A}_{0}||\mathcal{A}_{2}|\cos(\delta_{2}-\delta_{0})\,,$
$\displaystyle\langle P_{13}\rangle$ $\displaystyle\equiv\langle
P_{1}\rangle-\frac{14}{9}\langle
P_{3}\rangle\propto\frac{2}{\sqrt{3}}|\mathcal{A}_{0}||\mathcal{A}_{1}|\cos(\delta_{1}-\delta_{0})\,,$
(2)

with $\delta_{i}$ the phase of $\mathcal{A}_{i}$, i.e.,
$\mathcal{A}_{i}=|\mathcal{A}_{i}|e^{i\delta_{i}}$. As first proposed in Ref.
Du et al. (2018), we use the linear combination $\langle P_{13}\rangle$
instead of $\langle P_{1}\rangle$ and $\langle P_{3}\rangle$ individually,
since it only depends on the $S$-$P$-wave interference up to $\ell=2$ and is
particularly sensitive to the $S$-wave phase motion.

Figure 1:  Comparison of the predictions of Eq. (5) from UChPT (blue) and a
Breit-Wigner parametrization (green) for $\delta_{0}$ with the phase extracted
in Ref. Aaij et al. (2016) (red) and that using Eq. (3) (black). The bands
correspond to errors propagated from the input UChPT scattering amplitudes and
from the Breit-Wigner resonance parameters.

For $M_{D^{+}\pi^{-}}<2.2\,\mbox{GeV}$, even the $D$-wave can be neglected,
since the narrow tensor resonance $D_{s2}(2460)$ is sufficiently far away.
This can be verified from the data of the angular moments, i.e., Fig. 3 in
Ref. Aaij et al. (2016). Therefore, in this kinematic regime one obtains

$\displaystyle\cos(\delta_{0}-\delta_{1})=\sqrt{\frac{3}{10}}\frac{\langle
P_{13}\rangle}{\sqrt{\langle P_{2}\rangle}\sqrt{\langle
P_{0}\rangle-\frac{5}{2}\langle P_{2}\rangle}}$ (3)

for the $S$-$P$ phase difference. The $P$-wave is dominated by the vector
resonance $D^{*}(2007)^{0}$ below the $D^{+}\pi^{-}$ threshold with a width of
less than $60\,\mbox{keV}$ Rosner (2013); Guo (2019). The next vector $D^{*}$
resonance is far above this energy region. Thus, the phase of the $P$-wave
$\delta_{1}$ can be safely fixed to $180^{\circ}$ for the region we are
interested in. The $S$-wave $D\pi$ phase motion below 2.2 GeV can then be
extracted, see Fig. 1. For comparison, the phase motion of the $D^{+}\pi^{-}$
$S$-wave up to 2.4 GeV obtained in the LHCb analysis Aaij et al. (2016) (with
the phase at 2.4 GeV fixed to $180^{\circ}$) is also shown, which is fully in
line with the phase we extracted from Eq. (3) below 2.2 GeV.

Figure 2: The decay $B^{-}\to D^{+}\pi^{-}\pi^{-}$ via the coupled channel
$B^{-}\to D^{0}\pi^{0}\pi^{-}$. The filled square denotes the $D^{0}\pi^{0}\to
D^{+}\pi^{-}$ $T$-matrix element.

For $M_{D^{+}\pi^{-}}<2.4\,\mbox{GeV}$, the effect of the $\rho$ meson could
be significant via the coupled channel $B^{-}\to D^{0}\pi^{0}\pi^{-}$; see
Fig. 2. This follows directly from the large branching ratio
$\mathcal{B}(B^{-}\to D^{0}\rho^{-})=1.34\%$, which is an order larger than
$\mathcal{B}(B^{-}\to D^{+}\pi^{-}\pi^{-})=0.107\%$ Zyla et al. (2020). It is
therefore reasonable to assume that the decay $B^{-}\to D^{+}\pi^{-}\pi^{-}$
is dominated by the process $B^{-}\to D^{0}\rho^{-}\to D^{0}\pi^{0}\pi^{-}\to
D^{+}\pi^{-}\pi^{-}$. By virtue of soft-pion theorems, one has Lin and Chiang
(1983)

$\displaystyle\mathcal{A}(B^{-}\to D^{+}\pi^{-}\pi^{-})\big{|}_{p_{\pi^{-}}\to
0}$
$\displaystyle=\frac{1}{F_{\pi}}\mathcal{A}(B^{0}\to\bar{D}^{0}\pi^{0})\,,$
$\displaystyle\mathcal{A}(B^{-}\to D^{0}\pi^{0}\pi^{-})\big{|}_{p_{\pi^{0}}\to
0}$ $\displaystyle=-\frac{1}{F_{\pi}}\mathcal{A}(B^{-}\to D^{0}\pi^{-})\,,$
(4)

where $p_{\pi^{-}(\pi^{0})}$ is the momentum of the $\pi^{-}(\pi^{0})$, and
$F_{\pi}$ is the pion decay constant (in the chiral limit). From
$\mathcal{B}(B^{0}\to\bar{D}^{0}\pi^{0})=2.63\times 10^{-4}$ and
$\mathcal{B}(B^{-}\to D^{0}\pi^{-})=4.68\times 10^{-3}$, one concludes that at
low energies for the $D^{+}\pi^{-}(D^{0}\pi^{0})$ system, the amplitude of
$B^{-}\to D^{0}\pi^{0}\pi^{-}$ is much larger than that of $B^{-}\to
D^{+}\pi^{-}\pi^{-}$. Furthermore, isospin symmetry shows that for the decays
$B\to D\pi\pi$ with even relative angular momenta between the pions, the
amplitude for $B^{-}\to D^{+}\pi^{-}\pi^{-}$ is larger than that of $B^{-}\to
D^{0}\pi^{0}\pi^{-}$ by a factor of $2\sqrt{2}$ Savage and Wise (1989);
Niecknig and Kubis (2015); Du et al. (2018, 2019). As in addition even angular
momenta here imply isospin $I=2$ and therefore nonresonant partial waves, the
relative angular momentum of $\pi^{0}\pi^{-}$ in the decay $B^{-}\to
D^{0}\pi^{0}\pi^{-}$ is by far dominantly odd in the low-energy regime for
$D^{0}\pi^{0}$, and the $\rho^{-}$ plays a crucial role.

If we assume that the decay $B^{-}\to D^{+}\pi^{-}\pi^{-}$ is dominated by the
process in Fig. 2, the $D\pi$ $S$-wave part of the triangle diagram can be
estimated by the integral

$\displaystyle\mathcal{A}_{0}^{\text{trig}}(s)=\frac{1}{\pi}\int_{s_{\text{th}}}^{\infty}ds^{\prime}\frac{\hat{P}(s^{\prime})\rho(s^{\prime})T_{D^{0}\pi^{0}\to
D^{+}\pi^{-}}(s^{\prime})}{s^{\prime}-s},$ (5)

where $\hat{P}(s)$ is the production amplitude for $B^{-}\to D^{0}\rho^{-}\to
D^{0}\pi^{0}\pi^{-}$ projected to the $D^{0}\pi^{0}$ $s$-channel,
$\rho(s)=\sqrt{\lambda(s,M_{D}^{2},M_{\pi}^{2})}/{(16\pi s)}$ is the $D\pi$
phase space with $\lambda(a,b,c)=a^{2}+b^{2}+c^{2}-2ab-2ac-2bc$ the Källén
function, $T_{D^{0}\pi^{0}\to D^{+}\pi^{-}}(s)$ the $S$-wave scattering
amplitude for $D^{0}\pi^{0}\to D^{+}\pi^{-}$, and
$s_{\text{th}}=(M_{D}+M_{\pi})^{2}$. The expression for $\hat{P}(s)$ is the
same as $\hat{\mathcal{F}}_{0}^{1/2}(s)$ in Eq. (12) below.

The evaluation of Eq. (5) depends on the asymptotic behavior of the integrand,
which is divergent in general. We may estimate Eq. (5) using a cutoff at
$\sqrt{s_{\text{max}}}=\sqrt{q_{\text{max}}^{2}+M_{D}^{2}}+\sqrt{q_{\text{max}}^{2}+M_{\pi}^{2}}$,
where $q_{\text{max}}\approx 1\,\mbox{GeV}$ (another way is to introduce a
form factor, e.g., $e^{-(s-s_{\text{th}})/s_{0}}$ with
$s_{0}=\mathcal{O}(1\,\mbox{GeV})$ Szczepaniak (2015)). We evaluate Eq. (5) by
employing both the $D\pi$ scattering amplitude from UChPT Liu et al. (2013)
and that of a Breit-Wigner (BW) parametrization of the $D_{0}^{*}(2300)$ for
comparison, despite the deficiencies of the latter discussed in Ref. Du et al.
(2019); see also Ref. Gardner and Meißner (2002).

The results with $q_{\text{max}}=1\,\mbox{GeV}$ are shown in Fig. 1, where the
solid blue band and the green dashed band correspond to the $D\pi$ scattering
amplitudes from UChPT and BW, respectively. The obtained phase describes the
data perfectly for the UChPT amplitude, while the BW one fails. We have
checked that the obtained phases are insensitive to a variation of the cutoff
in a reasonable region, $q_{\text{max}}\in[0.8,1.2]\,\mbox{GeV}$.

Khuri-Treiman formalism.—While Eq. (5) provides a reasonable estimation of the
$S$-wave decay amplitude with a clear underlying physical picture, it does not
respect three-body unitarity. In order to check if the conclusion formulated
above is robust, we cure this deficiency by employing the Khuri-Treiman
equations Khuri and Treiman (1960), which are based on two-body elastic phase
shifts and explicitly generate the crossed-channel rescattering between final-
state particles. The formulas are constructed from dispersion relations for
the related crossed scattering processes and then analytically continued to
the decay region, referring to the continuation of the triangle graph Bronzan
and Kacser (1963).

We can write amplitudes for $\mathcal{A}_{+--}(B^{-}\to D^{+}\pi^{-}\pi^{-})$
and $\mathcal{A}_{00-}(B^{-}\to D^{0}\pi^{0}\pi^{-})$ in terms of single-
variable functions according to a reconstruction theorem Niecknig and Kubis
(2015, 2018),

$\displaystyle\mathcal{A}_{+--}$ $\displaystyle(s,$ $\displaystyle
t,u)=\mathcal{F}_{0}^{1/2}(s)+\frac{\kappa(s)}{4}z_{s}\mathcal{F}_{1}^{1/2}(s)$
$\displaystyle+$
$\displaystyle\frac{\kappa(s)^{2}}{16}(3z_{s}^{2}-1)\mathcal{F}_{2}^{1/2}(s)+(t\leftrightarrow
s)\,,$ $\displaystyle\mathcal{A}_{00-}$ $\displaystyle(s,$ $\displaystyle
t,u)=-\frac{1}{\sqrt{2}}\mathcal{F}_{0}^{1/2}(s)-\frac{\kappa(s)}{4\sqrt{2}}z_{s}\mathcal{F}_{1}^{1/2}(s)$
$\displaystyle-$
$\displaystyle\frac{\kappa(s)^{2}}{16\sqrt{2}}(3z_{s}^{2}-1)\mathcal{F}_{2}^{1/2}(s)+\frac{\kappa_{u}(u)}{4}z_{u}\mathcal{F}_{1}^{1}(u)\,,$

where the subindex $\ell$ and superindex $I$ of the single-variable amplitudes
$\mathcal{F}_{\ell}^{I}$ represent the angular momentum and isospin,
respectively, and only the $I<3/2$ and $\ell\leq 2$ terms are taken into
account. The Mandelstam variables of the $B$-meson decay $B^{-}(p_{B})\to
D(p_{D})\pi(p_{1})\pi^{-}(p_{2})$ are $s=(p_{B}-p_{2})^{2}$,
$t=(p_{B}-p_{1})^{2}$, and $u=(p_{B}-p_{D})^{2}$. The corresponding angles are
given by

$z_{s}\equiv\cos\theta_{s}=\frac{s(t-u)-\Delta}{\kappa(s)}\,,\quad
z_{u}\equiv\cos\theta_{u}=\frac{t-s}{\kappa_{u}(u)}\,,$ (7)

where
$\kappa(s)=\lambda^{1/2}(s,M_{D}^{2},M_{\pi}^{2})\lambda^{1/2}(s,M_{B}^{2},M_{\pi}^{2})$,
$\kappa_{u}(u)=\lambda^{1/2}(u,M_{B}^{2},M_{D}^{2})\sqrt{1-{4M_{\pi}^{2}}/{u}}$,
and $\Delta=(M_{B}^{2}-M_{\pi}^{2})(M_{D}^{2}-M_{\pi}^{2})$.

Since we are interested in the $s$-channel process, we use the index $A$ ($B$)
to label the two-body channels corresponding to $D^{+}\pi^{-}$ and
$D^{0}\pi^{0}$. The partial-wave decomposition for the decay amplitudes
$\mathcal{A}_{A}$ reads

$\displaystyle\mathcal{A}_{A}(s,z_{s})=\sum_{I,\ell}b_{I,\ell}^{A}P_{\ell}(z_{s})f_{\ell}^{I}(s)\,,$
(8)

with $b_{I,\ell}^{A}$ denoting Clebsch-Gordan coefficients. By comparing with
Eq. (1), it is easy to obtain
$\mathcal{A}_{\ell}(s)=(2\ell+1)^{-1/2}\sum_{I}b_{I,\ell}^{1}f_{\ell}^{I}(s)$.
We have the following partial-wave unitarity relation for elastic
rescattering:

$\displaystyle\text{disc}f_{\ell}^{I}(s)=2if_{\ell}^{I}(s)\sin\delta_{\ell}^{I}(s)e^{-i\delta_{\ell}^{I}(s)}\theta(s-s_{\text{th}})\,,$
(9)

where $\delta_{\ell}^{I}(s)$ is the elastic final-state scattering phase
shift. The discontinuities of $f_{\ell}^{I}$ and those of the single-variable
functions $\kappa^{\ell}\mathcal{F}_{\ell}^{I}$ coincide on the right-hand cut
by construction. Thus, one has

$\text{disc}\mathcal{F}_{\ell}^{I}(s)=2i\big{[}\mathcal{F}_{\ell}^{I}(s)+\hat{\mathcal{F}}_{\ell}^{I}(s)\big{]}\sin\delta_{\ell}^{I}(s)e^{-i\delta_{\ell}^{I}(s)}\theta(s-s_{\text{th}})\,,$
(10)

where the inhomogeneities $\hat{\mathcal{F}}_{\ell}^{I}(s)$ encode the left-
hand cut contributions and are free of discontinuities on the right-hand cut.
This discontinuity relation is solved by

$\mathcal{F}_{\ell}^{I}(s)=\Omega_{\ell}^{I}(s)\bigg{\\{}Q_{\ell}^{I}(s)+\frac{s^{n}}{\pi}\int_{s_{\text{th}}}^{\infty}\frac{ds^{\prime}}{s^{\prime
n}}\frac{\sin\delta_{\ell}^{I}(s^{\prime})\hat{\mathcal{F}}_{\ell}^{I}(s^{\prime})}{|\Omega_{\ell}^{I}(s^{\prime})|(s^{\prime}-s)}\bigg{\\}},$
(11)

where
$\Omega_{\ell}^{I}(s)=\exp\big{\\{}{s}/{\pi}\int_{s_{\text{th}}}^{\infty}ds^{\prime}{\delta_{\ell}^{I}(s^{\prime})}/[s^{\prime}(s^{\prime}-s)]\big{\\}}$
is the Omnès function Omnès (1958), $Q_{\ell}^{I}(s)$ is a polynomial at least
of degree $(n-1)$ (see discussion below), and the number of subtractions $n$
is chosen to guarantee the convergence of the dispersion integral.

The inhomogeneity $\hat{\mathcal{F}}_{\ell}^{I}$ is determined by the partial-
wave decomposition of Eq. (Where is the lightest charmed scalar meson?) as the
projection of the crossed-channel amplitudes onto the considered channel.
Around the $D\pi$ threshold in the $s$-channel, $\sqrt{t}\sim 5\,\mbox{GeV}$,
there is no resonance in the $t$-channel, and thus the interaction is supposed
to be very weak. The only possible significant crossed-channel effect is from
the $\rho$ meson through $B^{-}\to D^{0}\pi^{0}\pi^{-}$. The resulting
inhomogeneity for the $S$-wave $s$-channel amplitude is Niecknig and Kubis
(2015)

$\displaystyle\hat{\mathcal{F}}_{0}^{1/2}(s)=-\frac{1}{4\sqrt{2}}\int_{-1}^{1}dz_{s}(t-s)\mathcal{F}_{1}^{1}(u)\,.$
(12)

For technical details regarding this integral, see Ref. Kacser (1963) and the
Supplemental Material sup .

The full solution for the decay amplitudes can be obtained by solving a set of
coupled integral equations in terms of a few linearly independent complex
subtraction constants contained in $Q_{\ell}^{I}(s)$, which cannot be
determined a priori in the framework of dispersion theory. Since we are only
interested in the $D\pi$ low-energy regime and especially in its $S$-wave,
based on the large branching ratio of $B^{-}\to D^{0}\rho^{-}$, it is
reasonable to approximate $\mathcal{F}_{1}^{1}(u)$ in Eq. (Where is the
lightest charmed scalar meson?) by a BW function for the $\rho$ meson. In this
case, $\mathcal{F}_{1}^{1}(u)$ behaves as $u^{-1}$ for $u\to\infty$, thus
$\hat{\mathcal{F}}_{\ell}^{I}(s)$ in Eq. (12) approaches a constant as
$s\to\infty$. The number of the subtractions $n$ in Eq. (11) is then
determined by the asymptotic behavior of the scattering phase
$\delta_{\ell}^{I}(s)$. For the BW phase and that of UChPT taken from Ref. Liu
et al. (2013), one single subtraction is sufficient. Moreover, the phase of
UChPT, as well as that of the BW, is unreliable at high energies. Thus, the
dispersion integral will be evaluated up to a cutoff $\Lambda$, and the effect
of cutting off the integral may be absorbed into the polynomial
$Q_{\ell}^{I}(s)$. Explicitly, for an integral

$\displaystyle g(s)$
$\displaystyle=\int_{s_{\text{th}}}^{\infty}ds^{\prime}\frac{f(s^{\prime})}{s^{\prime}-s}=\int_{s_{\text{th}}}^{\Lambda}ds^{\prime}\frac{f(s^{\prime})}{s^{\prime}-s}+\int_{\Lambda}^{\infty}ds^{\prime}\frac{f(s^{\prime})}{s^{\prime}-s}$
$\displaystyle\approx
g_{0}+g_{1}s+\int_{s_{\text{th}}}^{\Lambda}ds^{\prime}\frac{f(s^{\prime})}{s^{\prime}-s}\,.$
(13)

For simplicity, we neglect the $I=3/2$ contribution since it contains no
resonances. Therefore, for the $S$-wave amplitude at low energies, Eq. (11)
can be written as

$\displaystyle\mathcal{F}^{1/2}_{0}(s)$
$\displaystyle=\Omega_{0}^{1/2}(s)\bigg{\\{}g_{0}+g_{1}\frac{s-M_{D}^{2}}{M_{D}^{2}}$
$\displaystyle\quad+\frac{s}{\pi}\int_{s_{\text{th}}}^{\Lambda}\frac{ds^{\prime}}{s^{\prime}}\frac{\sin\delta_{0}^{1/2}(s^{\prime})\hat{\mathcal{F}}_{0}^{1/2}(s^{\prime})}{|\Omega_{0}^{1/2}(s^{\prime})|(s^{\prime}-s)}\bigg{\\}}\,.$
(14)

The constants $g_{0}$ and $g_{1}$ have to be fixed by data.

We fit $\langle P_{0}\rangle$, $\langle P_{13}\rangle$, and $\langle
P_{2}\rangle$ of the decay $B^{-}\to D^{+}\pi^{-}\pi^{-}$ up to 2.4 GeV, which
is below the $D\eta$ and $D_{s}\bar{K}$ thresholds. To describe the angular
moments, one needs the explicit amplitudes for the $D\pi$ $P$\- and $D$-waves.
The $P$-wave can be safely parametrized as
$\delta_{1}^{1/2}(s)=\pi\,\theta\big{(}s-M_{D^{*0}}^{2}\big{)}$ as discussed
below Eq. (3). Consequently, the dispersion integral (11) for the $P$-wave can
be neglected since $\sin\pi=0$. The $D$-wave is dominated by the resonance
$D_{2}^{0}(2460)$ with a width of $47.5\,\mbox{MeV}$ Zyla et al. (2020), which
is above the region we are interested in. Thus, the $D$-wave phase is close to
$0$ below 2.4 GeV, and the corresponding dispersion integral can be neglected
as well. Therefore, for the $P$-($D$-)wave amplitudes
$\mathcal{F}_{1,2}^{1/2}$, we can use the same BW forms as those in the LHCb
analysis Aaij et al. (2016), which is equivalent to the corresponding Omnès
function multiplied by a polynomial. In the isobar model used in Ref. Aaij et
al. (2016), complex factors are introduced for each resonance BW function.
Without crossed-channel effects, these factors become real according to
Watson’s theorem Watson (1954). For the $P$\- and $D$-waves, as discussed
above, the Omnès representation should be a good approximation in the energy
region we are interested in, and the normalization factor is real. We also
consider a complex normalization and find the results unchanged.

Figure 3: Results of UChPT (blue) and the BW description (green), with the
best fits $\chi^{2}/\text{d.o.f.}=1.2$ and 2.0, respectively.

For the $D\pi$ $S$-wave, we employ both the scattering phase shifts from UChPT
Liu et al. (2013), which contains the $D_{0}^{*}(2100)$, and the BW for the
$D_{0}^{*}(2300)$. The fit results are shown in Fig. 3, where the blue and
green bands correspond to the best fits from UChPT and the BW, respectively.
While UChPT describes the data very well with $\chi^{2}/\text{d.o.f.}=1.2$,
the BW fails to reproduce the data with $\chi^{2}/\text{d.o.f.}=2.0$. The
difference of these two values is significant from the statistical point of
view: the corresponding $p$-values are 0.1 for the UChPT fit and $3\times
10^{-5}$ for the BW fit, respectively. Thus, the former can be accepted as a
good description of the data, while the latter is highly disfavored Press et
al. (1992). The error bands correspond to the $1\sigma$ uncertainties
propagated from the input phases. The borders of the band of $\langle
P_{13}\rangle$ for the BW are plotted in dotted and dashed curves to make it
evident that the data for $\langle P_{0}\rangle$ and $\langle P_{13}\rangle$
cannot be described by the BW phase. With the fitted parameters $g_{0}$ and
$g_{1}$, we obtain the $S$-wave phase of the decay amplitude for $B^{-}\to
D^{+}\pi^{-}\pi^{-}$ shown in Fig. 4, where the results corresponding to UChPT
and the BW are plotted as blue and green bands, respectively. As expected,
UChPT describes the $S$-wave phase extracted using Eq. (3) and that obtained
in Ref. Aaij et al. (2016) well up to 2.4 GeV. For the BW one, although the
error band is broad, either the low-energy or the high-energy region cannot be
described.

Conclusion.—The existence of the $D_{0}^{*}(2300)$ as given in the RPP is the
starting point of many theoretical analyses (see, e.g., Refs. Mehen and
Springer (2005); Colangelo et al. (2012); Alhakami (2016); Cheng and Yu
(2017)). The results obtained in this Letter show that the $D_{0}^{*}(2300)$,
whose resonance parameters were obtained using the BW parametrization from the
Belle Abe et al. (2004) and BaBar Aubert et al. (2009) analyses, is in
conflict with the much more precise LHCb data for $B^{-}\to
D^{+}\pi^{-}\pi^{-}$, which, however, can be well reproduced by the UChPT
amplitude containing the $D_{0}^{*}(2100)$.

We expect that the $D_{1}(2430)$ as given in the RPP Zyla et al. (2020) will
also be in conflict with future high-quality data of $B^{-}\to
D^{*+}\pi^{-}\pi^{-}$ from LHCb Aaij et al. (2020) and Belle-II, and that the
lightest $D_{1}$ meson is the $D_{1}(2250)$ predicted by UChPT Albaladejo et
al. (2017); Du et al. (2018).

The $D_{0}^{*}$ is analogous to the more famous $f_{0}(500)$ and
$K_{0}^{*}(700)$, whose masses have been significantly shifted from earlier
versions of the RPP due to improved data and improved theoretical analyses—for
recent discussions see Refs. Caprini et al. (2006); Descotes-Genon and
Moussallam (2006); Peláez (2016); Peláez and Rodas (2020a, b) and the review
on scalar mesons in the RPP Zyla et al. (2020). We expect a similar change in
all systems emerging from the scattering of a pion off an isospin-nonsinglet
hadron. The lightest resonance in that case should not be extracted from data
using the usual BW form—a parametrization accounting for chiral symmetry and
coupled channels is mandatory. The $D\pi$ $S$-wave phase extracted model-
independently here provides valuable information for further understanding
matter-field–Goldstone-boson scattering and the structure of positive-parity
heavy hadrons.

Figure 4: Phases obtained from Eq. (Where is the lightest charmed scalar
meson?) with scattering phase shifts from UChPT (blue) and BW (green).

Note added.—Recently, a lattice calculation also concluded that the
$D_{0}^{*}$ mass should be lower than the RPP value Gayer et al. (2021). The
authors found a mass of $(2196\pm 64)\,\mbox{MeV}$ with a pion mass of
$239\,\mbox{MeV}$, only $(77\pm 64)\,\mbox{MeV}$ above the $D\pi$ threshold.
Thus, our conclusion receives a strong support from lattice QCD calculations.

###### Acknowledgements.

This work is supported in part by the National Natural Science Foundation of
China (NSFC) under Grants No. 11835015, No. 12047503, and No. 11961141012, by
the NSFC and the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) through the funds provided to the Sino-German Collaborative
Research Center “Symmetries and the Emergence of Structure in QCD” (NSFC Grant
No. 12070131001, DFG Project-ID 196253076 – TRR110), by the Chinese Academy of
Sciences (CAS) under Grants No. XDB34030000 and No. QYZDB-SSW-SYS013, and by
the CAS Center for Excellence in Particle Physics (CCEPP). The work of U.G.M.
was also supported by the Chinese Academy of Sciences (CAS) President’s
International Fellowship Initiative (PIFI) (Grant No. 2018DM0034), by
VolkswagenStiftung (Grant No. 93562), and by the EU (Strong2020).

## References

  * Aubert et al. (2003) B. Aubert et al. (BaBar), Phys. Rev. Lett. 90, 242001 (2003), eprint hep-ex/0304021.
  * Besson et al. (2003) D. Besson et al. (CLEO), Phys. Rev. D 68, 032002 (2003), [Erratum: Phys. Rev. D 75, 119908 (2007)], eprint hep-ex/0305100.
  * Godfrey and Isgur (1985) S. Godfrey and N. Isgur, Phys. Rev. D 32, 189 (1985).
  * Godfrey and Moats (2016) S. Godfrey and K. Moats, Phys. Rev. D 93, 034035 (2016), eprint 1510.08305.
  * Ebert et al. (2010) D. Ebert, R. N. Faustov, and V. O. Galkin, Eur. Phys. J. C 66, 197 (2010), eprint 0910.5612.
  * Barnes et al. (2003) T. Barnes, F. E. Close, and H. J. Lipkin, Phys. Rev. D 68, 054006 (2003), eprint hep-ph/0305025.
  * van Beveren and Rupp (2003) E. van Beveren and G. Rupp, Phys. Rev. Lett. 91, 012003 (2003), eprint hep-ph/0305035.
  * Szczepaniak (2003) A. P. Szczepaniak, Phys. Lett. B 567, 23 (2003), eprint hep-ph/0305060.
  * Kolomeitsev and Lutz (2004) E. E. Kolomeitsev and M. F. M. Lutz, Phys. Lett. B 582, 39 (2004), eprint hep-ph/0307133.
  * Chen and Li (2004) Y.-Q. Chen and X.-Q. Li, Phys. Rev. Lett. 93, 232001 (2004), eprint hep-ph/0407062.
  * Guo et al. (2006) F.-K. Guo, P.-N. Shen, H.-C. Chiang, R.-G. Ping, and B.-S. Zou, Phys. Lett. B 641, 278 (2006), eprint hep-ph/0603072.
  * Guo et al. (2007) F.-K. Guo, P.-N. Shen, and H.-C. Chiang, Phys. Lett. B 647, 133 (2007), eprint hep-ph/0610008.
  * Gamermann et al. (2007) D. Gamermann, E. Oset, D. Strottman, and M. J. Vicente Vacas, Phys. Rev. D 76, 074016 (2007), eprint hep-ph/0612179.
  * Cheng and Hou (2003) H.-Y. Cheng and W.-S. Hou, Phys. Lett. B 566, 193 (2003), eprint hep-ph/0305038.
  * Maiani et al. (2005) L. Maiani, F. Piccinini, A. D. Polosa, and V. Riquer, Phys. Rev. D 71, 014028 (2005), eprint hep-ph/0412098.
  * Browder et al. (2004) T. E. Browder, S. Pakvasa, and A. A. Petrov, Phys. Lett. B 578, 365 (2004), eprint hep-ph/0307054.
  * Abe et al. (2004) K. Abe et al. (Belle), Phys. Rev. D 69, 112002 (2004), eprint hep-ex/0307021.
  * Link et al. (2004) J. M. Link et al. (FOCUS), Phys. Lett. B 586, 11 (2004), eprint hep-ex/0312060.
  * Liu et al. (2013) L. Liu, K. Orginos, F.-K. Guo, C. Hanhart, and U.-G. Meißner, Phys. Rev. D 87, 014508 (2013), eprint 1208.4535.
  * Mohler et al. (2013) D. Mohler, C. B. Lang, L. Leskovec, S. Prelovsek, and R. M. Woloshyn, Phys. Rev. Lett. 111, 222001 (2013), eprint 1308.3175.
  * Lang et al. (2014) C. B. Lang, L. Leskovec, D. Mohler, S. Prelovsek, and R. M. Woloshyn, Phys. Rev. D 90, 034510 (2014), eprint 1403.8103.
  * Moir et al. (2016) G. Moir, M. Peardon, S. M. Ryan, C. E. Thomas, and D. J. Wilson, JHEP 10, 011 (2016), eprint 1607.07093.
  * Bali et al. (2017) G. S. Bali, S. Collins, A. Cox, and A. Schäfer, Phys. Rev. D 96, 074501 (2017), eprint 1706.01247.
  * Cheung et al. (2020) G. K. C. Cheung, C. E. Thomas, D. J. Wilson, G. Moir, M. Peardon, and S. M. Ryan (2020), eprint 2008.06432.
  * Aaij et al. (2016) R. Aaij et al. (LHCb), Phys. Rev. D 94, 072001 (2016), eprint 1608.01289.
  * Albaladejo et al. (2017) M. Albaladejo, P. Fernandez-Soler, F.-K. Guo, and J. Nieves, Phys. Lett. B 767, 465 (2017), eprint 1610.06727.
  * Du et al. (2018) M.-L. Du, M. Albaladejo, P. Fernández-Soler, F.-K. Guo, C. Hanhart, U.-G. Meißner, J. Nieves, and D.-L. Yao, Phys. Rev. D 98, 094018 (2018), eprint 1712.07957.
  * Guo et al. (2018a) F.-K. Guo, C. Hanhart, U.-G. Meißner, Q. Wang, Q. Zhao, and B.-S. Zou, Rev. Mod. Phys. 90, 015004 (2018a), eprint 1705.00141.
  * Du et al. (2017) M.-L. Du, F.-K. Guo, U.-G. Meißner, and D.-L. Yao, Eur. Phys. J. C 77, 728 (2017), eprint 1703.10836.
  * Guo et al. (2009) F.-K. Guo, C. Hanhart, and U.-G. Meißner, Eur. Phys. J. A 40, 171 (2009), eprint 0901.1597.
  * Guo et al. (2015) Z.-H. Guo, U.-G. Meißner, and D.-L. Yao, Phys. Rev. D 92, 094008 (2015), eprint 1507.03123.
  * Guo et al. (2019) Z.-H. Guo, L. Liu, U.-G. Meißner, J. A. Oller, and A. Rusetsky, Eur. Phys. J. C 79, 13 (2019), eprint 1811.05585.
  * Guo et al. (2018b) X.-Y. Guo, Y. Heo, and M. F. M. Lutz, PoS LATTICE2018, 085 (2018b), eprint 1811.00478.
  * Guo et al. (2018c) X.-Y. Guo, Y. Heo, and M. F. M. Lutz, Phys. Rev. D 98, 014510 (2018c), eprint 1801.10122.
  * Aaij et al. (2014) R. Aaij et al. (LHCb), Phys. Rev. D 90, 072003 (2014), eprint 1407.7712.
  * Aaij et al. (2015a) R. Aaij et al. (LHCb), Phys. Rev. D 92, 032002 (2015a), eprint 1505.01710.
  * Aaij et al. (2015b) R. Aaij et al. (LHCb), Phys. Rev. D 91, 092002 (2015b), [Erratum: Phys. Rev. D 93, 119901 (2016)], eprint 1503.02995.
  * Aaij et al. (2015c) R. Aaij et al. (LHCb), Phys. Rev. D 92, 012012 (2015c), eprint 1505.01505.
  * Du et al. (2019) M.-L. Du, F.-K. Guo, and U.-G. Meißner, Phys. Rev. D 99, 114002 (2019), eprint 1903.08516.
  * Albaladejo et al. (2018) M. Albaladejo, P. Fernandez-Soler, J. Nieves, and P. G. Ortega, Eur. Phys. J. C 78, 722 (2018), eprint 1805.07104.
  * Meißner (2020) U.-G. Meißner, Symmetry 12, 981 (2020), eprint 2005.06909.
  * Zyla et al. (2020) P. A. Zyla et al. (Particle Data Group), PTEP 2020, 083C01 (2020).
  * Rosner (2013) J. L. Rosner, Phys. Rev. D 88, 034034 (2013), eprint 1307.2550.
  * Guo (2019) F.-K. Guo, Phys. Rev. Lett. 122, 202002 (2019), eprint 1902.11221.
  * Lin and Chiang (1983) W.-L. Lin and C.-C. Chiang, Lett. Nuovo Cim. 38, 503 (1983).
  * Savage and Wise (1989) M. J. Savage and M. B. Wise, Phys. Rev. D 39, 3346 (1989), [Erratum: Phys. Rev. D 40, 3127 (1989)].
  * Niecknig and Kubis (2015) F. Niecknig and B. Kubis, JHEP 10, 142 (2015), eprint 1509.03188.
  * Szczepaniak (2015) A. P. Szczepaniak, Phys. Lett. B 747, 410 (2015), eprint 1501.01691.
  * Gardner and Meißner (2002) S. Gardner and U.-G. Meißner, Phys. Rev. D 65, 094004 (2002), eprint hep-ph/0112281.
  * Khuri and Treiman (1960) N. N. Khuri and S. B. Treiman, Phys. Rev. 119, 1115 (1960).
  * Bronzan and Kacser (1963) J. B. Bronzan and C. Kacser, Phys. Rev. 132, 2703 (1963).
  * Niecknig and Kubis (2018) F. Niecknig and B. Kubis, Phys. Lett. B 780, 471 (2018), eprint 1708.00446.
  * Omnès (1958) R. Omnès, Nuovo Cim. 8, 316 (1958).
  * Kacser (1963) C. Kacser, Phys. Rev. 132, 2712 (1963).
  * (55) See Supplemental Material in the appendix for technical details, parameter values from the fits, and a numerical demonstration of the effect of chiral symmetry on resonance parameters extracted using an isobar model.
  * Watson (1954) K. M. Watson, Phys. Rev. 95, 228 (1954).
  * Press et al. (1992) W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, _Numerical Recipes in FORTRAN: The Art of Scientific Computing_ (Cambridge University Press, Cambridge, 1992).
  * Mehen and Springer (2005) T. Mehen and R. P. Springer, Phys. Rev. D 72, 034006 (2005), eprint hep-ph/0503134.
  * Colangelo et al. (2012) P. Colangelo, F. De Fazio, F. Giannuzzi, and S. Nicotri, Phys. Rev. D 86, 054024 (2012), eprint 1207.6940.
  * Alhakami (2016) M. H. Alhakami, Phys. Rev. D 93, 094007 (2016), eprint 1603.08848.
  * Cheng and Yu (2017) H.-Y. Cheng and F.-S. Yu, Eur. Phys. J. C 77, 668 (2017), eprint 1704.01208.
  * Aubert et al. (2009) B. Aubert et al. (BaBar), Phys. Rev. D 79, 112004 (2009), eprint 0901.1291.
  * Aaij et al. (2020) R. Aaij et al. (LHCb), Phys. Rev. D 101, 032005 (2020), eprint 1911.05957.
  * Caprini et al. (2006) I. Caprini, G. Colangelo, and H. Leutwyler, Phys. Rev. Lett. 96, 132001 (2006), eprint hep-ph/0512364.
  * Descotes-Genon and Moussallam (2006) S. Descotes-Genon and B. Moussallam, Eur. Phys. J. C 48, 553 (2006), eprint hep-ph/0607133.
  * Peláez (2016) J. R. Peláez, Phys. Rept. 658, 1 (2016), eprint 1510.00653.
  * Peláez and Rodas (2020a) J. R. Peláez and A. Rodas, Phys. Rev. Lett. 124, 172001 (2020a), eprint 2001.08153.
  * Peláez and Rodas (2020b) J. R. Peláez and A. Rodas (2020b), eprint 2010.11222.
  * Gayer et al. (2021) L. Gayer, N. Lang, S. M. Ryan, D. Tims, C. E. Thomas, and D. J. Wilson (2021), eprint 2102.04973.

Supplemental Material

##  Appendix A Technical details

The inhomogeneity for the $S$-wave $s$-channel amplitude is

$\hat{\mathcal{F}}_{0}^{1/2}(s)=-\frac{1}{4\sqrt{2}}\int_{-1}^{1}dz_{s}(t-s)\mathcal{F}_{1}^{1}(u)=-\frac{s}{4\sqrt{2}\kappa(s)}\int_{u_{-}(s)}^{u_{+}(s)}du(\Sigma_{0}-u-2s)\mathcal{F}_{1}^{1}(u),$
(15)

where $\Sigma_{0}=M_{B}^{2}+M_{D}^{2}+2M_{\pi}^{2}$, and the integral end
points and path need to be understood with the prescription $M_{B}^{2}\to
M_{B}^{2}+i\epsilon$. The inhomogeneities give the contributions from the
crossed channels, i.e., the left-hand cut contributions. The integration end
points $u_{\pm}(s)$ contain the only nontrivial $M_{B}$-dependent pieces. The
continuation following the $M_{B}^{2}\to M_{B}^{2}+i\epsilon$ prescription
gives

$\displaystyle u_{+}(s)$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}\frac{\Sigma_{0}s-s^{2}-\Delta+|\kappa_{s}(s)|+i\epsilon}{2s},&s\in[(M_{D}+M_{\pi})^{2},(M_{B}-M_{\pi})^{2}),\\\
\frac{\Sigma_{0}s-s^{2}-\Delta+i|\kappa_{s}(s)|}{2s},&s\in[(M_{B}-M_{\pi})^{2},(M_{B}+M_{\pi})^{2}),\\\
\frac{\Sigma_{0}s-s^{2}-\Delta-|\kappa_{s}(s)|}{2s},&s\in[(M_{B}+M_{\pi})^{2},\infty),\end{array}\right.$
(19) $\displaystyle u_{-}(s)$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}\frac{\Sigma_{0}s-s^{2}-\Delta-|\kappa_{s}(s)|+i\epsilon}{2s},&s\in[(M_{D}+M_{\pi})^{2},(M_{B}^{2}+M_{D}^{2}-2M_{\pi}^{2})/2),\\\
\frac{\Sigma_{0}s-s^{2}-\Delta-|\kappa_{s}(s)|-i\epsilon}{2s},&s\in\Big{[}(M_{B}^{2}+M_{D}^{2}-2M_{\pi}^{2})/2,(M_{B}-M_{\pi})^{2}\Big{)},\\\
\frac{\Sigma_{0}s-s^{2}-\Delta-i|\kappa_{s}(s)|}{2s},&s\in[(M_{B}-M_{\pi})^{2},(M_{B}+M_{\pi})^{2}),\\\
\frac{\Sigma_{0}s-s^{2}-\Delta+|\kappa_{s}(s)|}{2s},&s\in[(M_{B}+M_{\pi})^{2},\infty).\end{array}\right.$
(24)

The amplitude $\mathcal{F}_{1}^{1}(u)$ is parameterized by a Breit-Wigner (BW)
amplitude of the $\rho^{-}$ meson, since, as explained in the main text, the
process $B^{-}\to D^{0}\pi^{0}\pi^{-}$ is assumed to proceed through the
intermediate decay $B^{-}\to D^{0}\rho^{-}$,

$\displaystyle\mathcal{F}_{1}^{1}(u)=c_{S}R_{1}(u)\times
X_{1}(|\vec{p}_{u}|r_{\text{BW}})\times X_{1}(|\vec{q}_{u}|r_{\text{BW}}),$
(25)

where the $\vec{p}_{u}$ and $\vec{q}_{u}$ are the momenta of the $D^{0}$ and
$\pi^{+}$ ($\pi^{0}$) in the $\pi^{0}\pi^{-}$ rest frame, respectively. The
$\rho$ resonance lineshape is modelled by a relativistic BW function
$R_{1}(u)$ given by

$\displaystyle R_{1}(u)=\frac{1}{M_{\rho}^{2}-u-iM_{\rho}\Gamma({u})},$ (26)

with a mass-dependent decay width defined as

$\displaystyle\Gamma({u})=\Gamma_{\rho}\Big{(}\frac{|\vec{q}_{u}|}{|\vec{q}_{u}^{0}|}\Big{)}^{3}\frac{M_{\rho}}{\sqrt{u}}X_{1}^{2}(|\vec{q}_{u}|r_{\text{BW}}),$
(27)

where $q_{u}^{0}$ is the value of $q_{u}$ when $\sqrt{u}=M_{\rho}$, and
$\Gamma_{\rho}$ is the full width. The Blatt-Weisskopf barrier factors
$X_{L}(x)$ are given by

$X_{0}(x)=1,\qquad X_{1}(x)=\sqrt{\frac{x_{0}^{2}+1}{x^{2}+1}},\qquad
X_{2}(x)=\sqrt{\frac{x_{0}^{4}+3x_{0}^{2}+9}{x^{4}+3x^{2}+9}},$ (28)

where $x=|\vec{q}|r_{\text{BW}}$ or $|\vec{p}|r_{\text{BW}}$, and $x_{0}$ is
the value of $x$ evaluated at $\sqrt{u}=M_{\rho}$. The barrier radius,
$r_{\text{BW}}$, is taken to be $4.0\,\mbox{GeV}^{-1}$ to be consistent with
Ref. Aaij et al. (2016).

For the $P$-($D$-)wave amplitudes $\mathcal{F}_{1,2}^{1/2}(s)$, we use the
same BW amplitudes as those in the LHCb analysis Aaij et al. (2016), i.e.,

$\displaystyle\mathcal{F}_{J}^{1/2}(s)$ $\displaystyle=$ $\displaystyle
c_{J}R_{J}(s)\times X_{J}(|\vec{p}|r_{\text{BW}})\times
X_{J}(|\vec{q}|r_{\text{BW}})(-s)^{J},$ (29)

where $\vec{p}$ and $\vec{q}$ are the momenta of the $\pi^{-}$($p_{2}$) and
the $\pi$($p_{1}$), calculated in the $D\pi(p_{1})$ rest frame. $R_{J}(s)$ is
given by

$R_{J}(s)=\frac{1}{m_{J}^{2}-s-im_{J}\Gamma_{J}({s})},\qquad\Gamma_{J}({s})=\Gamma_{J}^{0}\Big{(}\frac{|{\vec{q}}|}{|\vec{q}_{0}|}\Big{)}^{2J+1}\Big{(}\frac{m_{J}}{\sqrt{s}}\Big{)}X_{J}^{2}(|\vec{q}|r_{\text{BW}}),$
(30)

where $\vec{q}_{0}$ is the value of $\vec{q}$ evaluated at the resonance mass
$\sqrt{s}=m_{J}$ with $m_{J}$ the mass of the $P$\- or $D$-wave resonance and
$\Gamma_{J}^{0}$ the corresponding width. The $P$-wave resonance $D^{*}$ is
outside the kinematically allowed region. Therefore, its mass $m_{1}$ is
modified to an effective mass $m_{1}^{\text{eff}}$ as in the LHCb analysis
Aaij et al. (2016),

$\displaystyle
m_{1}^{\text{eff}}=m^{\text{min}}+(m^{\text{max}}-m^{\text{min}})\Bigg{(}1+\text{tanh}\frac{m_{1}-\frac{m^{\text{min}}+m^{\text{max}}}{2}}{m^{\text{max}}-m^{\text{min}}}\Bigg{)},$
(31)

where $m^{\text{max}}$ and $m^{\text{min}}$ represent the upper and lower
limits of the kinematically allowed range.

##  Appendix B Parameters

The masses and widths of the $P$\- and $D$-wave resonances used in the
calculation are listed in Table 1 and are the same as those used in the LHCb
analysis Aaij et al. (2016). The UChPT amplitudes for the $D\pi$ interaction
and corresponding low energy constants (LECs) used in this work are taken from
Ref. Liu et al. (2013). For the BW parameterization of the $S$-wave $D\pi$
amplitude, the amplitude of Eq. (29) with $J=0$ is employed.

Table 1: Masses and widths of the vector and tensor charmed mesons as used in Ref. Aaij et al. (2016). Contribution | Spin | Mass (MeV) | Width (MeV)
---|---|---|---
$D^{\ast}(2007)^{0}$ | 1 | $2006.98$ | 2.1
$D_{2}^{\ast}(2460)^{0}$ | 2 | $2463.7$ | 47.0

The parameters determined in the best fits to data for the Khuri-Treiman
formalism, e.g., Eq. (Where is the lightest charmed scalar meson?), with input
from UChPT and the BW amplitudes, respectively, are collected in Table 2,
where only the statistical uncertainties are given. Each partial wave has an
overall normalization factor, e.g., $c_{1}$ and $c_{2}$ for $P$\- and
$D$-waves, respectively. Since we fit to the unnormalized angular moments, the
units of the fitted parameters contain some arbitrary overall normalization
factors.

Table 2: Parameters determined in the best fit to data for the Khuri-Treiman formalism with input from UChPT and the BW amplitudes, respectively, with statistical uncertainties only. Here, the $g_{i}=a_{i}e^{i\delta_{i}}$ refer to the subtraction constants in Eq. (Where is the lightest charmed scalar meson?), with $a_{i}$ real numbers. Parameters | UChPT | BW
---|---|---
$10^{3}\cdot a_{0}$ | $0.57\pm 0.14$ | $0.59\pm 0.25$
$\delta_{0}$ | $-0.38\pm 0.08$ | $-2.57\pm 0.05$
$10^{3}\cdot a_{1}$ | $1.15\pm 0.75$ | $0.95\pm 0.38$
$\delta_{1}$ | $1.56\pm 0.40$ | $-0.15\pm 0.05$
$10^{3}\cdot c_{S}$ | $1.36\pm 0.46$ | $0.00\pm 0.00$
$10^{4}\cdot c_{1}$ | $1.65\pm 0.06$ | $1.56\pm 0.05$
$10^{4}\cdot c_{2}$ [GeV-2] | $1.47\pm 0.07$ | $1.64\pm 0.06$

##  Appendix C Describing the data in a two-body isobar model

To make it explicit that the requirement of chiral symmetry has a sizeable
impact on the determination of resonance parameters using the BW
parameterization in the isobar model as argued in Ref. Du et al. (2019), here
we show two different fits:

  * •

model I: the $S$-wave $D\pi$ is parameterized as a relativistic BW amplitude,
Eq. (29);

  * •

model II: the $S$-wave $D\pi$ is parameterized as a BW amplitude times an
$E_{\pi}$ factor, as required by the chiral symmetry of QCD:

$\displaystyle\mathcal{F}_{0}^{\prime}(s)=c_{0}R_{0}(s)\times 2E_{\pi}/M_{D},$
(32)

with $E_{\pi}$ the energy of the pion in the $D\pi$ rest frame.

Note that in the isobar model, the crossed-channel effects are not explicitly
considered and the coefficients $c_{J}$ are complex to describe the relative
contribution of each intermediate process. Here, we do not use the resonance
parameters of $D_{0}^{*}$ in the Review of Particle Physics (RPP), but
determine the BW mass and width by fitting to the high quality LHCb data in
Ref. Aaij et al. (2016). The best fits to the angular moments are shown in
Fig. 5, where the orange (with $\chi^{2}/\text{d.o.f.}=1.8$) and the blue
curves ($\chi^{2}/\text{d.o.f.}=1.1$) correspond to model I and model II,
respectively. For model I, the fit results in the $D_{0}^{*}$ BW mass and
width parameters of $(2319\pm 9)\,\mbox{MeV}$ and $(420\pm 25)\,\mbox{MeV}$,
respectively. The fitted mass is consistent with the RPP average mass of
$(2300\pm 19)$ MeV. However, the width is much larger than the RPP average
value, i.e., $(274\pm 40)\,\mbox{MeV}$. Note that in Ref. Aaij et al. (2016),
a quasi-model-independent approach is used to describe the $D^{+}\pi^{-}$
$S$-wave. The $D^{+}\pi^{-}$ $S$-wave is extracted using cubic splines to
describe the magnitude and phase variation. Accordingly, in this reference no
pole parameters for the $D_{0}^{*}$ are given. When the $E_{\pi}$ factor is
considered, i.e., in model II, the fit quality is significantly improved with
a much smaller $\chi^{2}$. The fitted parameters are given in Table 3 with
only statistical uncertainties. In addition, the fit results in a BW mass of
$(2206\pm 4)\,\mbox{MeV}$, which is much lower than that given in the RPP,
while the width is $(341\pm 17)\,\mbox{MeV}$. This clearly shows that the BW
parameters for a broad resonance can be affected significantly by the chiral
symmetry constraint, as argued in Ref. Du et al. (2019). Nevertheless, the
modification of the BW in model II only applies in a small energy region
before the coupled-channel effects become important, and thus is neither
practical nor systematic. The resonance parameters would get further modified
after the $D\eta$ and $D_{s}\bar{K}$ coupled channels are taken into account.
The UChPT formalism used in the main text provides a theoretical framework
satisfying both chiral symmetry and coupled-channel unitarity.

Figure 5: The best fit results of model I (orange) and model II (blue) in the isobar model, with $\chi^{2}/{\rm d.o.f.}=1.8$ and 1.1, respectively. Table 3: Parameters determined by the fit to data in the two-body isobar model for model I and II with statistical uncertainties only, where $c_{i}=|c_{i}|e^{i\phi_{i}}$ with $\phi_{i}$ real numbers and the phase of $D$-wave is fixed to zero, i.e., $\phi_{2}=0$.  Parameters | model I | mdoel II
---|---|---
$10^{3}\cdot|c_{0}|$ [GeV2] | $3.39\pm 0.18$ | $7.28\pm 0.30$
$\phi_{0}$ | $1.08\pm 0.17$ | $0.98\pm 0.21$
$10^{4}\cdot|c_{1}|$ | $1.60\pm 0.05$ | $1.67\pm 0.05$
$\phi_{1}$ | $-0.08\pm 0.16$ | $0.30\pm 0.21$
$10^{4}\cdot|c_{2}|$ [GeV-2] | $1.71\pm 0.11$ | $1.43\pm 0.07$
$m_{\text{BW}}$ [GeV] | $2.319\pm 0.009$ | $2.206\pm 0.004$
$\Gamma_{\text{BW}}^{0}$ [GeV] | $0.420\pm 0.025$ | $0.341\pm 0.017$


frame=single, rulecolor=, numberstyle=, commentstyle=, basicstyle=,
keywordstyle=, showstringspaces=false, numbers=left, numbersep=8pt,
xleftmargin=1.95em, framexleftmargin=1.6em

# SUPPLEMENTAL MATERIAL TO: THZ FIELD-INDUCED SPIN DYNAMICS IN FERRIMAGNETIC
IRON GARNETS

## 1 Static magneto-optical characterization of TmIG sample

Figure 1: Measurements of the magneto-optical Faraday rotation using a
continuous-wave helium-neon laser ($\lambda=632.8$ nm) with both the external
field and the light’s wave-vector perpendicular to the sample plane. A
paramagnetic contribution $\sim 0.56$ deg/T from the cryostat glass windows
has been subtracted from the raw data. The data exhibits a large rotation and
demonstrates a weak easy-axis type of anisotropy with small coercive $<6$ mT
and saturation $<25$ mT field. No compensation point was observed above
nitrogen temperatures $>77$ K. Figure 2: Static polarization rotation
measurements with light at the wavelength of $800$ nm in the experimental
geometry (see Fig. 1(a) in article). The evolution of hysteresis loop form can
be due to temperature dependent anisotropy constants. Clearly, no
magnetization compensation point is observed in this temperature range. At T
$=6$ K, the saturated polarization rotation is $\sim\pm 1.65^{\circ}$, and
given the angle of the magnetic field of $10^{\circ}$, this has been used to
estimate the Faraday rotation ($\sim\pm 10^{\circ})$ when the magnetization is
along the sample normal.

Figure 3: Domain pattern seen by magneto-optical microscopy in transmission at
zero field. The typical ”labyrinth” type domains grow with decreasing
temperature, which indicates a growing role of easy axis anisotropy and
thulium magnetization [1]. When applying an external magnetic field along the
out-of-plane easy axis, the domains along this field expand and the sample
will be uniformly magnetized for relatively small field (see Suppl. Fig. 1)

## 2 Experimental setup

The experimental setup regarding THz generation by optical rectification in
lithium niobate is described in detail in [2, 3]. The THz path was purged with
nitrogen to avoid water absorption lines in the THz spectrum. A small part of
the initial $800$ nm beam is chopped out beforehand (ratio $1:100$) and is
brought to spectral and temporal overlap with the THz pump pulse. The focused
spot size of the probe beam is considerably smaller than that of the THz. The
waveform of the THz pulse was mapped using electro-optical sampling in a $50$
$\mu$m GaP [$110$] crystal seen in Supplemental Fig. 4.

Figure 4: THz waveform and corresponding Fourier spectrum measured by EO
sampling in GaP.

## 3 Supplemental Results

Figure 5: THz induced polarization rotation waveforms for two orthogonal THz
pump polarizations (two figures) and for several orientations of the probe
polarization. The angle depicted is the angle of the probe electric field with
respect to the experimental $x$-axis (Fig. 1(a) of the article). This data
implies that the THz induced signals are Faraday rotation (see main text).
Figure 6: Peak-to-peak amplitudes of THz induced waveforms as a function of
external magnetic field (a) and THz field (b) for two orthogonal THz pump
polarizations. Bending of the red dots at low THz fields is attributed to the
facts that the THz light is not perfectly linearly polarized and imperfections
of the wire grid polarizers.

Figure 7: THz induced Faraday rotation waveforms obtained at several angles of
the pump polarization angle $\alpha$ and applied external field of $\pm 250$
mT. Besides the previous figure, this is the only graph where a weaker THz
electric field of $160$ kV/cm has been used. The result shows that measuring
at $\alpha=\pm 45^{\circ}$, in between the fully symmetric $\alpha=0^{\circ}$
and fully antisymmetric orthogonal $\alpha=\pm 90^{\circ}$, results in a mix
of symmetry/antisymmetry. Moreover, the effects gradually become weaker
towards $\alpha=0^{\circ}$
($\mathbf{H}_{THz}\parallel\mathbf{M}_{\parallel}$).

Figure 8: FFT spectra of THz induced Faraday rotation for
$\mathbf{H}_{THz}\perp\mathbf{M}_{\parallel}$. The exchange-mode frequency at
about $375$ GHz shows softening similar to the case
$\mathbf{H}_{THz}\parallel\mathbf{M}_{\parallel}$ as is presented in article
Fig. 2. At lower temperatures another high frequency ($725$ GHz) appears. It
is known that crystal field transition may appear in this region [4], but if
it can be attributed to these transitions is yet unclear.  Figure 9:
Experimental waveforms of THz induced Faraday rotation as a function of
temperature. In both cases the same external fields (specified in the first
figure) have been applied to ensure saturation of static magnetization.

Figure 10: Preliminary data of THz-induced ferromagnetic resonance, used to
estimate the effective $g$-factor $g_{eff}\approx 6$. Figure 11: The dotted
line shows the experimentally calibrated THz magnetic field pulse, which has
been fitted using the Gaussian derivative function
$G^{\prime}(x)=-2A((x-d)/w)\exp\left[((x-d)/w)^{2}\right]$ with $A=404$ mT,
$d=1.17$ ps (variable, determines arrival time of pulse) and $w=0.2223$ ps
pulse-width.

Figure 12: Simulated dynamics of the iron magnetization $\mathbf{M}_{Fe}$
using LLG equations and plotted separately for the $x$, $y$ and $z$ components
where $z$ coincides with the sample out-of-plane axis. It shows how the
symmetry with respect to external field is exactly opposite when looking at
the $y$-component, to which we are not experimentally sensitive.

## 4 Equations of motion derived from Lagrangian formalism

We start from the following Lagrangian and Rayleigh dissipation functions,
which are equivalent to the LLG equations for a two-sublattice ferrimagnet
[5]:

$\displaystyle\mathcal{L}$ $\displaystyle=$ $\displaystyle T-\Phi$ (1)
$\displaystyle=$
$\displaystyle-\frac{M_{Fe}}{\gamma_{Fe}}\cos\theta_{Fe}\frac{\partial\phi_{Fe}}{\partial
t}-\frac{M_{R}}{\gamma_{R}}\cos\theta_{R}\frac{\partial\phi_{R}}{\partial
t}-\Phi$ $\displaystyle\mathcal{R}$ $\displaystyle=$
$\displaystyle\mathcal{R}_{Fe}+\mathcal{R}_{R},\text{ \ \ \ \ \ \ \ \
}\mathcal{R}_{Fe,R}=\frac{\alpha
M_{Fe,R}}{2\gamma_{Fe,R}}\big{(}\dot{\theta}^{2}_{Fe,R}+\sin^{2}\theta_{Fe,R}\dot{\phi}^{2}_{Fe,R}\big{)},$
(2)

where $\theta_{i}$ and $\phi_{i}$ the polar and azimuthal angles of the iron
(Fe) and rare-earth (R) sublattices in the experimental coordinate system with
the $x$-axis aligned to the external magnetic field $\mathbf{H}_{ext}$ (see
Fig. 13, here we ignore the $10^{\circ}$ inclination of the field for
simplification).

Figure 13: Coordinate system used for Lagrangian equation.

The thermodynamic potential used is:

$\displaystyle\Phi$ $\displaystyle=$
$\displaystyle-(\mathbf{M}_{Fe}+\mathbf{M}_{R})\cdot\mathbf{H}_{ef}-\Lambda\mathbf{M}_{Fe}\cdot\mathbf{M}_{R}-K_{Fe}\frac{(\mathbf{M}_{Fe}\cdot\mathbf{n})^{2}}{M_{Fe}^{2}}-K_{R}\frac{(\mathbf{M}_{R}\cdot\mathbf{n})^{2}}{M_{R}^{2}}.$
(3)

Here $\mathbf{n}=(0,0,1)$ is the directional vector of the easy axis of
anisotropy, $\Lambda<0$ the intersublattice exchange constant and $K_{Fe,R}>0$
the uniaxial anisotropy constants. The Euler-Lagrange equations w.r.t.
$\theta_{i}$, $\phi_{i}$ give rise to four (coupled) equations of motion (two
for each sublattice), and this is generally not easy and sometimes even
impossible to solve. Instead, in Ref. [5] an effective Lagrangian is obtained
by assuming the canting of the two sublattices are equal and are assumed to be
small. This approach generally works at field well below the exchange field
(small canting), and it is valid here as the static measurements indicate we
are well below the spin-flop field.

We introduce the usual definitions of the magnetization
$\mathbf{M}=\mathbf{M}_{Fe}+\mathbf{M}_{R}$ and antiferromagnetic (Néel)
vector $\mathbf{L}=\mathbf{M}_{Fe}-\mathbf{M}_{R}$. These two vectors are
parameterized using a set of angles $\theta,\epsilon$ and $\phi,\beta$ defined
as:

$\displaystyle\theta_{Fe}$ $\displaystyle=$
$\displaystyle\theta-\epsilon,\quad\theta_{R}=\pi-\theta-\epsilon,$ (4)
$\displaystyle\phi_{Fe}$ $\displaystyle=$
$\displaystyle\phi+\beta,\quad\phi_{R}=\pi+\phi-\beta.$ (5)

In the quasi-antiferromagnetic approximation [5], the canting angles are
assumed to be small $\epsilon\ll 1$, $\beta\ll 1$. In first order
approximation, the $\mathbf{M}$ and $\mathbf{L}$ are then naturally defined
as:

$\displaystyle\mathbf{M}$ $\displaystyle=$ $\displaystyle
m(\cos\theta,\sin\theta\cos\phi,\sin\theta\sin\phi)$ (6)
$\displaystyle\mathbf{L}$ $\displaystyle=$
$\displaystyle\mathcal{M}(\cos\theta,\sin\theta\cos\phi,\sin\theta\sin\phi)$
(7)

where $m\equiv M_{Fe}-M_{R}$ and $\mathcal{M}\equiv M_{Fe}+M_{R}$.
Substituting our new set of angles (4)-(5) into the Lagrangian (1) and
expanding up to quadratic terms in the small variables $\epsilon,\beta$ gives
for the kinetic energy part:

$\displaystyle\mathcal{L}$ $\displaystyle=$
$\displaystyle-\frac{m}{\gamma_{ef}}\dot{\phi}\cos\theta-\frac{\mathcal{M}}{\overline{\gamma}}\sin\theta\left(\dot{\phi}\epsilon+\beta\dot{\theta}\right)-\Phi$
(8)

where we defined:

$\frac{1}{\gamma_{ef}}\equiv\frac{M_{Fe}/\gamma_{Fe}-M_{R}/\gamma_{R}}{M_{Fe}-M_{R}}\
\ \ \ \ \ \ \ \ \text{ and }\ \ \ \ \ \ \ \ \ \
\frac{1}{\overline{\gamma}}\equiv\frac{M_{Fe}/\gamma_{Fe}+M_{R}/\gamma_{R}}{M_{Fe}+M_{R}}.$
(9)

The potential energy $\Phi$ can be expanded similarly. Here, we make the
simplification that both sublattices experience the same effective anisotropy
$K_{U}\equiv(K_{Fe}+K_{R})/2$ in which case the anisotropy terms can by
replaced by a single term $-K_{U}(\mathbf{l}\cdot\mathbf{n})^{2}$ where
$\mathbf{l}=\mathbf{L}/|\mathbf{L}|$. Furthermore, the effective field
$\mathbf{H}_{ef}$ in (3) consists of the static external field and the time-
dependent THz magnetic field
$\mathbf{H}_{ef}=\mathbf{H}_{ext}+\mathbf{H}_{THz}$. The external field is
chosen along the $x$-axis $\mathbf{H}_{ext}=(H_{0},0,0)$, while we assume the
THz magnetic field lies in the $x-y$ plane
$\mathbf{H}_{THz}\equiv(h_{x}(t),h_{y}(t),0$) (see Fig. 1 from the article).
Writing $\delta\equiv-4\Lambda M_{Fe}M_{R}$ and $H\equiv H_{ext}+h_{x}$, the
potential energy becomes after expanding in $\epsilon,\beta$:

$\displaystyle\Phi$ $\displaystyle=$ $\displaystyle-
mH\cos\theta-\mathcal{M}H\epsilon\sin\theta-
mh_{y}\sin\theta\cos\phi+\mathcal{M}h_{y}\beta\sin\theta\sin\phi$
$\displaystyle+$ $\displaystyle\mathcal{M}h_{y}\epsilon\cos\theta\cos\phi-
mh_{y}\cos\theta\sin\phi\
\epsilon\cdot\beta+\frac{\delta}{2}\left(\epsilon^{2}+\beta^{2}\sin^{2}\theta\right)-K_{U}\sin^{2}\theta\sin^{2}\phi.$

We will ignore the term containing $\epsilon\cdot\beta$ as it is very small,
from the quadratic terms only the ones proportional to the exchange constant
$\sim\delta$ survive. We exclude the variables $\epsilon$, $\beta$ by solving
the Euler-Lagrange equations
$\frac{d}{dt}\frac{\partial\mathcal{L}}{\partial\dot{\epsilon}}-\frac{\partial\mathcal{L}}{\partial\epsilon}=-\frac{\partial\mathcal{R}}{\partial\dot{\epsilon}}\approx
0$ and
$\frac{d}{dt}\frac{\partial\mathcal{L}}{\partial\dot{\beta}}-\frac{\partial\mathcal{L}}{\partial\beta}=-\frac{\partial\mathcal{R}}{\partial\dot{\beta}}\approx
0$, giving:

$\displaystyle\epsilon$ $\displaystyle=$
$\displaystyle\frac{\mathcal{M}}{\delta}\sin\theta\left(H-\frac{\dot{\phi}}{\overline{\gamma}}\right)-\frac{\mathcal{M}h_{y}}{\delta}\cos\theta\cos\phi,$
(11) $\displaystyle\beta\sin\theta$ $\displaystyle=$
$\displaystyle-\frac{\mathcal{M}}{\delta}\left(\frac{\dot{\theta}}{\overline{\gamma}}+h_{y}\sin\phi\right).$
(12)

Substituting in (8)-(4) and rearranging terms yields the effective Lagrangian
from the article:

$\displaystyle\begin{aligned}
\mathcal{L}_{eff}&=\frac{\mathcal{M}^{2}}{2\delta}\Bigg{[}\left(\Big{(}\frac{\dot{\phi}}{\overline{\gamma}}-H\Big{)}\sin\theta+h_{y}\cos\theta\cos\phi\right)^{2}+\left(\frac{\dot{\theta}}{\overline{\gamma}}+h_{y}\sin\phi\right)^{2}\Bigg{]}+m\Big{(}H-\frac{\dot{\phi}}{\gamma_{ef}}\Big{)}\cos\theta\\\
&\quad+mh_{y}\sin\theta\cos\phi+K_{U}\sin^{2}\theta\sin^{2}\phi.\end{aligned}$
(13)

The equations of motion are now determined by Euler-Lagrange equations:

$\displaystyle\frac{d}{dt}\Big{(}\frac{\partial\mathcal{L}_{eff}}{\partial\dot{\theta}}\Big{)}-\frac{\partial\mathcal{L}_{eff}}{\partial\theta}+\frac{\partial\mathcal{R}}{\partial\dot{\theta}}=0,$
(14)
$\displaystyle\frac{d}{dt}\Big{(}\frac{\partial\mathcal{L}_{eff}}{\partial\dot{\phi}}\Big{)}-\frac{\partial\mathcal{L}_{eff}}{\partial\phi}+\frac{\partial\mathcal{R}}{\partial\dot{\phi}}=0.$
(15)

We solve these equations and linearize them around the equilibrium (ground-
state) equilibrium values $\theta_{0}$ and $\phi_{0}$, which are found by
minimizing (3) yielding $\phi_{0}=\pi/2$ and $\theta_{0}$ depending on the
ratio of external field to anisotropy (i.e. when $H_{ext}=0$ we have
$\theta_{0}=\pi/2$ while $\theta_{0}=0$ when $H_{ext}\gg H_{anis}$).
Linearizing around these values, i.e. $\theta=\theta_{0}+\theta_{l}$ and
$\phi=\phi_{0}+\phi_{l}$ with $\theta_{l},\phi_{l}\ll 1$, the first equation
(14) gives:

$\displaystyle\ddot{\theta}_{l}+\frac{\alpha\mathcal{M}\overline{\gamma}}{\chi_{\perp}}\dot{\theta}_{l}+\Big{(}-\overline{\gamma}^{2}H^{2}\cos
2\theta_{0}+\frac{m\overline{\gamma}^{2}H}{\chi_{\perp}}\cos\theta_{0}-\frac{2K_{U}\overline{\gamma}^{2}}{\chi_{\perp}}\cos
2\theta_{0}\Big{)}\theta_{l}+\Big{(}\overline{\gamma}H\sin
2\theta_{0}-\frac{m\overline{\gamma}^{2}}{\gamma_{ef}\chi_{\perp}}\sin\theta_{0}\Big{)}\dot{\phi}_{l}$
$\displaystyle+\Big{(}-\overline{\gamma}^{2}Hh_{y}\cos
2\theta_{0}+\frac{m\overline{\gamma}^{2}h_{y}}{\chi_{\perp}}\cos\theta_{0}\Big{)}\phi_{l}=-\overline{\gamma}\dot{h}_{y}+\frac{\overline{\gamma}^{2}H^{2}}{2}\sin
2\theta_{0}-\frac{m\overline{\gamma}^{2}H}{\chi_{\perp}}\sin\theta_{0}+\frac{\overline{\gamma}^{2}K_{U}}{\chi_{\perp}}\sin
2\theta_{0}.$ (16)

where we introduced the notation
$\chi_{\perp}\equiv\frac{\mathcal{M}^{2}}{\delta}$. Similarly for the second
Euler-Lagrange equation (15):

$\displaystyle\ddot{\phi}_{l}+\frac{\alpha\mathcal{M}\overline{\gamma}}{\chi_{\perp}}\dot{\phi}_{l}+\phi_{l}\Big{(}-\overline{\gamma}\dot{h}_{y}\cot\theta_{0}+\overline{\gamma}^{2}h_{y}^{2}+\frac{2K_{U}\overline{\gamma}^{2}}{\chi_{\perp}}\Big{)}+\dot{\theta}_{l}\Big{(}-2\overline{\gamma}H\cot\theta_{0}+\frac{m\overline{\gamma}^{2}}{\gamma_{ef}\chi_{\perp}\sin\theta_{0}}\Big{)}\
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ $ $\displaystyle+\
\theta_{l}\Big{(}-2\overline{\gamma}\dot{h}_{x}\cot\theta_{0}+\overline{\gamma}^{2}Hh_{y}(1-\cot^{2}\theta_{0})+\frac{\overline{\gamma}^{2}mh_{y}}{\chi_{\perp}}\frac{\cos\theta_{0}}{\sin^{2}\theta_{0}}\Big{)}=\overline{\gamma}\dot{h}_{x}+\overline{\gamma}^{2}Hh_{y}\cot\theta_{0}-\frac{m\overline{\gamma}^{2}h_{y}}{\chi_{\perp}}\frac{1}{\sin\theta_{0}}.$
(17)

These equations can be drastically simplified by noting that
$\frac{1}{\chi_{\perp}}$ is proportional the the exchange constant and is
therefore relatively large. Also the field derivative term $\gamma\dot{h}_{i}$
is strong, while terms proportional to $\sim\gamma h_{x,y}$ within brackets
are driving terms proportional to the response and thus negligible. Equations
(16)-(17) are then given in approximation:

$\ddot{\theta}_{l}+\frac{\alpha\mathcal{M}\overline{\gamma}}{\chi_{\perp}}\dot{\theta}_{l}+\Big{(}-\overline{\gamma}^{2}H^{2}\cos
2\theta_{0}+\frac{m\overline{\gamma}^{2}H}{\chi_{\perp}}\cos\theta_{0}-\frac{2K_{U}\overline{\gamma}^{2}}{\chi_{\perp}}\cos
2\theta_{0}\Big{)}\theta_{l}+\Big{(}\overline{\gamma}H\sin
2\theta_{0}-\frac{m\overline{\gamma}^{2}}{\gamma_{ef}\chi_{\perp}}\sin\theta_{0}\Big{)}\dot{\phi}_{l}\\\
=-\overline{\gamma}\dot{h}_{y}-\frac{m\overline{\gamma}^{2}H}{\chi_{\perp}}\sin\theta_{0}+\frac{\overline{\gamma}^{2}K_{U}}{\chi_{\perp}}\sin
2\theta_{0},$ (18)
$\ddot{\phi}_{l}+\frac{\alpha\mathcal{M}\overline{\gamma}}{\chi_{\perp}}\dot{\phi}_{l}+\frac{2K_{U}\overline{\gamma}^{2}}{\chi_{\perp}}\phi_{l}+\Big{(}-2\overline{\gamma}H\cot\theta_{0}+\frac{m\overline{\gamma}^{2}}{\gamma_{ef}\chi_{\perp}\sin\theta_{0}}\Big{)}\dot{\theta}_{l}=\overline{\gamma}\dot{h}_{x}+\overline{\gamma}^{2}Hh_{y}\cot\theta_{0}-\frac{m\overline{\gamma}^{2}h_{y}}{\chi_{\perp}}\frac{1}{\sin\theta_{0}}.$
(19)

The large field derivatives of the THz field $\dot{h}_{x,y}$ appear as a
dominant driving force in these equations of motion. Interestingly, only the
$y$-component $\dot{h}_{y}$ appears in the equation of motion for
$\theta_{l}$, which we use here to understand the qualitative difference in
dependencies on THz pump polarization assuming the field-derivative driving
force is dominant.

In the experiment we saturate the magnetization with the external field at a
small angle $\theta_{0}\approx 0$ (thus $\theta_{0}\approx\pi$ for
$\mathbf{-H}_{ext}$). Given that $\phi_{0}=\pi/2$, we have that the
modulations in the magnetization $z$-component are
$M_{z}(t)=M\sin\phi\sin\theta\sim\pm\theta_{l}$ for $\pm\mathbf{H}_{ext}$
external field. The experiment reveals we are only sensitive to $M_{z}(t)$, so
the detectable Faraday rotation modulations should also be proportional
$\sim\pm\theta_{l}(t)$. When $\mathbf{H}_{THz}\perp\mathbf{M}$,
$\dot{h}_{x}=h_{x}\neq 0$ and we have a strong non-zero driving force in (18),
it explains why we see immediate strong oscillations in $M_{z}$. Because the
driving term has the same sign for both external field polarities
$\pm\mathbf{H}_{ext}$, the forced oscillations must be sensitive to the
polarity of the external magnetic field
$\ddot{M}_{z}(t=0)=\frac{d^{2}}{dt^{2}}\sin\big(\theta_{0}+\theta_{l}(t)\big{missing})\Bigr{\rvert}_{t=0}\sim\pm\ddot{\theta}_{l}(t=0)\sim\mp\gamma\dot{h}_{y}$
(as $\theta_{0}=0,$ $\pi$ for $\pm\mathbf{H}_{ext}$) i.e. this is an
$\mathbf{H}$-odd effect. After the THz pulse has left the sample, the system
of equations resembles those for a harmonic oscillator in $2$D, meaning the
subsequent free oscillations will have opposite phases in the cases of
opposite polarities of the external magnetic field.

Meanwhile by a similar argument, it is clear why a strong response is absent
in $M_{z}$ when $\mathbf{H}_{THz}\parallel\mathbf{M}_{\parallel}$
($\dot{h}_{y}=h_{y}=0$) as in this case only the equation of motion for the
in-plane dynamics $\phi_{l}(t)$ (Eq. (17)), to which we are not sensitive, has
an initial non-zero driving force $\gamma\dot{h}_{x}$ while $\theta_{l}(t)$
does not. Detectable oscillations in $\theta_{l}(t)$ are instead only driven
by cross-terms like
$-\frac{m\overline{\gamma}^{2}}{\gamma_{ef}\chi_{\perp}}\sin\theta_{0}\dot{\phi}_{l}$
(Eq. (18)). Here it is important that the ground state $\theta_{0}$ is not
exactly equal to $0$ and $\pi$, i.e. $\sin(\theta_{0})=\pm\rho$ for
$\pm\mathbf{H}_{ext}$ with $\rho>0$ small constant (due to experimental
canting of external field, otherwise no dynamics in this case is observed as
was also seen in the experiment and simulations). Thus for opposite external
field polarities, the driving force in $\theta_{l}$ has opposite sign
$\mp\frac{m\gamma}{\chi_{\perp}}\rho\dot{\phi}_{l}$, contrary to the previous
case. This means that subsequent oscillations are now expected to be even with
$\mathbf{H}_{ext}$, in accordance with what was seen experimentally. Because
these field-even oscillations are a secondary result from primary in-plane
oscillations $\phi(t)$, it also explains why the observed effects are
relatively weak when $\mathbf{H}_{THz}\parallel\mathbf{M}_{\parallel}$
compared to $\mathbf{H}_{THz}\perp\mathbf{M}$.

The eigenfrequencies in the article have been found by solving the coupled set
of equations:

$\displaystyle\ddot{\theta}_{l}+\frac{2K_{U}\overline{\gamma}^{2}}{\chi_{\perp}}\theta_{l}-\frac{m\overline{\gamma}^{2}}{\gamma_{ef}\chi_{\perp}}\dot{\phi}_{l}$
$\displaystyle=$ $\displaystyle 0,$ (20)
$\displaystyle\ddot{\phi}_{l}+\frac{2K_{U}\overline{\gamma}^{2}}{\chi_{\perp}}\phi_{l}+\frac{m\overline{\gamma}^{2}}{\gamma_{ef}\chi_{\perp}}\dot{\theta}_{l}$
$\displaystyle=$ $\displaystyle 0.$ (21)

Assuming $\theta_{l},\phi_{l}\sim\exp(i\omega t)$ the frequencies can be
solved by the equation

$\begin{vmatrix}-\omega^{2}+\omega_{K}^{2}&-i\omega\omega_{ex}\\\
i\omega\omega_{ex}&-\omega^{2}+\omega_{K}^{2}\end{vmatrix}=0,$ (22)

where $\omega_{K}^{2}=2K_{U}\overline{\gamma}^{2}/\chi_{\perp}$ and
$\omega_{ex}=m\overline{\gamma}^{2}/\gamma_{ef}\chi_{\perp}$ in the case of
weak anisotropy $\omega_{ex}\gg\omega_{K}$, and we obtain:

$\omega=\pm\frac{\omega_{ex}}{2}\pm\sqrt{\frac{\omega_{ex}^{2}}{4}+\omega_{K}^{2}}\approx\pm\frac{\omega_{ex}}{2}\pm(\frac{\omega_{ex}}{2}+\omega_{K}^{2}/\omega_{ex}).$
(23)

Thus we obtain two approximate absolute frequencies:

$\displaystyle\omega_{1}$ $\displaystyle=$
$\displaystyle\omega_{K}^{2}/\omega_{ex}=\gamma_{ef}\frac{2K_{U}}{m},$ (24)
$\displaystyle\omega_{2}$ $\displaystyle\approx$
$\displaystyle\omega_{ex}=\frac{m\overline{\gamma}^{2}}{\gamma_{ef}\chi_{\perp}}\approx|\Lambda|(|\gamma_{R}|M_{Fe}-|\gamma_{Fe}|M_{R})$
(25)

where in the last approximation we used that
$M_{Fe}M_{R}\approx(M_{Fe}+M_{R})^{2}/4$ and
$(M_{Fe}/\gamma_{Fe}+M_{R}/\gamma_{R})^{2}\approx(M_{Fe}+M_{R})^{2}/\gamma_{Fe}\gamma_{R}$
to recover the approximate Kaplan-Kittel expression of exchange resonance.

## References

  * [1] A. Kalashnikova, V. Pavlov, A. Kimel, A. Kirilyuk, T. Rasing, and R. Pisarev, “Magneto-optical study of holmium iron garnet Ho3Fe5O12,” Low Temperature Physics, vol. 38, 09 2012.
  * [2] E. A. Mashkovich, K. A. Grishunin, R. V. Mikhaylovskiy, A. K. Zvezdin, R. V. Pisarev, M. B. Strugatsky, P. C. M. Christianen, T. Rasing, and A. V. Kimel, “Terahertz optomagnetism: Nonlinear THz excitation of GHz spin waves in antiferromagnetic ${\mathrm{FeBO}}_{3}$,” Phys. Rev. Lett., vol. 123, p. 157202, Oct 2019.
  * [3] M. Sajadi, M. Wolf, and T. Kampfrath, “Terahertz-field-induced optical birefringence in common window and substrate materials,” Opt. Express, vol. 23, pp. 28985–28992, Nov 2015.
  * [4] A. J. Sievers and M. Tinkham, “Far infrared spectra of rare-earth iron garnets,” Phys. Rev., vol. 129, pp. 1995–2004, Mar 1963.
  * [5] M. D. Davydova, K. A. Zvezdin, A. V. Kimel, and A. K. Zvezdin, “Ultrafast spin dynamics in ferrimagnets with compensation point,” Journal of Physics: Condensed Matter, vol. 32, p. 01LT01, oct 2019.