$\displaystyle\left.\hphantom{\mathcal{V}(y)=}{}+\int_{C(a,b)}\mathcal{Q}(\eta)\,\mathrm{d}\eta+\int_{b}^{y}\mathcal{Q}(\eta)\,\mathrm{d}\eta\right).$
(4.2)

Here, $C(a,b)$ is any contour in the complex plane from $a$ to $b$ that avoids
all poles of $\mathcal{Q}$. Now rewriting (3.48) for $y>b$ in the form

$\displaystyle\mathcal{V}(y)=-\left(\frac{b}{6}\right)^{1/2+\mathrm{i}p}\exp\left(\int_{b}^{y}\left(\frac{1}{2}+\mathrm{i}p\right)\frac{\mathrm{d}\eta}{\eta}\right)(1+o(1)),$
$\displaystyle
y\to+\infty,\qquad-\frac{1}{2}<\mathrm{Re}(\mathrm{i}q)<\frac{1}{2},$

we compare with (4.2) and eliminate $\mathcal{V}(y)$ to find

$\displaystyle\exp\left(\int_{-\infty}^{a}\left[\mathcal{Q}(\eta)-\mathrm{i}\sqrt{-\frac{\eta}{3}}+\left(\frac{1}{4}+\mathrm{i}\frac{p}{2}\right)\frac{1}{\eta}\right]\,\mathrm{d}\eta\right.$
$\displaystyle\left.\qquad\quad{}+\int_{C(a,b)}\mathcal{Q}(\eta)\,\mathrm{d}\eta+\int_{b}^{y}\left[\mathcal{Q}(\eta)-\left(\frac{1}{2}+\mathrm{i}p\right)\frac{1}{\eta}\right]\,\mathrm{d}\eta\right)$
$\displaystyle\qquad{}=\frac{\sqrt{2\pi}\mathrm{e}^{2\mathrm{i}(-a/3)^{3/2}}}{\tau
p\Gamma(\mathrm{i}p)}\mathrm{e}^{-\mathrm{i}\pi(1/2+\mathrm{i}p)/2}\left(b\sqrt{-\frac{a}{3}}\right)^{1/2+\mathrm{i}p}(1+o(1)),$
$\displaystyle\qquad\quad{}\
y\to+\infty,\quad-\frac{1}{2}\mathrm{Re}(\mathrm{i}q)<\frac{1}{2}.$ (4.3)

Taking the limit $y\to+\infty$ in (4.3) yields the total integral identity

$\displaystyle\exp\left(\int_{-\infty}^{a}\left[\mathcal{Q}(y)-\mathrm{i}\sqrt{-\frac{y}{3}}+\left(\frac{1}{4}+\mathrm{i}\frac{p}{2}\right)\frac{1}{y}\right]\,\mathrm{d}y\right.$
$\displaystyle\left.\qquad\quad{}+\int_{C(a,b)}\mathcal{Q}(y)\,\mathrm{d}y+\int_{b}^{+\infty}\left[\mathcal{Q}(y)-\left(\frac{1}{2}+\mathrm{i}p\right)\frac{1}{y}\right]\,\mathrm{d}y\right)$
$\displaystyle\qquad{}=\frac{\sqrt{2\pi}\mathrm{e}^{2\mathrm{i}(-a/3)^{3/2}}}{\tau
p\Gamma(\mathrm{i}p)}\mathrm{e}^{-\mathrm{i}\pi(1/2+\mathrm{i}p)/2}\left(b\sqrt{-\frac{a}{3}}\right)^{1/2+\mathrm{i}p},\qquad-\frac{1}{2}<\mathrm{Re}(\mathrm{i}q)<\frac{1}{2},$

in which the integral over the interval $(b,+\infty)$ is not generally
absolutely convergent, but it necessarily makes sense as an improper integral.
Suppose that we take for $C(a,b)$ a real path with infinitesimal semicircular
indentations in the lower half-plane centered at each real pole of
$\mathcal{Q}$. It is well-known that all poles of $\mathcal{Q}$ are simple and
have residue $\pm 1$. Therefore,

$\displaystyle\int_{C(a,b)}\mathcal{Q}(y)\,\mathrm{d}y=\mathrm{P.V.}\int_{a}^{b}\mathcal{Q}(y)\,\mathrm{d}y+\mathrm{i}\pi(N_{+}-N_{-}),$

where “P.V.” denotes the Hadamard principal value and $N_{\pm}$ is the number
of real poles of $\mathcal{Q}$ of residue $\pm 1$.

Setting $A:=-(\tfrac{2}{3})^{1/3}b<0$ and $B:=-(\tfrac{2}{3})^{1/3}a>0$ so all
real poles of
$u^{-}_{\mathrm{TT}}(x;\alpha)=-\big{(}\tfrac{3}{2}\big{)}^{1/3}\mathcal{Q}(-\big{(}\tfrac{3}{2}\big{)}^{1/3}x)$
lie in the interval $(A,B)$, we have established the total integral formula
(1.10) for the increasing tritronquée solution
$u=u^{-}_{\mathrm{TT}}(x;\alpha)$ of (1.1). This completes the proof of
Theorem 1.3.

##  5 Proof of Theorem 1.4

When $p>0$, which also implies that $\tau\in\mathbb{R}$, Riemann–Hilbert
Problem 2.1 has a unique solution for every $y\in\mathbb{R}$ as a consequence
of Zhou’s vanishing lemma [29]. Indeed, let us define a modified unknown
$\widetilde{\mathbf{W}}$ as $\widetilde{\mathbf{W}}:=\mathbf{W}$ for
$|\zeta|<1$ and
$\widetilde{\mathbf{W}}:=\mathbf{W}\zeta^{\mathrm{i}p\sigma_{3}}$ for
$|\zeta|>1$. This transformation removes the jump from the negative real axis
for $|\zeta|>1$, modifies the remaining jump matrices for $|\zeta|>1$ by
conjugation, and introduces a new diagonal jump across the unit circle. For
the latter, it is convenient in the setting of [29] to take the upper and
lower semicircles both to be oriented from left to right. The jump conditions
for $\widetilde{\mathbf{W}}(\zeta;y)$ are illustrated in Fig. 8.

Figure 8: The jump contour and jump matrix
$\widetilde{\mathbf{V}}^{\mathrm{PII}}(\zeta;y)$ for
$\widetilde{\mathbf{W}}(\zeta;y)$. On the segments of the jump contour within
the unit disk, the jump matrix is exactly as shown in Fig. 3.

Now $\widetilde{\mathbf{W}}(\zeta;y)$ satisfies the conditions of an identity-
normalized (as $\zeta\to\infty$) Riemann–Hilbert problem whose jump matrix is
easily checked to satisfy the cyclic condition at each self-intersection point
that is necessary for consistency. Moreover, when $p>0$ and $y\in\mathbb{R}$,
the jump matrix satisfies
$\widetilde{\mathbf{V}}^{\mathrm{PII}}(\zeta;y)=\widetilde{\mathbf{V}}^{\mathrm{PII}}(\zeta^{*};y)^{\dagger}$
for all non-real $\zeta$ and that
$\widetilde{\mathbf{V}}^{\mathrm{PII}}(\zeta;y)+\widetilde{\mathbf{V}}^{\mathrm{PII}}(\zeta;y)^{\dagger}$
is positive-definite for all real $\zeta$ in the jump contour. The vanishing
lemma [29] therefore guarantees the existence of
$\widetilde{\mathbf{W}}(\zeta;y)$, which also yields $\mathbf{W}(\zeta;y)$ by
inverting the substitution made for $|\zeta|>1$. Uniqueness of the solution is
standard, as is the fact that the solution satisfies
$\det(\mathbf{W}(\zeta;y))=1$.

The existence of a solution of Riemann–Hilbert Problem 2.1 for every real $y$
when $p>0$ actually implies via analytic Fredholm theory that the functions
$\mathcal{U}(y)$ and $\mathcal{V}(y)$ defined by (2.6) are analytic for
$y\in\mathbb{R}$. Moreover, these functions are Schwarz reflections of each
other: $\mathcal{V}(y)=-\mathcal{U}(y^{*})^{*}$. The relation (2.8) then
implies that $\mathcal{U}(y)$ and $\mathcal{V}(y)$ cannot have any real zeros
unless $p=0$ (in which case they vanish identically because
$\mathbf{W}(\zeta;y)\equiv\mathbb{I}$), because $\mathcal{U}(y)$ and
$\mathcal{V}(y)$ would have to vanish simultaneously. The same argument shows
that $\mathcal{U}(y)$ and $\mathcal{V}(y)$ cannot have any real critical
points for positive real $p$. Therefore, when $p>0$, $\mathcal{P}(y)$ and
$\mathcal{Q}(y)$ defined by (2.9) are nonvanishing analytic functions for
$y\in\mathbb{R}$, and they are related by the symmetry
$\mathcal{Q}(y)=\mathcal{P}(y^{*})^{*}$. This implies that when $p>0$, the
increasing tritronquée solution
$u=u^{-}_{\mathrm{TT}}\big{(}x;\tfrac{1}{2}+\mathrm{i}p\big{)}$ of (1.1) for
$\alpha=\tfrac{1}{2}+\mathrm{i}p$ is globally analytic and nonvanishing for
real $x$. Applying the symmetries (1.5) for $x\in\mathbb{R}$ then shows that
$u_{\mathrm{TT}}^{-}\big{(}x,-\tfrac{1}{2}+\mathrm{i}p\big{)}$,
$u_{\mathrm{TT}}^{+}\big{(}x,\tfrac{1}{2}-\mathrm{i}p\big{)}$, and
$u_{\mathrm{TT}}^{+}\big{(}x,-\tfrac{1}{2}-\mathrm{i}p\big{)}$ are global
nonvanishing solutions for $x\in\mathbb{R}$ whenever $p>0$. This completes the
proof of Theorem 1.4.

## Appendix A. Parabolic cylinder function parametrix

Consider the contour shown in Fig. 9 and the indicated jump matrix defined
thereon.

Figure 9: The jump contour in the $\lambda$-plane and jump matrix
$\mathbf{V}^{\mathrm{PC}}$.

######  Riemann–Hilbert Problem A.1 (parabolic cylinder parametrix).

Let $p,\tau\in\mathbb{C}$ be related by $\tau^{2}=\mathrm{e}^{2\pi p}-1$. Seek
a $2\times 2$ matrix-valued function
$\mathbf{P}(\lambda)=\mathbf{P}(\lambda;p,\tau)$ with the following
properties.

  * Analyticity: $\mathbf{P}(\lambda)$ is analytic for $\lambda$ in the five sectors $S_{0}\colon|\arg(\lambda)|<\tfrac{1}{4}\pi$, $S_{1}\colon\tfrac{1}{4}\pi<\arg(\lambda)<\tfrac{3}{4}\pi$, $S_{-1}\colon-\frac{3}{4}\pi<\arg(\lambda)<-\tfrac{1}{4}\pi$, $S_{2}\colon\tfrac{3}{4}\pi<\arg(\lambda)<\pi$, and $S_{-2}\colon-\pi<\arg(\lambda)<-\tfrac{3}{4}\pi$. It takes continuous boundary values on the excluded rays and at the origin from each sector.

  * Jump conditions: $\mathbf{P}_{+}(\lambda)=\mathbf{P}_{-}(\lambda)\mathbf{V}^{\mathrm{PC}}(\lambda)$, where $\mathbf{V}^{\mathrm{PC}}(\lambda)$ is the matrix function defined on the jump contour shown in Fig. 9.

  * Normalization: $\mathbf{P}(\lambda)\lambda^{\mathrm{i}p\sigma_{3}}\to\mathbb{I}$ as $\lambda\to\infty$ uniformly in all directions.

To solve Riemann–Hilbert Problem A.1, first note that there can be at most one
solution, and if it exists it must have unit determinant. This problem has no
solution if $p\in\mathrm{i}\mathbb{Z}\setminus\\{0\\}$ for the same reason as
in the case of Riemann–Hilbert Problem 2.1, and its solution is
$\mathbf{P}(\lambda)\equiv\mathbb{I}$ if $p=0$.

Therefore we now assume that $p\not\in\mathrm{i}\mathbb{Z}$, and consider the
matrix
$\mathbf{Q}(\lambda):=\mathbf{P}(\lambda)\mathrm{e}^{-\mathrm{i}\lambda^{2}\sigma_{3}/2}$.
This change of variables removes the exponentials
$\mathrm{e}^{\pm\mathrm{i}\lambda^{2}}$ from the jump conditions, from which
it follows that $\mathbf{Q}^{\prime}(\lambda)$ and $\mathbf{Q}(\lambda)$
satisfy exactly the same jump conditions. Appealing to invertibility and the
specified regular behavior of $\mathbf{P}$ near $\lambda=0$ we see that
$\mathbf{Q}^{\prime}(\lambda)\mathbf{Q}(\lambda)^{-1}$ is an entire function
of $\lambda$. If $\mathbf{P}(\lambda)$ has an asymptotic expansion as
$\lambda\to\infty$ consistent with the normalization condition and hence of
the form
$\mathbf{P}(\lambda)=\big{(}\mathbb{I}+\lambda^{-1}\mathbf{P}^{[1]}+\cdots\big{)}\lambda^{-\mathrm{i}p\sigma_{3}}$,
the additional assumption that the series is differentiable term-by-term with
respect to $\lambda$ shows by a Liouville argument that
$\mathbf{Q}^{\prime}(\lambda)\mathbf{Q}(\lambda)^{-1}=-\mathrm{i}\lambda\sigma_{3}+\mathrm{i}\big{[}\sigma_{3},\mathbf{P}^{[1]}\big{]}$;
therefore $\mathbf{Q}(\lambda)$ is a matrix solution of the first-order system

$\displaystyle\frac{\mathrm{d}\mathbf{Q}}{\mathrm{d}\lambda}=\begin{bmatrix}-\mathrm{i}\lambda&r\\\
s&\mathrm{i}\lambda\end{bmatrix}\mathbf{Q},\qquad
r:=2\mathrm{i}P^{[1]}_{12},\qquad s:=-2\mathrm{i}P^{[1]}_{21}.$ (A.1)

Eliminating the second row of $\mathbf{Q}$ and rescaling by
$t=\sqrt{2}\mathrm{e}^{-\mathrm{i}\pi/4}\lambda$ shows that the first row
elements of $\mathbf{Q}$ satisfy Weber’s equation for parabolic cylinder
functions [26, equation (12.2.2)]

$\displaystyle\frac{\mathrm{d}^{2}Q_{1k}}{\mathrm{d}t^{2}}-\left(\frac{1}{4}t^{2}+a\right)Q_{1k}=0,\qquad
a:=\frac{1}{2}(1+\mathrm{i}rs),\qquad k=1,2.$

This differential equation has particular solutions $U(a,\pm t)$ and
$U(-a,\pm\mathrm{i}t)$, and the special function $U(\cdot,\cdot)$ is described
in detail in [26, Section 12]. In each of the five sectors, the general
solution can be written as a linear combination of a suitable “numerically
satisfactory” fundamental pair:

$\displaystyle
Q_{1k}(\lambda)=\begin{cases}rA^{(0)}_{k}U(a,t)+rB^{(0)}_{k}U(-a,\mathrm{i}t),&\lambda\in
S_{0},\\\ rA^{(1)}_{k}U(a,t)+rB^{(1)}_{k}U(-a,-\mathrm{i}t),&\lambda\in
S_{1},\\\ rA^{(-1)}_{k}U(a,-t)+rB^{(-1)}_{k}U(-a,\mathrm{i}t),&\lambda\in
S_{-1},\\\ rA^{(2)}_{k}U(a,-t)+rB^{(2)}_{k}U(-a,-\mathrm{i}t),&\lambda\in
S_{2},\\\ rA^{(-2)}_{k}U(a,-t)+rB^{(-2)}_{k}U(-a,-\mathrm{i}t),&\lambda\in
S_{-2}.\end{cases}$

The first row of the matrix differential equation (A.1) then gives the
elements of the second row explicitly in terms of those of the first and their
derivatives; however the derivatives can be eliminated using [26, equations
(12.8.2)–(12.8.3)], and we therefore obtain

$\displaystyle
Q_{2k}(\lambda)=\sqrt{2}\mathrm{e}^{-\mathrm{i}\pi/4}\begin{cases}-A^{(0)}_{k}U(a-1,t)+\mathrm{i}\big{(}a-\tfrac{1}{2}\big{)}B^{(0)}_{k}U(1-a,\mathrm{i}t),&\lambda\in
S_{0},\\\
-A^{(1)}_{k}U(a-1,t)-\mathrm{i}\big{(}a-\tfrac{1}{2}\big{)}B^{(1)}_{k}U(1-a,-\mathrm{i}t),&\lambda\in
S_{1},\\\
A^{(-1)}_{k}U(a-1,-t)+\mathrm{i}\big{(}a-\tfrac{1}{2}\big{)}B^{(-1)}_{k}U(1-a,\mathrm{i}t),&\lambda\in
S_{-1},\\\
A^{(2)}_{k}U(a-1,-t)-\mathrm{i}\big{(}a-\tfrac{1}{2}\big{)}B^{(2)}_{k}U(1-a,-\mathrm{i}t),&\lambda\in
S_{2},\\\
A^{(-2)}_{k}U(a-1,-t)-\mathrm{i}\big{(}a-\tfrac{1}{2}\big{)}B^{(-2)}_{k}U(1-a,-\mathrm{i}t),&\lambda\in
S_{-2}.\end{cases}$

In each sector, the asymptotic expansion [26, equation (12.9.1)]:

$\displaystyle
U\left(a,t\right)\sim\mathrm{e}^{-\frac{1}{4}t^{2}}t^{-a-\frac{1}{2}}\sum_{j=0}^{\infty}(-1)^{j}\frac{{\big{(}\frac{1}{2}+a\big{)}_{2j}}}{j!(2t^{2})^{j}},\qquad
t\to\infty,\qquad|\arg(t)|<\frac{3}{4}\pi$ (A.2)

can be used to obtain the asymptotic behavior of the matrix elements of
$\mathbf{Q}(\lambda)$. Since the normalization condition on
$\mathbf{P}(\lambda)=\mathbf{Q}(\lambda)\mathrm{e}^{\mathrm{i}\lambda^{2}\sigma_{3}/2}=\mathbf{Q}(\lambda)\mathrm{e}^{-t^{2}\sigma_{3}/4}$
forbids exponential growth, it is necessary that $A_{1}^{(i)}=0$ and
$B_{2}^{(i)}=0$ for all sector indices $i=0,\pm 1,\pm 2$. From the jump
conditions it follows that the second column of $\mathbf{Q}$ must match
between sectors $S_{0}$ and $S_{1}$, as well as between sectors $S_{-1}$ and
$S_{-2}$. This implies that $A_{2}^{(0)}=A_{2}^{(1)}$ and that
$A_{2}^{(-1)}=A_{2}^{(-2)}$. Similarly, the first column of $\mathbf{Q}$ must
match between sectors $S_{0}$ and $S_{-1}$, as well as between sectors $S_{1}$
and $S_{2}$. This implies that $B_{1}^{(0)}=B_{1}^{(-1)}$ and that
$B_{1}^{(1)}=B_{1}^{(2)}$. Lastly, the jump condition between sectors $S_{2}$
and $S_{-2}$ implies that $B_{1}^{(-2)}=\mathrm{e}^{-2\pi
p}B_{1}^{(2)}=\mathrm{e}^{-2\pi p}B_{1}^{(1)}$ and that
$A_{2}^{(2)}=\mathrm{e}^{-2\pi p}A_{2}^{(-2)}=\mathrm{e}^{-2\pi
p}A_{2}^{(-1)}$. Therefore, if $\mathbf{Q}^{(i)}(\lambda)$ denotes the
restriction of $\mathbf{Q}(\lambda)$ to $\lambda\in S_{i}$, we have so far
found that

$\displaystyle\mathbf{Q}^{(0)}(\lambda)=\begin{bmatrix}rB_{1}^{(0)}U(-a,\mathrm{i}t)&rA_{2}^{(0)}U(a,t)\\\
\sqrt{2}\mathrm{e}^{\mathrm{i}\pi/4}\big{(}a-\tfrac{1}{2}\big{)}B_{1}^{(0)}U(1-a,\mathrm{i}t)&\sqrt{2}\mathrm{e}^{3\mathrm{i}\pi/4}A_{2}^{(0)}U(a-1,t)\end{bmatrix},\qquad\lambda\in
S_{0},$ (A.3)
$\displaystyle\mathbf{Q}^{(1)}(\lambda)=\begin{bmatrix}rB_{1}^{(1)}U(-a,-\mathrm{i}t)&rA_{2}^{(0)}U(a,t)\\\
\sqrt{2}\mathrm{e}^{-3\mathrm{i}\pi/4}\big{(}a-\tfrac{1}{2}\big{)}B_{1}^{(1)}U(1-a,-\mathrm{i}t)&\sqrt{2}\mathrm{e}^{3\mathrm{i}\pi/4}A_{2}^{(0)}U(a-1,t)\end{bmatrix},\qquad\lambda\in
S_{1},$
$\displaystyle\mathbf{Q}^{(-1)}(\lambda)=\begin{bmatrix}rB_{1}^{(0)}U(-a,\mathrm{i}t)&rA_{2}^{(-1)}U(a,-t)\\\
\sqrt{2}\mathrm{e}^{\mathrm{i}\pi/4}\big{(}a-\tfrac{1}{2}\big{)}B_{1}^{(0)}U(1-a,\mathrm{i}t)&\sqrt{2}\mathrm{e}^{-\mathrm{i}\pi/4}A_{2}^{(-1)}U(a-1,-t)\end{bmatrix},\qquad\lambda\in
S_{-1},$
$\displaystyle\mathbf{Q}^{(2)}(\lambda)=\begin{bmatrix}rB_{1}^{(1)}U(-a,-\mathrm{i}t)&\mathrm{e}^{-2\pi
p}rA_{2}^{(-1)}U(a,-t)\\\
\sqrt{2}\mathrm{e}^{-3\mathrm{i}\pi/4}\big{(}a-\tfrac{1}{2}\big{)}B_{1}^{(1)}U(1-a,-\mathrm{i}t)&\mathrm{e}^{-2\pi
p}\sqrt{2}\mathrm{e}^{-\mathrm{i}\pi/4}A_{2}^{(-1)}U(a-1,-t)\end{bmatrix},$
$\displaystyle\hphantom{\mathbf{Q}^{(2)}(\lambda)=}{}\ \lambda\in S_{2},$
$\displaystyle\mathbf{Q}^{(-2)}(\lambda)=\begin{bmatrix}\mathrm{e}^{-2\pi
p}rB_{1}^{(1)}U(-a,-\mathrm{i}t)&rA_{2}^{(-1)}U(a,-t)\\\ \mathrm{e}^{-2\pi
p}\sqrt{2}\mathrm{e}^{-3\mathrm{i}\pi/4}\big{(}a-\tfrac{1}{2}\big{)}B_{1}^{(1)}U(1-a,-\mathrm{i}t)&\sqrt{2}\mathrm{e}^{-\mathrm{i}\pi/4}A_{2}^{(-1)}U(a-1,-t)\end{bmatrix},$
$\displaystyle\hphantom{\mathbf{Q}^{(-2)}(\lambda)=}{}\ \lambda\in S_{-2}.$
(A.4)

Now taking $\lambda\to\infty$ in each sector we find that the matrix
$\mathbf{P}(\lambda)=\mathbf{Q}(\lambda)\mathrm{e}^{\mathrm{i}\lambda^{2}\sigma_{3}/2}$
satisfies the required normalization condition provided that

$\displaystyle a=\frac{1}{2}-\mathrm{i}p\qquad\text{which implies that}\quad
rs=-2p,$ (A.5)

and that the remaining unknown coefficients are subject to the following:

$\displaystyle rB_{1}^{(0)}=\mathrm{e}^{-\pi
p/4}\mathrm{e}^{\mathrm{i}p\ln(2)/2},\qquad$ $\displaystyle
A_{2}^{(0)}=\frac{\mathrm{e}^{-3\mathrm{i}\pi/4}}{\sqrt{2}}\mathrm{e}^{-\pi
p/4}\mathrm{e}^{-\mathrm{i}p\ln(2)/2},$ $\displaystyle
rB_{1}^{(1)}=\mathrm{e}^{3\pi p/4}\mathrm{e}^{\mathrm{i}p\ln(2)/2},\qquad$
$\displaystyle
A_{2}^{(-1)}=\frac{\mathrm{e}^{\mathrm{i}\pi/4}}{\sqrt{2}}\mathrm{e}^{3\pi
p/4}\mathrm{e}^{-\mathrm{i}p\ln(2)/2}.$ (A.6)

Only the values of the constants $r$ and $s$ remain undetermined. It is
interesting to observe that if one tries to use the definitions in (A.1) to
find $r$ and $s$ from the formulæ (A.3)–(A.4) after substituting from (A.5)
and (A.6) and using the asymptotic expansion (A.2) for the off-diagonal
entries, one finds simply the truisms $r=r$ and $s=s$. These constants must
therefore be determined from the only other information we have not yet used:
the nontrivial jumps captured by the first (resp. second) column of the jump
matrix $\mathbf{V}^{\mathrm{PC}}$ for
$\arg(\lambda)=\tfrac{1}{4}\pi,-\tfrac{3}{4}\pi$ (resp. for
$\arg(\lambda)=-\tfrac{1}{4}\pi,\tfrac{3}{4}\pi$). In particular, the jump
condition for $\arg(\lambda)=\tfrac{1}{4}\pi$ requires that

$\displaystyle Q^{(1)}_{11}(\lambda)=Q^{(0)}_{11}(\lambda)+\tau
Q^{(0)}_{12}(\lambda),$

which given all available information is equivalent to the condition

$\displaystyle\mathrm{e}^{\mathrm{i}\pi(a+1/2)/2}U(-a,-\mathrm{i}t)+\mathrm{e}^{-\mathrm{i}\pi(1+1/2)/2}U(-a,\mathrm{i}t)$
$\displaystyle\qquad{}=\tau
r\frac{\mathrm{e}^{-\mathrm{i}\pi/4}}{\sqrt{2}}\mathrm{e}^{-\pi
p/2}\mathrm{e}^{-\mathrm{i}p\ln(2)}U(a,t),\qquad
a=\frac{1}{2}-\mathrm{i}p,\qquad\arg(t)=0.$

But it is easy to check that this matches the connection formula [26, equation
(12.2.18)] for the parabolic cylinder function $U(\cdot,\cdot)$ provided that

$\displaystyle
r=r(p,\tau):=2\mathrm{e}^{\mathrm{i}\pi/4}\sqrt{\pi}\frac{\mathrm{e}^{\pi
p/2}\mathrm{e}^{\mathrm{i}p\ln(2)}}{\tau\Gamma(\mathrm{i}p)}.$ (A.7)

Using (A.5) one then deduces that

$\displaystyle
s=s(p,\tau)=-\frac{2p}{r(p,\tau)}=\frac{\mathrm{e}^{3\mathrm{i}\pi/4}}{\sqrt{\pi}}\tau
p\Gamma(\mathrm{i}p)\mathrm{e}^{-\pi p/2}\mathrm{e}^{-\mathrm{i}p\ln(2)}.$
(A.8)

At last all free parameters have been determined, and it is straightforward to
verify with the use of connection formulæ for $U(\cdot,\cdot)$ that the other
jump conditions that we have not yet enforced are automatically satisfied.

######  Remark A.2.

In fact, the unenforced jump conditions cannot contain any further
information, because the cyclic product of the jump matrices for $\mathbf{Q}$
is the identity (as must be the case for consistency at the origin). This
matrix identity modulo unit determinant amounts to three conditions, which
determine all of the off-diagonal entries of $\mathbf{V}^{\mathrm{PC}}$ given
any one of them.

This completes the solution of Riemann–Hilbert Problem A.1. We pause to point
out that the normalization condition on
$\mathbf{P}(\lambda)=\mathbf{P}(\lambda;p,\tau)$ indeed holds in every sense
we assumed at the start in order to derive the differential equation (A.1),
and therefore in particular from the asymptotic expansion (A.2),

$\displaystyle\mathbf{P}(\lambda;p,\tau)\lambda^{\mathrm{i}p\sigma_{3}}=\mathbb{I}+\frac{1}{2\mathrm{i}\lambda}\begin{bmatrix}0&r(p,\tau)\\\
-s(p,\tau)&0\end{bmatrix}+\begin{bmatrix}O\big{(}\lambda^{-2}\big{)}&O\big{(}\lambda^{-3}\big{)}\\\
O\big{(}\lambda^{-3}\big{)}&O\big{(}\lambda^{-2}\big{)}\end{bmatrix},\qquad\lambda\to\infty.$
(A.9)

In fact, (A.2) shows that the diagonal (resp. off-diagonal) elements have
explicit asymptotic expansions in descending even (resp. odd) powers of
$\lambda$. The full expansion of
$\mathbf{P}(\lambda)\lambda^{\mathrm{i}p\sigma_{3}}$ as $\lambda\to\infty$ is
the same in all five sectors.

In the special case that $p>0$, which also implies that $\tau\in\mathbb{R}$,
we may use the identity for the modulus of the gamma function on the imaginary
axis (see [26, equation (5.4.3)]) to deduce that $s=-r^{*}$.

### Acknowledgements

The author’s work was supported by the National Science Foundation under
grants DMS-1513054 and DMS-1812625. The author thanks Thomas Bothner, Deniz
Bilman, and Liming Ling for useful discussions.

## References

  * [1] Baik J., Buckingham R., DiFranco J., Its A., Total integrals of global solutions to Painlevé II, Nonlinearity 22 (2009), 1021–1061, arXiv:0810.2586.
  * [2] Bertola M., Tovbis A., Universality for the focusing nonlinear Schrödinger equation at the gradient catastrophe point: rational breathers and poles of the tritronquée solution to Painlevé I, Comm. Pure Appl. Math. 66 (2013), 678–752, arXiv:1004.1828.
  * [3] Bilman D., Ling L., Miller P.D., Extreme superposition: rogue waves of infinite order and the Painlevé-III hierarchy, arXiv:1806.00545.
  * [4] Bothner T., Transition asymptotics for the Painlevé II transcendent, Duke Math. J. 166 (2017), 205–324, arXiv:1502.03402.
  * [5] Bothner T., Miller P.D., Rational solutions of the Painlevé-III equation: large parameter asymptotics, arXiv:1808.01421.
  * [6] Bothner T., Miller P.D., Sheng Y., Rational solutions of the Painlevé-III equation, Stud. Appl. Math. 141 (2018), 626–679, arXiv:1801.04360.
  * [7] Boutroux P., Recherches sur les transcendantes de M. Painlevé et l’étude asymptotique des équations différentielles du second ordre, Ann. Sci. École Norm. Sup. 30 (1913), 255–375.
  * [8] Boutroux P., Recherches sur les transcendantes de M. Painlevé et l’étude asymptotique des équations différentielles du second ordre (suite), Ann. Sci. École Norm. Sup. 31 (1914), 99–159.
  * [9] Buckingham R.J., Miller P.D., The sine-Gordon equation in the semiclassical limit: critical behavior near a separatrix, J. Anal. Math. 118 (2012), 397–492, Corrigenda, J. Anal. Math. 119 (2013), 403–405, arXiv:1106.5716.
  * [10] Buckingham R.J., Miller P.D., Large-degree asymptotics of rational Painlevé-II functions: noncritical behaviour, Nonlinearity 27 (2014), 2489–2578, arXiv:1310.2276.
  * [11] Buckingham R.J., Miller P.D., Large-degree asymptotics of rational Painlevé-II functions: critical behaviour, Nonlinearity 28 (2015), 1539–1596, arXiv:1406.0826.
  * [12] Claeys T., Grava T., Painlevé II asymptotics near the leading edge of the oscillatory zone for the Korteweg–de Vries equation in the small-dispersion limit, Comm. Pure Appl. Math. 63 (2010), 203–232, arXiv:0812.4142.
  * [13] Clarkson P.A., On Airy solutions of the second Painlevé equation, Stud. Appl. Math. 137 (2016), 93–109, arXiv:1510.08326.
  * [14] Deift P.A., Zhou X., Asymptotics for the Painlevé II equation, Comm. Pure Appl. Math. 48 (1995), 277–337.
  * [15] Dubrovin B., Grava T., Klein C., On universality of critical behavior in the focusing nonlinear Schrödinger equation, elliptic umbilic catastrophe and the tritronquée solution to the Painlevé-I equation, J. Nonlinear Sci. 19 (2009), 57–94, arXiv:0704.0501.
  * [16] Flaschka H., Newell A.C., Monodromy- and spectrum-preserving deformations. I, Comm. Math. Phys. 76 (1980), 65–116.
  * [17] Fokas A.S., Its A.R., Kapaev A.A., Novokshenov V.Yu., Painlevé transcendents: the Riemann–Hilbert approach, Mathematical Surveys and Monographs, Vol. 128, Amer. Math. Soc., Providence, RI, 2006.
  * [18] Gambier B., Sur les équations différentielles du second ordre et du premier degré dont l’intégrale générale est a points critiques fixes, Acta Math. 33 (1910), 1–55.
  * [19] Jimbo M., Miwa T., Monodromy preserving deformation of linear ordinary differential equations with rational coefficients. II, Phys. D 2 (1981), 407–448.
  * [20] Joshi N., Mazzocco M., Existence and uniqueness of tri-tronquée solutions of the second Painlevé hierarchy, Nonlinearity 16 (2003), 427–439, math.CA/0212117.
  * [21] Kapaev A.A., Asymptotic expressions for the second Painlevé functions, Theoret. and Math. Phys. 77 (1988), 1227–1234.
  * [22] Kapaev A.A., Global asymptotics of the second Painlevé transcendent, Phys. Lett. A 167 (1992), 356–362.
  * [23] Kapaev A.A., Scaling limits in the second Painlevé transcendent, J. Math. Sci. 83 (1997), 38–61.
  * [24] Lu B.-Y., Miller P.D., Universality near the gradient catastrophe point in the semiclassical sine-Gordon equation, in preparation.
  * [25] Novokshenov V.Yu., Tronquée solutions of the Painlevé II equation, Theoret. and Math. Phys. 172 (2012), 1136–1146.
  * [26] Olver F.W.J., Olde Daalhuis A.B., Lozier D.W., Schneider B.I., Boisvert R.F., Clark C.W., Miller B.R., Saunders B.V. (Editors), NIST digital library of mathematical functions, Release 1.0.17, 2017, available at http://dlmf.nist.gov/.
  * [27] Olver S., RHPackage, Version 0.4, 2011, available at https://github.com/dlfivefifty/RHPackage.
  * [28] Tracy C.A., Widom H., Level-spacing distributions and the Airy kernel, Comm. Math. Phys. 159 (1994), 151–174, hep-th/9211141.
  * [29] Zhou X., The Riemann–Hilbert problem and inverse scattering, SIAM J. Math. Anal. 20 (1989), 966–986.


††thanks: These authors contributed equally to this study.††thanks: These
authors contributed equally to this study.††thanks: These authors contributed
equally to this study.

# Chern numbers of topological phonon band crossing determined with inelastic
neutron scattering

Zhendong Jin  Biaoyan Hu  Yiran Liu  International Center for Quantum
Materials, School of Physics, Peking University, Beijing 100871, China  Yangmu
Li  Beijing National Laboratory for Condensed Matter Physics, and Institute of
Physics, Chinese Academy of Sciences, Beijing 100190, China  Condensed Matter
Physics and Materials Science Division, Brookhaven National Laboratory, Upton,
New York 11973, USA  University of Chinese Academy of Sciences, Beijing
100049, China  Tiantian Zhang  Beijing National Laboratory for Condensed
Matter Physics, and Institute of Physics, Chinese Academy of Sciences, Beijing
100190, China  Kazuki Iida  Kazuya Kamazawa  Neutron Science and Technology
Centre, Comprehensive Research Organisation for Science and Society (CROSS),
Tokai, Ibaraki 319-1106, Japan  A. I. Kolesnikov  M. B. Stone  Neutron
Scattering Division, Oak Ridge National Laboratory, Oak Ridge, Tennessee
37831, USA  Xiangyu Zhang  Haiyang Chen  Yandong Wang  State Key Laboratory
for Advance Metals and Materials, University of Science and Technology
Beijing, Beijing 10083, China  I. A. Zaliznyak  J. M. Tranquada  Condensed
Matter Physics and Materials Science Division, Brookhaven National Laboratory,
Upton, New York 11973, USA  Chen Fang<EMAIL_ADDRESS>Beijing National
Laboratory for Condensed Matter Physics, and Institute of Physics, Chinese
Academy of Sciences, Beijing 100190, China  Yuan Li<EMAIL_ADDRESS>International Center for Quantum Materials, School of Physics, Peking
University, Beijing 100871, China  Collaborative Innovation Center of Quantum
Matter, Beijing 100871, China

###### Abstract

Topological invariants in the band structure, such as Chern numbers, are
crucial for the classification of topological matters and dictate the
occurrence of exotic properties, yet their direct spectroscopic determination
has been largely limited to electronic bands. Here, we use inelastic neutron
scattering in conjunction with ab initio calculations to identify a variety of
topological phonon band crossings in MnSi and CoSi single crystals. We find a
distinct relation between the Chern numbers of a band-crossing node and the
scattering intensity modulation in momentum space around the node. Given
sufficiently high resolution, our method can be used to determine arbitrarily
large Chern numbers of topological phonon band-crossing nodes.

##  I Introduction

Ever since the discovery of topological quantum numbers in quantum Hall states
[1, 2], the concept of band topology has shed light on the exploration and
classification of crystalline materials [3, 4, 5, 6, 7]. Topological
insulators, semimetals, and superconductors are extensively studied, both as
novel phases of matter and for their potential applications. Unlike
conventional phases of matter described by symmetry in the Landau paradigm,
topological phases are classified by topological invariants, which do not
change over adiabatic deformations of the band structure.

An important topological invariant is called the Chern number, which is
associated with a mapping from a two-dimensional (2D) closed surface in
reciprocal space to the Hilbert space of Bloch states. The Chern number
characterizes the topological structure of such mapping and has observable
consequences. In the gapped energy spectrum of 2D quantum Hall systems,
nonzero Chern numbers correspond to the number of edge states which lead to
the quantization of the Hall conductance [8, 3]. In three-dimensional (3D)
Weyl semimetals, Weyl nodes act as monopoles of Berry flux and have nonzero
Chern numbers (defined by the mapping from their enclosing surface in momentum
space to the Hilbert space), which dictate the number of Fermi-arc surface
states [9, 10, 11] and the quantized magnitude of circular photogalvanic
effect [12, 13]. As band topology is independent of the statistics of the
constituent quasiparticles, similar phenomena are also expected in bosonic
systems. For example, topological photonic and acoustic bands and their
corresponding surface states have been found in artificial structures [14, 15,
16, 17, 18]. In natural crystals, a variety of Dirac and Weyl nodes have been
predicted and/or observed in phonon [19, 20, 21, 22, 23, 24, 25] and magnon
bands [26, 27, 28, 29, 30]. Previous experiments on these systems mainly
focused on the bulk dispersion relation near the band-crossing points, rather
than on topological invariants such as the Chern numbers, partly because it is
difficult to measure phonon and magnon surface states and transport behaviors.

Figure 1:  (a) Cubic primitive cell of $M$Si ($M$ = Mn, Co). (b) and (c)
Representative INS intensities of MnSi and CoSi, respectively, along a high-
symmetry momentum trajectory, measured at $T$ = 40 K. Data measured with
$E_{\rm i}=57$ meV and 90 meV are combined after proper intensity
normalization. (d) The Brillouin zone (BZ), with the momentum trajectory
marked in red. (e) and (f) $\cal S(\mathbf{Q},\omega)$ calculated from the
fitted-force-constant model along the same trajectories as in (b) and (c).
Solid lines in (b), (c), (e), and (f) indicate the calculated phonon
dispersions.

It is experimentally possible to determine Chern numbers in phonon and magnon
bands, if one can measure the topological structure of wave functions, i.e.,
eigenvectors of motion, in momentum space. Scattering methods such as
inelastic neutron scattering (INS) and X-ray scattering are suitable for this
purpose, because their dynamical structure factor $\cal S(\mathbf{Q},\omega)$
is related to the excitations’ eigenvectors [31, 32]. As the eigenvectors vary
strongly in the vicinity of Weyl and Dirac points, the observed intensities
are expected to undergo strong and distinct modulations, which can reflect the
topological structure. Such modulations have been recently reported in several
topological magnon systems around Dirac points [33, 34] and nodal lines [35,
36].

Here, we report an INS study of $M$Si ($M=$ Co, Mn) single crystals, which
host multiple types of topological phonon band crossing nodes [21, 37, 22]. By
comparing the observed INS intensities with our fitted model based on density
functional perturbation theory (DFPT) calculations, we verify the
theoretically predicted coexistence of two-fold quadruple Weyl points, three-
fold spin-1 Weyl points, and four-fold charge-2 Dirac points in $M$Si. We
further explore the spectroscopic signatures of topological structures near
the topological band crossing points and show, theoretically and in some cases
with comparison to the experimental data, that the number of intensity extrema
on a momentum sphere enclosing the band-crossing node equals the Chern number
of the node. Our result demonstrates the capability of INS for direct Chern-
number determination.

This paper is organized as follows: In Section II, we describe the INS
experiment and the fitted-force-constant model based on DFPT. In Sec. III, we
show phonon dispersions, both in a global view and close to topological band
crossing nodes. In Sec. IV, we discuss INS spectroscopic features near the
band-crossing nodes and investigate their relation with the Chern numbers. In
Sec. V, we make a brief discussion and a summary.

##  II Experiment and calculation methods

###  II.1 INS experiment

High-quality single crystals of MnSi and CoSi were grown by a traveling-
solvent floating zone method. The INS experiments were performed on the
4SEASONS spectrometer at MLF, J-PARC, Japan and the SEQUOIA spectrometer at
SNS, ORNL, USA [38, 39]. A total of 33 (28) grams of MnSi (CoSi) twin-free
single crystals with a mosaic spread of $\lesssim 1.3^{\circ}$ full-width at
half-maximum (FWHM) were used for the experiments (Fig. S5 [40, 41]). The INS
data shown in this paper were collected with incident neutron energies $E_{\rm
i}$ = 57 and 90 meV at a fixed temperature of $T$ = 40 K, and analyzed with
the utsusemi and horace software [42, 43]. As MnSi and CoSi share the
$B$20-type structure belonging to the non-centrosymmetric space group
$P2_{1}3$ [Fig.1(a)], data from equivalent momenta have been symmetrized and
averaged accordingly in order to improve counting statistics. Intensities are
presented in absolute scattering cross sections by using the incoherent
elastic scattering of the sample for normalization [44]. To best visualize the
phonon cross sections, we present coherent scattering signals from the sample
only, whereas intensities arising from incoherent scattering of the sample and
from the aluminum sample holders have been subtracted as background (Fig. S6
[40]).

Figure 2:  (a) INS spectra near a spin-1 Weyl point in CoSi, plotted along a
R-$\Gamma$-R trajectory. Colored dashed lines correspond to energy cuts in
(b), which are fitted with a sum of Gaussian profiles on a linear background,
by assuming a total of one and three peak(s) at and away from the
$\Gamma$-point, respectively, and under the constraint that equivalent
$\mathbf{q}$ positions must have the same energies and peak widths. (c) and
(d) Phonon intensities from INS experiment and the fitted-force-constant
model, respectively, plotted along a X-R-X trajectory near a charge-2 Dirac
point in CoSi. (e) and (f) Similar to (c) and (d), but for MnSi and along an
M-R-M trajectory. (g) Energy cuts at a series of $\mathbf{Q}$ points using the
same data as in (e) and (f). The INS and fitted model calculated data are
displayed by circles and dashed lines, respectively. Open squares denote peak
positions estimated from two-peak fits to the data (the fits are not shown),
which form an approximate linear band crossing. Data in (b) and (g) are offset
for clarity.

###  II.2 Fitted-force-constant model

The phonon force constant matrices of CoSi and MnSi were calculated with the
Vienna ab intio simulation package (vasp) [45, 46, 47, 48, 49] using the DFPT
method. The calculations were done with the Perdew-Burke-Ernzerhof (PBE) type
exchange-correlation functional, under the generalized gradient approximation
(GGA) [50]. The kinetic energy cutoff was set to 400 eV. Integrations over the
Brillouin zone were performed with Monkhorst-Pack $\mathbf{Q}$-point grids
(equivalent to $12\times 12\times 12$ grid for CoSi and $15\times 15\times 15$
grid for MnSi). Lattice constants and atomic positions were relaxed until
residual forces drop below 0.001 eV/Å. The relaxed lattice constants were 4.35
Å for CoSi and 4.42 Å for MnSi, which are slightly smaller than our
experimental values, $4.43\pm 0.02$ Å for CoSi and $4.56\pm 0.01$ Å for MnSi.
After obtaining the band dispersion $\omega(\mathbf{Q})$, the coherent
dynamical structure factors ${\cal S}_{\text{coh}}^{(i)}(\mathbf{Q},\omega)$
for all vibration modes $(i)$ were written as [31, 32, 51]

$\displaystyle\begin{split}{\cal
S}_{\text{coh}}^{(i)}(\mathbf{Q},\omega)&=\frac{(2\pi)^{3}}{v_{0}}\sum_{\mathbf{G},\mathbf{q}}\frac{|F^{(i)}_{\text{coh}}(\mathbf{Q})|^{2}}{2\omega^{(i)}(\mathbf{q})}\delta(\mathbf{Q}-\mathbf{q}-\mathbf{G})\delta(\omega-\omega^{(i)}(\mathbf{q})),\\\
F^{(i)}_{\rm coh}({\bf
Q})&=\sum_{d}\frac{b_{d,\text{coh}}}{\sqrt{m_{d}}}\mathbf{Q}\cdot\bm{\xi}^{(i)}_{d}(\mathbf{q})e^{i\mathbf{Q}\cdot\mathbf{r}_{d}},\end{split}$
(1)

where $\mathbf{Q}=\mathbf{q}+\mathbf{G}$ is the total momentum transfer with
$\mathbf{G}$ being a reciprocal lattice vector, and $m_{d}$, $\mathbf{r}_{d}$,
$b_{d,\text{coh}}$, $\bm{\xi}^{s}_{d}(\mathbf{q})$ denote the mass, position,
coherent scattering length and eigen vector of the $d$th atom in the primitive
cell.

Notably, the original DFPT results have deviations from the INS data,
including a global rescaling in energies, slight distortions in the
dispersions and scattering intensity (Fig. S9 [40]). Such deviations may stem
from inaccuracy in the calculated crystal structure and force constants. To
improve the accuracy of our model, we have performed parametric fits on the
leading force constants while discarding the weaker interactions. The energy
values of all phonon branches at high-symmetry points ($\Gamma$, X, M, and R)
are extracted from the experimental spectra (Table S1 [40]) and used in the
fitting. By adopting about ten pairwise interactions, we are able to reproduce
features in the experimental dispersions and intensities with satisfactory
accuracy over many BZs (Fig. S9 [40]). All model calculation results in the
main text are obtained with the optimized parameters listed in Table S2 [40].

##  III Dispersion and topological band crossing

To begin, we present in Fig. 1(b) and (c) representative INS spectra of MnSi
and CoSi, respectively, along high-symmetry lines in the irreducible Brillouin
zone (BZ) [Fig. 1(d)]. A wealth of phonon scattering signals are observed.
Strong contaminations below 10 meV are due to multiple scattering, and
residual aluminum phonon scattering are observed below 25 meV (especially in
the CoSi data). Overall, the phonon INS signals compare very well with our
model calculations in Fig. 1(e) and (f). The model for MnSi is particularly
satisfactory. This suggests that our fitted-force-constant model provides a
good approximate representation of the phonons.

Having established our fitted model, we now zoom into the topological band
crossings. We start from the $\Gamma$-point at the BZ center. According to
group theoretical analysis, phonons at the $\Gamma$-point are irreducible
representations of the tetrahedral $T(23)$ point group: $\Gamma=2A+2E+6T$,
where $A$, $E$, and $T$ represent singly, doubly, and triply degenerated
modes, respectively. Importantly, all triply degenerate phonons at the
$\Gamma$-point are spin-1 Weyl points protected by the point-group symmetry
[21], and all doubly degenerate phonons are quadruple Weyl points protected by
the extra time-reversal symmetry $\mathcal{T}$ [37, 52]. We have additionally
verified the irreducible representations of the BZ-center phonons using
polarized Raman spectroscopy (Fig. S7 [40]).

Based on the above information, we zoom into the spin-1 Weyl point at the
highest energy (52.5 meV) in CoSi. The INS intensities along a R-$\Gamma$-R
momentum trajectory are displayed in Fig. 2(a). This band crossing has
relatively large dispersion velocities and is far away from other bands, yet
still, the fact that a total of three branches are involved in the crossing
makes them challenging to resolve experimentally. By making energy line cuts
and fitting the intensity profiles systematically [Fig. 2(b)], we find that
the INS data do support a crossing of three bands. For comparison, fitting the
spectra obtained symmetrically away from $\Gamma$ with only two or fewer peaks
does not yield a consistent description (Fig. S8(a) [40]). We thus conclude
that the spin-1 Weyl point in CoSi is experimentally confirmed. In MnSi,
however, all three-fold degenerate modes at the $\Gamma$-point are too close
in energy to other phonon branches [Fig.1(b) and (e)], precluding a similar
confirmation to be made.

We next turn to the R-point at the BZ corner, where all band crossings are
four-fold degenerate “charge-2 Dirac points” ensured by the crystallographic
and $\mathcal{T}$ symmetries [21]. A charge-2 Dirac point is the direct sum of
two identical spin-1/2 Weyl points. Because the bands remain two-fold
degenerate along the $\overline{\rm RX}$ and $\overline{\rm RM}$ directions,
we expect to observe only two linearly dispersing branches along these
directions, which can be regarded as a key signature of the Dirac points. INS
spectra consistent with such understanding are displayed near 41 meV for MnSi
and 52 meV for CoSi in Fig. 2(c) and (e), respectively. The band crossings are
also reproduced in our model calculations [Fig. 2(d) and (f)]. Energy cuts at
a series of successive $\mathbf{Q}$ positions [Fig. 2(e)] further confirm the
approximate linear band crossing [Fig. 2(g)], and a similar case for CoSi is
displayed in Fig. S8(b) [40]. For MnSi, the fitted-force-constant model is
furthermore able to account for the scattering cross sections quite
accurately, as seen from the colored dashed lines in Fig. 2(g), which actually
represent model calculated intensities rather than peak fitting. This
quantitative agreement suggests that it is possible to use the intensity
information, from the experiment and/or the calculation, to elucidate the
Chern number of a topological band crossing, which is our next subject.

##  IV Detection of Chern numbers by neutron scattering

###  IV.1 General theoretical scheme for two-fold Weyl points

In this section, we will derive the explicit relation among phonon
eigenvectors, Chern numbers, and the INS dynamical structure factor, using
low-energy effective models near topological band crossing points. We will
first present theoretical considerations for two-fold Weyl points. Then, we
will use the general formulism in our specific analyses of four-fold charge-2
Dirac points ($2\times$ two-fold Weyl points) and quadruple Weyl points.

To begin with, for a two-fold Weyl point, the effective Hamiltonian can be
written as a $2\times 2$ Hermitian matrix

$\displaystyle H_{2\times
2}(\mathbf{q})=\sum_{i=x,y,z}{f_{i}}(\mathbf{q})\cdot{\sigma_{i}}+f_{0}(\mathbf{q})\sigma_{0},$
(2)

where $\mathbf{q}$ is momentum measured from the Weyl point and ${\sigma_{i}}$
are the Pauli matrices. In this notation, the eigenvector of one of the bands
(e.g., the upper band) is represented by a spinor
$\bm{\xi}^{\text{up}}(\mathbf{q})=(\psi_{1}(\mathbf{q}),\psi_{2}(\mathbf{q}))^{\text{T}}$.
We can further define a pseudospin quantity $\mathbf{S}({\mathbf{q}})$

$\displaystyle
S_{i}(\mathbf{q})\equiv\frac{f_{i}(\mathbf{q})}{|\mathbf{f}(\mathbf{q})|}=\left\langle\xi^{\text{up}}\left|\sigma_{i}\right|\xi^{\text{up}}\right\rangle,$
(3)

whose direction is represented by a point on the Bloch sphere
$\mathbb{S}_{\mathrm{B}}$. If we consider a surface $\mathbb{S}_{\mathbf{q}}$
enclosing the Weyl point in momentum space, on which a gap always exists
between the upper and lower bands, a wrapping number can be used to
characterize the mapping from $\mathbb{S}_{\mathbf{q}}$ to
$\mathbb{S}_{\mathrm{B}}$. This number is the Chern number $\pm C$ of the Weyl
node [37], which does not depend on the shape of $\mathbb{S}_{\mathbf{q}}$.
Specifically, $\mathbf{S}({\mathbf{q}})$ will take every possible directions
on $\mathbb{S}_{\mathrm{B}}$ at least $|C|$ times as $\mathbf{q}$ moves around
$\mathbb{S}_{\mathbf{q}}$. We will see examples of this in Figs. 3(b) and
4(a).

Next, we show that the pseudospin texture $\mathbf{S}({\mathbf{q}})$ can leave
distinct signatures in the INS dynamical structure factors. As we can see from
Eq. (1), the dynamic structure factor ${\cal
S}_{\text{coh}}^{(i)}(\mathbf{Q},\omega)$ is sensitive to the inner product of
$\bm{\xi}^{(i)}_{d}(\mathbf{q})$ and $\mathbf{Q}$. Taking acoustic phonons for
example, the INS intensity vanishes if the polarization vector $\bm{\xi}$
(same for all atoms for acoustic phonons) lie perpendicular to $\mathbf{Q}$,
and reaches maximum when the two vectors are parallel. For optical branches,
similar conclusions do exist by generalizing the real-space polarization
vectors to abstract phonon eigenvectors in the Hilbert space. The main
difference is that the eigenvector now has $3N$ distinct components (for
describing the collective vibration of $N$ atoms in the primitive cell) which
are generically complex numbers. In the close vicinity to a Weyl node,
assuming that $|\mathbf{Q}|\gg|\mathbf{q}|$ and $|\mathbf{q}|\ll 2\pi/a$, so
that $\mathbf{Q}\approx\mathbf{G}$, the only fast varying term in the formula
is the phonon eigenvector $\bm{\xi}^{\text{up}}(\mathbf{q})$. Using the
pseudospin quantities, it is straightforward to show that the INS intensity of
the upper band is approximately [40]

$\displaystyle{\cal S}_{\text{coh}}^{\rm up}(\mathbf{Q},\omega)$
$\displaystyle\propto|\mathbf{V}|\left(1+\frac{\mathbf{S}(\mathbf{q})\cdot\mathbf{V}(\mathbf{G})}{\left|\mathbf{S}(\mathbf{q})\right|\left|\mathbf{V}(\mathbf{G})\right|}\right),$
(4)

where $\mathbf{V}(\mathbf{G})$ is a constant vector that does not sensitively
depend on $\mathbf{q}$ but varies between BZs. Both the pseudospin
$\mathbf{S}$ and the vector $\mathbf{V}$ are vectors in the effective two-band
Hilbert space, which ultimately encode the phonon eigenvectors and their
sampling by INS in the BZ of $\mathbf{G}$. It then becomes clear that the INS
intensity is related to a projection of $\mathbf{S}(\mathbf{q})$ along
$\mathbf{V}$. In particular, ${\cal S}_{\text{coh}}^{\rm
up}(\mathbf{Q},\omega)$ would reach maximum when $\mathbf{S}(\mathbf{q})$ and
$\mathbf{V}$ are parallel, and vanish when they are anti-parallel.

Combining Eq. (4) with our former arguments about the wrapping (Chern) number,
we come to the following explicit statement: On a momentum surface
$\mathbb{S}_{\mathbf{q}}$ that encloses the Weyl node, there are at least
$|C|$ momenta where pseudospin $\mathbf{S}$ is parallel to $\mathbf{V}$, where
the INS intensity of the upper band ${\cal S}_{\text{coh}}^{\rm
up}(\mathbf{Q},\omega)$ reaches maximum. Similarly, there are also at least
$|C|$ momenta where the intensity reaches zero. In most practical cases where
the lowest-order $k\cdot p$ theory holds, the total numbers of maxima and
zeros are simply $|C|$. The approximation of
$\mathbf{Q}=\mathbf{q}+\mathbf{G}\approx\mathbf{G}$ can also be released since
$\mathbf{Q}$ has trivial topology on the surface $\mathbb{S}_{\mathbf{q}}$ as
long as the origin of $\mathbf{Q}$ is not enclosed. In summary, a Weyl node
serves as a singular point of the pseudospin in momentum space, and the
contrasting INS intensity distribution around it reveals its Chern number.

###  IV.2 The charge-2 Dirac point

Figure 3:  (a) Phonon dispersions near the charge-2 Dirac point at 29.6 meV in
MnSi. Solid color indicates energy integration range (30–32 meV). (b)
Pseudospin texture on a $\mathbf{q}$-sphere enclosing the R-point. Color
indicates energy-integrated INS intensities near
$\mathbf{Q}=(2.5,\,2.5,\,2.5)$. (c) INS intensity along the red circle
$\overset{\frown}{ABA}$ in (b), fitted with a cosine function (red line, see
text). (d) and (e) INS intensities averaged over solid $\mathbf{q}$-spheres
enclosing $(2.5,\,2.5,\,2.5)$ and $(2.5,\,3.5,\,3.5)$ respectively, with a
radius of 0.2 r.l.u. and a cone-smoothing width of $\pm 20^{\circ}$. (e) and
(f) Fitted-force-constant model calculations according to Eq. (4), for the
same $\mathbf{q}$-spheres as in (d) and (e), respectively.

Now that we have linked the pseudospin texture with the INS intensity
distribution, we next use it to analyze the charge-2 Dirac point in the phonon
bands of $M$Si at the R-point of the BZ. Although being four-fold degenerate,
the effective Hamiltonian near the band crossing
$H_{4}(\mathbf{q})\propto\left(\begin{smallmatrix}\mathbf{q}\cdot\bm{\sigma}&0\\\
0&\mathbf{q}\cdot\bm{\sigma}\end{smallmatrix}\right)$ is the direct sum of
that of two identical spin-1/2 Weyl points, each with Chern number $C=\pm
1$[21]. The general statement between pseudospin and INS intensity can be
easily generalized in this case with only minor modifications [40]. Namely, as
long as the INS intensities of the two upper bands are considered as a whole,
there will be exactly one maximum and one minimum on the enclosing
$\mathbb{S}_{\mathbf{q}}$, i.e., resembling that of a regular spin-1/2 Weyl
point. The only difference is an extra constant term in the intensity, so that
the minimum is finite rather than zero [40].

For a concrete example, we inspect the INS intensities near a charge-2 Dirac
point (at the BZ R-point) at about 29.6 meV in MnSi. The phonon dispersion
nearby is schematically shown in Fig. 3(a). On a small $\mathbf{q}$-sphere
around the R-point, the pseudospin texture is visualized in Fig. 3(b), where
purple arrows indicate the pseudospin directions $\mathbf{S}(\mathbf{q})$. The
outward hedgehog configuration of the arrows indicates that
$\mathbf{S}(\mathbf{q})$ takes every direction once on the Bloch sphere, i.e.,
the Chern number is $|C|=1$.

From an experimental perspective, a $\mathbf{q}$-sphere as small as possible
should be used to extract the INS intensity, in order to avoid an overlap with
neighboring bands. But the sphere cannot be too small or thin, as otherwise
the counting statistics would be too low. In our case, we find that binning
INS data over a finite solid sphere produces satisfactory results. In Fig.
3(d), intensities integrated over [30, 32] meV around R = (2.5, 2.5, 2.5) are
displayed as a function of the polar and azimuthal angles. The radius of the
solid sphere is set to be 0.2 reciprocal lattice units (r.l.u.), and the
displayed intensity at each angle represents the average over a cone volume
within a $20^{\circ}$ half-apex angle.

The above INS result agrees nicely with our fitted-force-constant model
calculations [Fig. 3(f)]. To rationalize their characteristics with our
effective model in Eq. (4), we note that due to the three-fold rotational
symmetry, the vector ${\mathbf{V}}$ for the Dirac point at R = (2.5, 2.5, 2.5)
is along the $[111]$ direction. Consequently, the intensity will reach its
maximum and minimum along the $[111]$ and $[\bar{1}\bar{1}\bar{1}]$ directions
(or vice versa). Moreover, on a great circle $\overset{\frown}{ABA}$ [red
dashed lines in Fig. 3(b) and (d)] that passes through the $[111]$ and
$[\bar{1}\bar{1}\bar{1}]$ directions, the intensity is expected to have a
cosine behavior, which we also confirm experimentally [Fig. 3(c)].

Besides for the R-point (2.5, 2.5, 2.5) where the symmetry is high, the same
analysis can be performed at other R-points, such as (2.5, 3.5, 3.5). The
results [Fig. 3(e) and (g)] show that the extrema are no longer located along
the diagonal direction, due to the different orientation of $\mathbf{V}$.
Nevertheless, the fact that the intensity exhibits one minimum and one maximum
indicates that the underlying Chern number is $|C|=1$.

###  IV.3 The quadruple Weyl point

We further explore the quadruple Weyl points at the $\Gamma$-point, which are
protected by time-reversal symmetry $\mathcal{T}=K\sigma_{x}$ and have an
unusual Chern number of $\pm 4$. The effective Hamiltonian can be written as
[37]

$H_{2\times
2}(\mathbf{q})=-\begin{pmatrix}Aq_{x}q_{y}q_{z}&B\left(q_{x}^{2}+\omega
q_{y}^{2}+\omega^{2}q_{z}^{2}\right)\\\
B\left(q_{x}^{2}+\omega^{2}q_{y}^{2}+\omega
q_{z}^{2}\right)&-Aq_{x}q_{y}q_{z}\end{pmatrix},$ (5)

where $\omega=\exp(-2\pi i/3)$, and $A$ and $B$ are real constants. We have
omitted the kinetic energy term $f_{0}(\mathbf{q})$ that has nothing to do
with the band topology. The pseudospin texture $\mathbf{S}(\mathbf{q})$ around
the quadruple Weyl point is shown in Fig. 4(a). The $S_{z}$ component of the
pseudospin is illustrated by the colors on the spherical $\mathbf{q}$-surface
in Fig. 4(a) and plotted on a flat map in Fig. 4(b). The direction of
$\mathbf{S}(\mathbf{q})$ is almost always in the $q_{x}q_{y}$ plane except
near $[111]$ and equivalent diagonal directions. As indicated by the blue and
red arrows in Fig. 4(a), the eight diagonal directions can be sub-divided into
two tetrahedron-vertex sets with opposite $S_{z}$, and the total wrapping
number (Chern number) is $|C|=4$.

Figure 4:  (a) The pseudospin texture on a $\mathbf{q}$-sphere enclosing
$\Gamma$, with the colors indicating the magnitude of $S_{z}$ component. The
orientations of $\mathbf{S}$ at 8 diagonal directions are marked especially by
blue and red arrows, showing tetrahedral symmetry. (b) $S_{z}$ component of
the pseudospin as a function of polar and azimuthal angles on the sphere. (c)
Simulated INS intensity near $\mathbf{Q}=(2,4,3)$, showing similar patterns to
(b) since $\mathbf{V}$ is almost in the $z$ direction. All simulations are
done with $|\mathbf{q}|=0.05$ r.l.u. and for the upper band near the 40.1 meV
band-crossing point in MnSi.

An experimentally unfavorable aspect of the quadruple Weyl point is in the
weak dispersion – the energy splitting between the two bands only increases as
a quadratic (in most directions) or even a cubic (along $[111]$ and its
equivalents) function of $\mathbf{q}$. This makes it difficult to separately
measure ${\cal S}(\mathbf{Q},\omega)$ of one of the bands. Moreover, all the
quadruple Weyl points in MnSi and CoSi turn out to be close to other phonon
branches. With the energy resolution of our INS experiment, we are unable to
extract the pseudospin texture of the quadruple Weyl points from the INS data.
Nevertheless, our fitted numerical model have no such resolution limits. In
Fig. 4(c), we calculate ${\cal S}(\mathbf{Q},\omega)$ of the upper band of the
quadruple Weyl point at 40.1 meV in MnSi, on a 0.05 r.l.u.
$|\mathbf{q}|$-sphere surrounding $\Gamma=(2,4,3)$. The BZ center $(2,4,3)$ is
chosen because the $\mathbf{V}$ vector lies very close to the $z$ direction,
such that the INS intensities on the $\mathbf{q}$-sphere exhibit a similar
distribution as the $S_{z}$ component in Fig. 4(b), with four maxima and four
minima approximately along the $\langle{111}\rangle$ diagonal directions. This
“virtual” measurement suggests that INS has in principle the capability to
reveal the pseudospin’s wrapping behavior (and hence the Chern number)
associated with the novel quadruple Weyl points. It can be realized in
experiments with higher resolution and/or in materials where the quadruple
Weyl points are far away from other phonons.

##  V Discussion and conclusion

Our study demonstrates the capability of INS to measure Chern numbers of
topological phonon band crossing. While scattering intensity can vary with
momentum even on non-topological bands, the unique spectroscopic
characteristics of Weyl and Dirac points lie in the fact that they are
singular points for eigenvectors, hence the intensity distribution around them
exhibits abrupt variations: Even on an infinitely small enclosing momentum
surface, the intensity modulations around the topological nodes are still
present, whereas in the case of topologically trivial band crossings, the
intensity distribution would approach a constant as the enclosing surface
shrinks into a point. Consequently, in order to determine Chern numbers from
INS experiments, it is better to study the close vicinity of the band crossing
points, i.e., using a small $\mathbf{q}$-surface to both avoid other bands and
ensure that the intensity modulations arise solely from the topology. We also
note that the intensity modulations we discuss here are universal for two-fold
phonon Weyl points of any Chern number, as well as for some of the four-fold
Dirac points. According to group theory [52], phonon Weyl points exist in
crystals of certain space groups, and they can be found at specific high-
symmetry points in the BZ. As long as they are far away from other bands and
have a relatively large group velocity, intensity modulations should be
observable by INS.

The detection of wave functions (or, vibrational eigenvectors) is not
restricted to INS experiments on phonons, and it can be a common capability of
many spectroscopic methods. In polarization-dependent angle-resolved
photoemission spectroscopy (ARPES), changes in the signal intensity has been
suggested to reflect the wave functions of Dirac electrons [53]. Related
measurements have also been proposed for resonant inelastic x-ray scattering
(RIXS) [54]. While the wave-function texture of the quasiparticles is at the
origin of all the spectroscopic observables, the specific interactions between
the experimental probes and the quasiparticles may add further complexity to
the interpretation of experiments. For instance, INS measurements of phonons
involve neutron collisions with nuclei, giving rise to the
$\mathbf{Q}\cdot\bm{\xi}$ term in the scattering cross section; and because
$\mathbf{Q}$ uniquely determines $\mathbf{V}$ in Eq. (4), the associated
projection of $\bm{\xi}$ always allows for the determination of the Chern
number. On the contrary, for magnon bands, dipole-dipole interactions between
neutrons and magnetic moments lead to a
$\mathbf{Q}\times(\mathbf{Q}\times\mathbf{S})$ term in the cross section, and
the component of $\mathbf{S}$ parallel to $\mathbf{Q}$ is missing from the
detection. As a result, the number of extrema in the intensity modulation on
an enclosing momentum surface may change between different choices of the
measurement BZ, rendering it necessary to have extra knowledge about the
magnetic system in order to correctly infer the topological invariant. In the
cases of ARPES and RIXS, complexity may arise because the different
polarization channels have to be considered together. To this end, our INS
measurement of phonons may be regarded as a demonstration of principles that
motivates further studies.

In conclusion, we have performed a comprehensive study on the topological
phonon band crossings in MnSi and CoSi. Both the band dispersions and the
coherent dynamical structure factors are experimentally resolved with high
precision, yielding results that compare well with model calculations. The
existence of spin-1 Weyl points and charge-2 Dirac points in the phonon bands
are verified by the dispersions. Combining experiments and model calculations,
we further demonstrate the capability of INS for unambiguously determining
Chern numbers of band crossing nodes. Our general theoretical scheme based on
effective Hamiltonians suggests that related methods can be found in the study
of other topological quasiparticles, as well as in the use of other
spectroscopic methods.

###### Acknowledgements.

The work at Peking University was supported by the National Key R&D Program of
China (No. 2018YFA0305602) and the National Natural Science Foundation of
China (Nos. 12061131004 and 11888101). The work at Brookhaven National
Laboratory was supported by Office of Basic Energy Sciences (BES), Division of
Materials Sciences and Engineering, U.S. Department of Energy (DOE), under
contract DE-SC0012704. Part of this research was performed at the MLF, J-PARC,
Japan, under a user program (proposal Nos. 2018A0193, 2018B0201, 2019A0085). A
portion of this research used resources at Spallation Neutron Source, a DOE
Office of Science User Facility operated by the Oak Ridge National Laboratory.

## References

  * Klitzing  _et al._ [1980] K. v. Klitzing, G. Dorda, and M. Pepper, New method for high-accuracy determination of the fine-structure constant based on quantized Hall resistance, Phys. Rev. Lett. 45, 494 (1980).
  * Thouless  _et al._ [1982] D. J. Thouless, M. Kohmoto, M. P. Nightingale, and M. den Nijs, Quantized Hall conductance in a two-dimensional periodic potential, Phys. Rev. Lett. 49, 405 (1982).
  * Qi and Zhang [2011] X.-L. Qi and S.-C. Zhang, Topological insulators and superconductors, Rev. Mod. Phys. 83, 1057 (2011).
  * Burkov [2016] A. A. Burkov, Topological semimetals, Nat. Mater 15, 1145 (2016).
  * Bansil  _et al._ [2016] A. Bansil, H. Lin, and T. Das, Colloquium: Topological band theory, Rev. Mod. Phys. 88, 021004 (2016).
  * Chiu  _et al._ [2016] C.-K. Chiu, J. C. Y. Teo, A. P. Schnyder, and S. Ryu, Classification of topological quantum matter with symmetries, Rev. Mod. Phys. 88, 035005 (2016).
  * Armitage  _et al._ [2018] N. P. Armitage, E. J. Mele, and A. Vishwanath, Weyl and Dirac semimetals in three-dimensional solids, Rev. Mod. Phys. 90, 015001 (2018).
  * Roth  _et al._ [2009] A. Roth, C. Brüne, H. Buhmann, L. W. Molenkamp, J. Maciejko, X.-L. Qi, and S.-C. Zhang, Nonlocal transport in the quantum spin Hall state, Science 325, 294 (2009).
  * Wan  _et al._ [2011] X. Wan, A. M. Turner, A. Vishwanath, and S. Y. Savrasov, Topological semimetal and Fermi-arc surface states in the electronic structure of pyrochlore iridates, Phys. Rev. B 83, 205101 (2011).
  * Yang  _et al._ [2015] L. X. Yang, Z. K. Liu, Y. Sun, H. Peng, H. F. Yang, T. Zhang, B. Zhou, Y. Zhang, Y. F. Guo, M. Rahn, D. Prabhakaran, Z. Hussain, S.-K. Mo, C. Felser, B. Yan, and Y. L. Chen, Weyl semimetal phase in the non-centrosymmetric compound TaAs, Nat. Phys 11, 728 (2015).
  * Lv  _et al._ [2015] B. Q. Lv, H. M. Weng, B. B. Fu, X. P. Wang, H. Miao, J. Ma, P. Richard, X. C. Huang, L. X. Zhao, G. F. Chen, Z. Fang, X. Dai, T. Qian, and H. Ding, Experimental discovery of Weyl semimetal TaAs, Phys. Rev. X 5, 031013 (2015).
  * de Juan  _et al._ [2017] F. de Juan, A. G. Grushin, T. Morimoto, and J. E. Moore, Quantized circular photogalvanic effect in Weyl semimetals, Nat. Commun 8, 15995 (2017).
  * Flicker  _et al._ [2018] F. Flicker, F. de Juan, B. Bradlyn, T. Morimoto, M. G. Vergniory, and A. G. Grushin, Chiral optical response of multifold fermions, Phys. Rev. B 98, 155145 (2018).
  * Haldane and Raghu [2008] F. D. M. Haldane and S. Raghu, Possible realization of directional optical waveguides in photonic crystals with broken time-reversal symmetry, Phys. Rev. Lett. 100, 013904 (2008).
  * Lu  _et al._ [2013] L. Lu, L. Fu, J. D. Joannopoulos, and M. Soljačić, Weyl points and line nodes in gyroid photonic crystals, Nat. Photonics 7, 294 (2013).
  * Lu  _et al._ [2015] L. Lu, Z. Wang, D. Ye, L. Ran, L. Fu, J. D. Joannopoulos, and M. Soljačić, Experimental observation of Weyl points, Science 349, 622 (2015).
  * Xiao  _et al._ [2015] M. Xiao, W.-J. Chen, W.-Y. He, and C. T. Chan, Synthetic gauge flux and Weyl points in acoustic systems, Nat. Phys 11, 920 (2015).
  * Ge  _et al._ [2018] H. Ge, X. Ni, Y. Tian, S. K. Gupta, M.-H. Lu, X. Lin, W.-D. Huang, C. T. Chan, and Y.-F. Chen, Experimental observation of acoustic Weyl points and topological surface states, Phys. Rev. Applied 10, 014017 (2018).
  * He  _et al._ [2018] H. He, C. Qiu, L. Ye, X. Cai, X. Fan, M. Ke, F. Zhang, and Z. Liu, Topological negative refraction of surface acoustic waves in a Weyl phononic crystal, Nature 560, 61 (2018).
  * Li  _et al._ [2018a] F. Li, X. Huang, J. Lu, J. Ma, and Z. Liu, Weyl points and Fermi arcs in a chiral phononic crystal, Nat. Phys. 14, 30 (2018a).
  * Zhang  _et al._ [2018] T. Zhang, Z. Song, A. Alexandradinata, H. Weng, C. Fang, L. Lu, and Z. Fang, Double-Weyl phonons in transition-metal monosilicides, Phys. Rev. Lett. 120, 016401 (2018).
  * Miao  _et al._ [2018] H. Miao, T. T. Zhang, L. Wang, D. Meyers, A. H. Said, Y. L. Wang, Y. G. Shi, H. M. Weng, Z. Fang, and M. P. M. Dean, Observation of double Weyl phonons in parity-breaking FeSi, Phys. Rev. Lett. 121, 035302 (2018).
  * Li  _et al._ [2018b] J. Li, Q. Xie, S. Ullah, R. Li, H. Ma, D. Li, Y. Li, and X.-Q. Chen, Coexistent three-component and two-component Weyl phonons in TiS, ZrSe, and HfTe, Phys. Rev. B 97, 054305 (2018b).
  * Xia  _et al._ [2019] B. W. Xia, R. Wang, Z. J. Chen, Y. J. Zhao, and H. Xu, Symmetry-protected ideal type-II Weyl phonons in CdTe, Phys. Rev. Lett. 123, 065501 (2019).
  * Li  _et al._ [2021] H. Li, T. Zhang, A. Said, Y. Fu, G. Fabbris, D. G. Mazzone, J. Zhang, J. Lapano, H. N. Lee, H. C. Lei, M. P. M. Dean, S. Murakami, and H. Miao, Observation of a chiral wave function in the twofold-degenerate quadruple weyl system BaPtGe, Phys. Rev. B 103, 184301 (2021).
  * Mook  _et al._ [2016] A. Mook, J. Henk, and I. Mertig, Tunable magnon Weyl points in ferromagnetic pyrochlores, Phys. Rev. Lett. 117 (2016).
  * Li  _et al._ [2017] K. Li, C. Li, J. Hu, Y. Li, and C. Fang, Dirac and nodal line magnons in three-dimensional antiferromagnets, Phys. Rev. Lett. 119, 247202 (2017).
  * Yao  _et al._ [2018] W. Yao, C. Li, L. Wang, S. Xue, Y. Dan, K. Iida, K. Kamazawa, K. Li, C. Fang, and Y. Li, Topological spin excitations in a three-dimensional antiferromagnet, Nat. Phys 14, 1011 (2018).
  * Bao  _et al._ [2018] S. Bao, J. Wang, W. Wang, Z. Cai, S. Li, Z. Ma, D. Wang, K. Ran, Z.-Y. Dong, D. L. Abernathy, S.-L. Yu, X. Wan, J.-X. Li, and J. Wen, Discovery of coexisting Dirac and triply degenerate magnons in a three-dimensional antiferromagnet, Nat. Commun 9, 2591 (2018).
  * McClarty [2022] P. A. McClarty, Topological magnons: A review, Annual Review of Condensed Matter Physics 13, 171 (2022).
  * Lovesey [1984] S. W. Lovesey,  _Theory of neutron scattering from condensed matter_ (Oxford University Press, 1984).
  * Squires [2012] G. L. Squires,  _Introduction to the Theory of Thermal Neutron Scattering_, 3rd ed. (Cambridge University Press, 2012).
  * Shivam  _et al._ [2017] S. Shivam, R. Coldea, R. Moessner, and P. McClarty, Neutron scattering signatures of magnon Weyl points (2017), arXiv:1712.08535 .
  * Elliot  _et al._ [2021] M. Elliot, P. A. McClarty, D. Prabhakaran, R. D. Johnson, H. C. Walker, P. Manuel, and R. Coldea, Order-by-disorder from bond-dependent exchange and intensity signature of nodal quasiparticles in a honeycomb cobaltate, Nat. Commun 12, 3936 (2021).
  * Scheie  _et al._ [2022a] A. Scheie, P. Laurell, P. A. McClarty, G. E. Granroth, M. B. Stone, R. Moessner, and S. E. Nagler, Spin-exchange hamiltonian and topological degeneracies in elemental gadolinium, Phys. Rev. B 105, 104402 (2022a).
  * Scheie  _et al._ [2022b] A. Scheie, P. Laurell, P. A. McClarty, G. E. Granroth, M. B. Stone, R. Moessner, and S. E. Nagler, Dirac magnons, nodal lines, and nodal plane in elemental gadolinium, Phys. Rev. Lett. 128, 097201 (2022b).
  * Zhang  _et al._ [2020] T. Zhang, R. Takahashi, C. Fang, and S. Murakami, Twofold quadruple Weyl nodes in chiral cubic crystals, Phys. Rev. B 102, 125148 (2020).
  * Kajimoto  _et al._ [2011] R. Kajimoto, M. Nakamura, Y. Inamura, F. Mizuno, K. Nakajima, S. Ohira-Kawamura, T. Yokoo, T. Nakatani, R. Maruyama, K. Soyama, K. Shibata, K. Suzuya, S. Sato, K. Aizawa, M. Arai, S. Wakimoto, M. Ishikado, S.-i. Shamoto, M. Fujita, H. Hiraka, K. Ohoyama, K. Yamada, and C.-H. Lee, The Fermi Chopper Spectrometer 4SEASONS at J-PARC, J. Phys. Soc. Jpn. 80, SB025 (2011).
  * Granroth  _et al._ [2010] G. E. Granroth, A. I. Kolesnikov, T. E. Sherline, J. P. Clancy, K. A. Ross, J. P. C. Ruff, B. D. Gaulin, and S. E. Nagler, SEQUOIA: A newly operating chopper spectrometer at the SNS, J. Phys. Conf. Ser. 251, 012058 (2010).
  * [40] See supplementary materials for additional measurement data and analysis, .
  * Jin  _et al._ [2022] Z. Jin, Y. Li, Z. Hu, B. Hu, Y. Liu, K. Iida, K. Kamazawa, M. B. Stone, A. I. Kolesnikov, D. L. Abernathy, X. Zhang, H. Chen, Y. Wang, C. Fang, B. Wu, I. A. Zaliznyak, J. M. Tranquada, and Y. Li, Magnetic molecular orbitals in MnSi (2022), arXiv:2206.13699 .
  * Inamura  _et al._ [2013] Y. Inamura, T. Nakatani, J. Suzuki, and T. Otomo, Development status of software “Utsusemi” for chopper spectrometers at MLF, J-PARC, J. Phys. Soc. Jpn. 82, SA031 (2013).
  * Ewings  _et al._ [2016] R. Ewings, A. Buts, M. Le, J. van Duijn, I. Bustinduy, and T. Perring, Horace: Software for the analysis of data from single crystal spectroscopy experiments at time-of-flight neutron instruments, Nucl. Instrum. Methods Phys. Res. Sect. 834, 132 (2016).
  * Xu  _et al._ [2013] G. Xu, Z. Xu, and J. M. Tranquada, Absolute cross-section normalization of magnetic neutron scattering data, Rev. Sci. Instrum. 84, 083906 (2013).
  * Gonze and Lee [1997] X. Gonze and C. Lee, Dynamical matrices, Born effective charges, dielectric permittivity tensors, and interatomic force constants from density-functional perturbation theory, Phys. Rev. B 55, 10355 (1997).
  * Kresse and Hafner [1993] G. Kresse and J. Hafner, Ab initio molecular dynamics for liquid metals, Phys. Rev. B 47, 558 (1993).
  * Kresse and Hafner [1994] G. Kresse and J. Hafner, Ab initio molecular-dynamics simulation of the liquid-metal–amorphous-semiconductor transition in germanium, Phys. Rev. B 49, 14251 (1994).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set, Comput. Mater. Sci. 6, 15 (1996).
  * Kresse and Furthmüller [1996] G. Kresse and J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B 54, 11169 (1996).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996).
  * Shirane  _et al._ [2002] G. Shirane, S. M. Shapiro, and J. M. Tranquada,  _Neutron Scattering with a Triple-Axis Spectrometer: Basic Techniques_ (Cambridge University Press, 2002).
  * Liu  _et al._ [2020] Q.-B. Liu, Y. Qian, H.-H. Fu, and Z. Wang, Symmetry-enforced Weyl phonons, npj Comput. Mater 6, 95 (2020).
  * Hwang  _et al._ [2011] C. Hwang, C.-H. Park, D. A. Siegel, A. V. Fedorov, S. G. Louie, and A. Lanzara, Direct measurement of quantum phases in graphene via photoemission spectroscopy, Phys. Rev. B 84, 125422 (2011).
  * Kourtis [2016] S. Kourtis, Bulk spectroscopic measurement of the topological charge of Weyl nodes with resonant x rays, Phys. Rev. B 94, 125132 (2016).

Supplemental Material for
“Chern numbers of topological phonon band crossing determined with inelastic
neutron scattering”

##  VI Details in experiment

Figure S5: (a) Left side: photograph of CoSi single crystals co-aligned on an
aluminum sample holder. Right side: representative X-ray Laue pattern taken on
a natural surface of one CoSi single crystal along $\langle 100\rangle$
direction. (b) Neutron diffraction intensities of the (1, 1, 2) Bragg
reflection recorded when rotating the entire sample array. The solid line is a
Gaussian fit to the data, with full-width-half-maximum (FWHM) of 1.04∘. (c)
and (d) $\mathbf{Q}$-scan profiles of CoSi indicating that the entire sample
is nearly single-domained, since the (0, 2, 1) reflection is over 25 times
more intense than that of (0, 1, 2). Figure S6: (a) Raw INS spectra of CoSi
along the chosen trajectory of high symmetry directions. (b)
$\mathbf{Q}$-averaged INS intensity over the entire $(2,3,3)$ BZ. The
contributions of calculated coherent, incoherent and fitted aluminum (with
background) scattering intensities are shown respectively. (c) Modified INS
spectra after subtracting incoherent signals. (d) Best-estimated coherent INS
spectra shown in the main text, obtained by subtracting both incoherent and
aluminum signals. Figure S7: (a)–(f) Raman spectra of the optical phonons in
CoSi and MnSi. The $A,E,T$ irreps at $\Gamma$ are identified by utilizing
different polarizations, so as to double-check the degeneracy of each band-
crossing. X’Y’ directions are in 45∘ angle with XY. Figure S8: (a) Gaussian
fits using the same data in Fig. 2(a) in the main text, assuming two/one peak
for the green/red data points. The fits yield three different peak positions
illustrated by the colored dash lines. (b) Energy cuts of the INS data at a
series of $\mathbf{Q}$ points in Fig. 2(c) of main text. The open squares
denote the peak positions obtained from two-peak fits, under the constraint
that equivalent $\mathbf{q}$ positions must have same energies and peak
widths. As indicated by the blue cross, the peak positions are further fitted
collectively into an X-shape linear dispersion near the R point, in
consistance with a charge-2 Dirac point. Data in (a) and (b) are offset for
clarity.

##  VII The fitted-force-constant model based on DFPT

First of all, the force constant matrices in MnSi and CoSi are calculated by
DFPT with vasp. From the force constants, phonon dispersions and eigenmodes
can be obtained by solving a classical, harmonic “spring-ball” model. Although
capturing most features in the phonon dispersions qualitatively, the ab initio
results are not precise enough in predicting certain energy and intensity
values.

To solve this, we first extract the energy values of all phonon branches at
high symmetry points $\Gamma$, X, M, R in the BZ, by visually inspecting INS
spectra and taking average over many BZs. The results are listed in Table S1.
Retrieved phonon energies of MnSi are also shown by red dots in Fig. S9(a) and
(b), together with the intensity spectra. The original phonon dispersions
calculated directly from ab initio outputs are displayed in Fig. S9(b). It is
clear that phonon frequencies are overestimated globally, and $\sim 10\%$
variations are possible at some points in the BZ. The errors still exist
partially after performing a global rescale in energy, which accounts for the
inaccuracy in relaxed lattice constants.

MnSi |  $\Gamma$ |  0.0T0 |  24.3E |  24.8T |  30.2T |  33.6A |  39.4T |  40.0E |  41.4T |  49.7A |  55.9T |  55.9T |  55.9T
---|---|---|---|---|---|---|---|---|---|---|---|---|---
|  $\rm X$ |  17.5 |  19.0 |  24.5 |  26.5 |  28.7 |  31.5 |  34.4 |  36.8 |  40.2 |  42.9 |  53.2 |  55.1
|  $\rm M$ |  20.6 |  21.7 |  24.5 |  25.7 |  30.6 |  31.4 |  33.6 |  36.4 |  39.0 |  42.0 |  53.7 |  54.3
|  $\rm R$ |  21.7 |  24.9 |  29.8 |  34.4 |  40.9 |  53.6 |  |  |  |  |  |
CoSi |  $\Gamma$ |  0T |  23.2A |  25.7E |  28.0T |  30.5T |  38.2T |  38.7E |  43.3T |  45.5A |  52.5T |  |
|  $\rm X$ |  17.0 |  18.5 |  23.2 |  26.3 |  29.7 |  31.3 |  35.5 |  37.6 |  40.9 |  45.7 |  49.2 |  54.2
|  $\rm M$ |  18.0 |  20.3 |  22.9 |  25.5 |  29.5 |  32.1 |  33.9 |  36.4 |  40.6 |  45.9 |  50.3 |  53.7
|  $\rm R$ |  20.7 |  22.8 |  31.5 |  33.5 |  42.8 |  52.7 |  |  |  |  |  |
TABLE S1: Phonon energies (in meV) of all bands at high symmetry points
($\Gamma$, X, M, R), obtained from experimental INS spectra. The superscripts
$A,E,T$ stand for 1, 2, 3-fold degeneracy at $\Gamma$, determined by Raman
spectroscopy. Errors at all points are roughly $0.5$ meV due to the limitation
of resolution. Figure S9: (a) INS spectra of phonons in MnSi. Red dots show
the energy values at $\Gamma$, X, M, R. The black lines illustrate phonon
dispersions simulated with the fitted-force-constant model, in good agreement
with INS results. (b) Unadjusted phonon dispersions from DFPT outputs, with
significant difference with experiments. (c)–(e) Measured and calculated INS
intensities at selected $\mathbf{Q}$-points, showing remarkable accordance.

We found that such deviations can be greatly reduced by adjusting force
constants and performing parametric fits. First, the force constant matrices
between atom $d$ and $d^{\prime}$, $\mathbb{K}_{jj^{\prime}dd^{\prime}}$ are
categorized into equivalent classes according to space group symmetry, and
arranged by their matrix norm $|\mathbb{K}|$ (“strength of springs”). Then,
each of the dominating $N$ classes are multiplied by an adjustable variable
$a_{i}$, while all other insignificant force constants are set to zero. In
this way, the dynamical matrix $\mathbb{D}(\mathbf{q})$ is constructed by

$\displaystyle\mathbb{D}(\mathbf{q})=a_{1}\cdot\sum(\mathbb{K}\text{ in class
\\#1})+a_{2}\cdot\sum(\text{class \\#2})+\cdots+a_{N}\cdot\sum(\text{class
\\#N})$ (S6)

and diagonalized to get the phonon dispersions. Finally, a $\chi^{2}$-fit is
performed on $\\{a_{i}\\}$ to minimize the difference between calculated and
experimental energy datasets. The $\chi^{2}$ is defined as

$\displaystyle\chi^{2}=\sum_{\\{E_{i}\\}}\left(\frac{E_{i}^{\text{calc}}-E_{i}^{\text{exp}}}{\Delta
E_{i}}\right)^{2}$ (S7)

and minimized by the conjugate descend method. All $\Delta E_{i}$ are assumed
to be $0.5$ meV. We also reassured the algorithm’s convergence by testing with
different initial values.

Detailed information of the parametric fits can be found in Table S2. It is
surprising that only very few independent matrices are needed ($N=9$ for MnSi
and $N=11$ for CoSi) to greatly reduce $\chi^{2}$ and reproduce experimental
data to a remarkable extent, either globally [Fig.S9(a)] or at certain
$\mathbf{Q}$-points [Fig.S9(c)-(e)]. Since we set all other force constants to
zero, we are actually trying to find out which types of “springs” are most
important in real crystals. Generally, atom pairs with greater distances have
smaller force constants, i.e. weaker interactions; but it is not always the
case. Ranking by their optimized strengths, fitting parameters $\\{a_{i}\\}$
and atom distances of the leading force constants are listed together in Table
S2.

Material |  $N$ |  $\chi^{2}$ | Optimized bond strengths $|\mathbb{K}|_{i}$ (fitting parameters $\\{a_{i}\\}$)
---|---|---|---
MnSi |  0 |  1462 | ——– (DFPT output)
MnSi |  1 |  112.2 |  0.857 (global rescale)
MnSi |  9 |  56.9 |  17.61 (0.896) |  Si–Si |  0.00 Å |  1.11 (0.719) |  Mn–Mn |  2.79 Å
|  |  |  14.60 (0.810) |  Mn–Mn |  0.00 Å |  1.10 (0.782) |  Si–Si |  2.82 Å
|  |  |  6.93 (0.949) |  Mn–Si |  2.31 Å |  0.79 (1.358) |  Mn–Mn |  4.14 Å
|  |  |  4.25 (0.796) |  Mn–Si |  2.40 Å |  0.51 (1.387) |  Mn–Mn |  4.57 Å
|  |  |  2.12 (0.956) |  Mn–Si |  2.54 Å |  |  |
CoSi |  0 |  465.6 | ——– (DFPT output)
CoSi |  1 |  379.7 |  0.960 (global rescale)
CoSi |  11 |  70.6 |  16.17 (0.989) |  Si–Si |  0.00 Å |  1.37 (1.011) |  Co–Si |  2.44 Å
|  |  |  12.88 (0.909) |  Co–Co |  0.00 Å |  0.88 (1.323) |  Co–Co |  3.97 Å
|  |  |  6.12 (1.477) |  Co–Si |  2.33 Å |  0.62 (0.625) |  Si–Si |  3.87 Å
|  |  |  3.54 (1.038) |  Co–Si |  2.34 Å |  0.52 (1.235) |  Si–Si |  4.47 Å
|  |  |  1.95 (0.815) |  Co–Co |  2.74 Å |  0.35 (0.695) |  Si–Si |  4.45 Å
|  |  |  1.88 (1.057) |  Si–Si |  2.75 Å |  |  |
TABLE S2: Results of the $\chi^{2}$-fits and the optimized parameter sets for
MnSi and CoSi. The optimized strengths (norms of force constant matrice
$|\mathbb{K}|_{i}$, in eV/Å2) of different types of bonds are listed in a
descending order. The fitting parameters $\\{a_{i}\\}$ are shown in blue and
in parentheses. Corresponding bond types and bond lengths are also listed.

##  VIII Determination of Chern numbers from inelastic neutron scattering
intensity

###  VIII.1 Generic two-fold Weyl points

Aside from some unimportant constant coefficients, the coherent dynamical
structure factor of the $i$th phonon mode,
$S_{\text{coh}}^{(i)}(\mathbf{Q},\omega)$, is dominated by the polarization
vector $\bm{\xi}^{(i)}_{d}(\mathbf{q})$[31, 32].

$\displaystyle
S_{\text{coh}}^{(i)}(\mathbf{Q},\omega)\propto|F^{(i)}_{\text{coh}}(\mathbf{Q})|^{2}=\left|\sum_{d}\frac{b_{d,\text{coh}}}{\sqrt{m_{d}}}\mathbf{Q}\cdot\bm{\xi}^{(i)}_{d}(\mathbf{q})e^{i\mathbf{Q}\cdot\mathbf{r}_{d}}\right|^{2}.$
(S8)

In this formula, $\mathbf{Q}=\mathbf{q}+\mathbf{G}$ is the total momentum
transfer, and $m_{d},\mathbf{r}_{d},b_{d,\text{coh}}$ denote the mass,
position and coherent scattering length [51] of the $d$th atom in a single
unit cell. The most important term
$\mathbf{Q}\cdot\bm{\xi}^{(i)}_{d}(\mathbf{q})$ is decided by the interaction
between neutrons and point-like nuclei.

For an arbitrary two-fold band crossing at $\mathbf{q}=0$, we assume that the
two phonon eigenmodes at $\mathbf{q}=0$ are $\bm{\xi}^{(i=1,2)}$. Under this
basis, the effective Hamiltonian of a two-fold Weyl point is denoted by a
$2\times 2$ Hermitian matrix $H_{2\times
2}(\mathbf{q})=\mathbf{f}(\mathbf{q})\cdot\bm{\sigma}+f_{0}(\mathbf{q})\sigma_{0}$.
Then, eigenvector of the upper branch at finite $\mathbf{q}$ is

$\displaystyle\bm{\xi}^{\text{up}}(\mathbf{q})=\psi_{1}(\mathbf{q})\bm{\xi}^{(1)}+\psi_{2}(\mathbf{q})\bm{\xi}^{(2)}.$
(S9)

In the close vicinity of the Weyl point, we consider the limit
$|\mathbf{q}|\ll|\mathbf{G}|$ and $|\mathbf{q}|\ll 2\pi/a$, so that the upper
band intensity is approximately

$\displaystyle S_{\text{coh}}^{\text{up}}(\mathbf{Q},\omega)$
$\displaystyle\propto\left|\sum_{d}\frac{b_{d,\text{coh}}}{\sqrt{m_{d}}}\mathbf{G}\cdot\bm{\xi}^{\text{up}}_{d}(\mathbf{q})e^{i\mathbf{G}\cdot\mathbf{r}_{d}}\right|^{2}$
$\displaystyle=\left|\psi_{1}(\mathbf{q})W^{(1)}(\mathbf{G})+\psi_{2}(\mathbf{q})W^{(2)}(\mathbf{G})\right|^{2},$
(S10)

where $W^{(i=1,2)}$ are defined as

$\displaystyle
W^{(i)}(\mathbf{G})\equiv\sum_{d}\frac{b_{d,\text{coh}}}{\sqrt{m_{d}}}\mathbf{G}\cdot\bm{\xi}^{(i)}_{d}e^{i\mathbf{G}\cdot\mathbf{r}_{d}}.$
(S11)

From the expressions above, it is clear that
$S_{\text{coh}}^{\text{up}}(\mathbf{Q},\omega)$ is nothing but the inner
product of two complex vectors: spinor
$\bm{\xi}^{\text{up}}(\mathbf{q})=\left(\psi_{1}(\mathbf{q}),\psi_{2}(\mathbf{q})\right)^{\text{T}}$
and
$\mathbf{W}(\mathbf{G})=\left(W^{(1)}(\mathbf{G})^{*},W^{(2)}(\mathbf{G})^{*}\right)^{\text{T}}$,
a vector independent of $\mathbf{q}$ but varying among different BZs. From
this formula, the INS intensity would reach maximum when the two vectors are
parallel, i.e.

$\displaystyle\left(\psi_{1}(\mathbf{q}),\psi_{2}(\mathbf{q})\right)=c\left(W^{(1)}(\mathbf{G})^{*},W^{(2)}(\mathbf{G})^{*}\right)$
(S12)

and reach zero when they are perpendicular, i.e.

$\displaystyle\left(\psi_{1}(\mathbf{q}),\psi_{2}(\mathbf{q})\right)=c\left(-W^{(2)}(\mathbf{G}),W^{(1)}(\mathbf{G})\right).$
(S13)

Following the definition of pseudospin operators in the main text, the spinor
$\bm{\xi}^{\text{up}}(\mathbf{q})$ is represented by a unit 3-dimensional real
vector $\mathbf{S}(\mathbf{q})$ on the Bloch sphere

$\displaystyle
S_{i}(\mathbf{q})\equiv\frac{f_{i}(\mathbf{q})}{|\mathbf{f}(\mathbf{q})|}=\left\langle\xi^{\text{up}}\left|\sigma_{i}\right|\xi^{\text{up}}\right\rangle.$
(S14)

Similarly, the complex vector $\mathbf{W}(\mathbf{G})$ can also be represented
by a 3-dimensional real vector $\mathbf{V}(\mathbf{G})$

$\displaystyle
V_{i}(\mathbf{G})=\left\langle\mathbf{W}(\mathbf{G})\left|\sigma_{i}\right|\mathbf{W}(\mathbf{G})\right\rangle,$
(S15)

whose modulus $|\mathbf{V}|$ is not necessarily unitary. In this notation, the
scattering intensity is proportional to the inner product

$\displaystyle S_{\text{coh}}^{\text{up}}(\mathbf{Q},\omega)$
$\displaystyle\propto\left|\psi_{1}(\mathbf{q})W^{(1)}(\mathbf{G})+\psi_{2}(\mathbf{q})W^{(2)}(\mathbf{G})\right|^{2}$
$\displaystyle=\frac{1}{2}(|\mathbf{V}|+\mathbf{S}(\mathbf{q})\cdot\mathbf{V}(\mathbf{G}))$
$\displaystyle\propto|\mathbf{V}|\left(1+\frac{\mathbf{S}(\mathbf{q})\cdot\mathbf{V}(\mathbf{G})}{|\mathbf{S}(\mathbf{q})||\mathbf{V}(\mathbf{G})|}\right).$
(S16)

Now we arrive at the formula described in the main text.

Next, we consider a closed surface $\mathbb{S}_{\mathbf{q}}$ in
$\mathbf{q}$-space that encloses the origin. For two-fold Weyl points, the
upper (lower) band is associated with non-zero Chern number $\pm C$, defined
by the total Berry flux going in or out of the surface. If momentum
$\mathbf{q}$ runs over the surface, the pseudospin-1/2 spinor
$\left(\psi_{1}(\mathbf{q}),\psi_{2}(\mathbf{q})\right)^{T}$ will wrap the
$\mathbb{CP}^{1}$-space exactly $|C|$ times, or equivalently, the pseudospin
$\mathbf{S}$ will wrap the Bloch sphere $\mathbb{S}_{\mathrm{B}}$ exactly
$|C|$ times. Therefore, we conclude that there should be at least $|C|$
momenta on the $\mathbf{q}$-sphere where intensity reaches maximum, and at
least $|C|$ momenta where intensity reaches zero. As mentioned in the main
text, the approximation of $\mathbf{Q}=\mathbf{q}+\mathbf{G}\approx\mathbf{G}$
can also be released since $\mathbf{Q}$ has trivial topology on the surface
$\mathbb{S}_{\mathbf{q}}$ as long as the origin of $\mathbf{Q}$ is not
enclosed.

###  VIII.2 Charge-2 Dirac points

We next generalize this deduction to the Charge-2 Dirac points at R. Although
being four-fold degenerated, the Hamiltonian
$H_{4}(\mathbf{q})\propto\left(\begin{smallmatrix}\mathbf{q}\cdot\bm{\sigma}&0\\\
0&\mathbf{q}\cdot\bm{\sigma}\end{smallmatrix}\right)$ is the direct sum of two
identical spin-1/2 Weyl points, each with Chern number $C=\pm 1$ [21]. So
eigenvectors of the two upper branches are

$\displaystyle\begin{split}\bm{\xi}^{\text{up}}(\mathbf{q})&=\psi_{1}(\mathbf{q})\bm{\xi}^{(1)}+\psi_{2}(\mathbf{q})\bm{\xi}^{(2)},\\\
\bm{\xi}^{\text{up}\prime}(\mathbf{q})&=\psi_{1}(\mathbf{q})\bm{\xi}^{(3)}+\psi_{2}(\mathbf{q})\bm{\xi}^{(4)}.\end{split}$
(S17)

Again we assume $|\mathbf{q}|\ll|\mathbf{G}|$ and $|\mathbf{q}|\ll 2\pi/a$,
and the total intensity of the two upper bands is

$\displaystyle S_{\text{coh}}^{\text{up}}(\mathbf{Q},\omega)$
$\displaystyle\propto\left|\psi_{1}(\mathbf{q})V^{(1)}(\mathbf{G})+\psi_{2}(\mathbf{q})V^{(2)}(\mathbf{G})\right|^{2}+\left|\psi_{1}(\mathbf{q})V^{(3)}(\mathbf{G})+\psi_{2}(\mathbf{q})V^{(4)}(\mathbf{G})\right|^{2}$
$\displaystyle=\left(\psi_{1}(\mathbf{q}),\psi_{2}(\mathbf{q})\right)^{*}\cdot\mathbb{V}(\mathbf{G})\cdot\begin{pmatrix}\psi_{1}(\mathbf{q})\\\
\psi_{2}(\mathbf{q})\end{pmatrix},$ (S18)

where $\mathbb{V}(\mathbf{G})$ is a $2\times 2$ Hermitian matrix

$\displaystyle\begin{pmatrix}|V^{(1)}|^{2}+|V^{(3)}|^{2}&V^{(1)}V^{(2)*}+V^{(3)}V^{(4)*}\\\
V^{(1)*}V^{(2)}+V^{(3)*}V^{(4)}&|V^{(2)}|^{2}+|V^{(4)}|^{2}\end{pmatrix}$
(S19)

It can be written as
$\mathbb{V}(\mathbf{G})=\mathbf{V}(\mathbf{G})\cdot\bm{\sigma}+V_{0}$, where
$V_{0}$ is equal to or larger than $|\mathbf{V}|$ because
$\mathbb{V}(\mathbf{G})$ is semi positive definite. Again, it is also
represented by a 3-dimensional real vector $\mathbf{V}(\mathbf{G})$, but with
an extra constant $V_{0}$. In this case the scattering intensity would be

$\displaystyle S_{\text{coh}}^{\text{up}}(\mathbf{Q},\omega)$
$\displaystyle\propto
2V_{0}+|\mathbf{V}|+\mathbf{S}(\mathbf{q})\cdot\mathbf{V}(\mathbf{G})$
$\displaystyle=2V_{0}+|\mathbf{V}|\left(1+\frac{\mathbf{S}(\mathbf{q})\cdot\mathbf{V}(\mathbf{G})}{|\mathbf{S}(\mathbf{q})||\mathbf{V}(\mathbf{G})|}\right).$
(S20)

where $V_{0}$ is a non-negative constant.

In the case of charge-2 Dirac points, each spin-1/2 copy in the Hamiltonian
has Chern number $\pm 1$. Similar to the arguments above, the pseudospin
$\mathbf{S}$ will wrap the Bloch sphere exactly $|C|$ times on a
$\mathbf{q}$-surface that encloses the origin. Therefore, the total INS
intensity of the two upper bands will have exactly one maximum and two minimum
on the $\mathbf{q}$-sphere. The only difference due to the two spin-1/2
components is that the smallest intensity may be larger than zero.