Topological Signatures in Nodal Semimetals through Neutron Scattering

Topological nodal semimetals are known to host a variety of fascinating electronic properties due to the topological protection of the band-touching nodes. Neutron scattering, despite its power in probing elementary excitations, has not been routinely applied to topological semimetals, mainly due to the lack of an explicit connection between the neutron response and the signature of topology. In this work, we theoretically investigate the role that neutron scattering can play to unveil the topological nodal features: a large magnetic neutron response with spectral non-analyticity can be generated solely from the nodal bands. A new formula for the dynamical structure factor for generic topological nodal metals is derived. For Weyl semimetals, we show that the locations of Weyl nodes, the Fermi velocities and the signature of chiral anomaly can all leave hallmark neutron spectral responses. Our work offers a neutron-based avenue towards probing bulk topological materials.

For some time, the discovery of new topological materials was largely conducted on a case-by-case basis, as was the case for the prototypical Bi 2 Se 3 TI family (with members Sb 2 Te 3 and Bi 2 Te 3 ) and the TaAs WSM family (with members TaP, NbP, and NbAs). Recent theoretical progress in symmetry-based indicators [42][43][44] and topological quantum chemistry [45] has profoundly accelerated the theoretical search of nontrivial topology within materials. Through symmetry-based indicators, bandstructures are interpreted in terms of elementary basis states whereby the topologically-trivial atomic insulator states are removed allowing the nontrivial topological bands to be revealed. On the other hand, topological quantum chemistry provides a holistic description of The fact that approximately 30% of all 3D stoichiometric materials may carry nontrivial topology brings us to the question of experimental verification. In contrast to the extensive developments on the theoretical forefront, the experimental determination of nontrivial topology remains an arduous challenge. Angular-resolved photoemission spectroscopy (ARPES) plays a significant role in probing topology by directly measuring the bandstructures and afterward, comparing with ab initio calculations, but a few difficulties exist [50][51][52]. Fundamentally, only materials that can cleave without difficulty can be probed by ARPES. This limits the types of materials that can be measured with this technique, since many materials do not form the required atomically-flat and well-ordered surface following cleavage. On top of this, the signals from photoelectrons prohibit the application of a magnetic field to study field-driven phenomena. From a practical viewpoint, the high vacuum and the atomically flat surface pose high technical barriers for sample preparation. Other experimental probes, such as scanning tunneling microscopy (STM), pose an even higher technical impediment due to the necessity of protecting against minuscule vibrations, while quantum transport, despite being routinely carried out, is generally considered as an indirect method to infer bulk topology in TNMs [53,54]. Given the large number of topological material candidates, an alternative experimental approach that can address the aforementioned challenges of existing techniques, but still capable of capturing the non-trivial topology, would be highly desirable.
In this work, we demonstrate that neutron scattering is a powerful, yet unexplored, tool to capture the topological nature of TNMs (Fig. 1). Through quantum singlebody calculations, we derive a new dynamical structure factor formula for generic TNMs, showing that they can have a large magnetic neutron scattering response. This solely arises from the nodal bands which potentially carry nontrivial topology. This contrasts with the conventional Fermi liquid in which the magnetic signal from itinerant electrons is negligibly small [55]. In particular, a closedform expression for the dynamical structure factor is obtained for prototypical WSMs, whereby we show that the Weyl node locations in the Brillouin zone and the magnitude of the Fermi velocities will leave hallmark signatures in neutron spectra through their momentum space separation and their chemical potential. Moreover, when both electrical and magnetic fields are applied, we show that the chiral anomaly can leave measurable fingerprints on the spectra. These features retain for a realistic WSM containing multiple Weyl nodes and are generalizable to other types of TNMs which exhibit nodal lines or surfaces. Given that neutron scattering can directly measure bulk samples and is compatible with external fields, our proposed work offers a new avenue towards measuring TNMs with advantages such as the ability to probe a broader repertoire of TNMs which are only sensitive to topological bands, and to study possible field-driven topological phase transitions in these exotic materials. The rest of the paper is organized as follows. In Section II, we begin by deriving a formula of the dynamical structure factor for a generic TNM and subsequently, in Section III, we expand on this to compute this quantity for a toy model of Weyl fermions. A model in the case of a nodal line semimetal is provided in Section IV. We then calculate the corresponding neutron response function for realistic WSMs with multiple pairs of Weyl nodes in the Brillouin zone through which the neutron response can reveal their location (Section V) and offer an approach to experimentally observe a signature of chiral anomaly (Section VI). Effects of finite mass in DSMs and finite temperature on the neutron response are investigated in Sections VII and VIII, respectively. We holistically discuss our results in Section IX on the large magnetic neutron scattering response resulting from nontrivial topology in TNMs with a WSM example and possible candidate materials. We conclude in Section X with a general discussion on how our results fit in with current experimental techniques used to probe TNMs. Appendices contain further details of calculations.

II. DYNAMICAL STRUCTURE FACTOR FOR GENERIC TOPOLOGICAL NODAL METALS
We start with a general linear-quadratic Hamiltonian which can be used to describe TNMs where ψ † (r) (ψ(r)) creates (annihilates) an electron at position r; α, β denote general spin and band indices; i, j = 1, 2, 3 label spatial directions; is the reduced Planck constant; v αβ i is the nodal term of the generalized Fermi velocities; a αβ i is the node shift; and b αβ ij is a conventional quadratic term.
To describe the influence of a magnetic field, we adopt the minimal coupling prescription i ∂ i → i ∂ i + eA i in Eq. (1). Within the linear-response regime, the current operator in momentum space can be written as (2) where V is the volume, k is the wavevector and Q ≡ k i − k f is the momentum transfer between initial and final states which is also the neutron momentum transfer.
Following the derivation sketched in Appendix A, the dimensionless magnetization obtained from the current operator in the case of a TNM can be written as where µ B is the Bohr magneton in Gaussian units and Q = |Q|. The corresponding magnetic dynamical structure factor tensor can be written into the following form where ω is the frequency of the neutron and M l (Q, t) = e iHt/ M l (Q, 0)e −iHt/ is the neutron magnetization operator in the Heisenberg picture. By comparison, the magnetization operator comprising of localized moments and itinerant electrons is written as [56] where the sum runs over all local sites N and s N , k e,N are the local spin and electron wavevector, respectively. For non-relativistic electrons, we have J = −e k/m e where m e is the electron mass. Our result in Eq. (3) for TNMs thus resembles the magnetization for itinerant electrons in Eq. (5). However, there is one major difference between Eqs. (3) and (5) that enables TNMs to induce a large neutron scattering signal. In conventional Fermi liquids (second term of Eq. (5)), the total magnetization that arises from different itinerant electrons k e,N cancels out, leading to a small neutron magnetic scattering signal [55]. To understand this, taking the Q → 0 limit where magnetic signal is large, the magnetization is M ∼ i N k e,N /Q, and the magnetic signal from an electron with momentum k e,N cancels out the −k e,N contribution from another electron. In contrast, for TNMs, topology plays a role to protect the gapless nodal features, in which different charge carriers can contribute to the magnetization in an additive manner (v αβ i term in Eq. (2)). This leads to a large non-cancellation effect of the corresponding magnetization and thus, to a large neutron signal. For instance, in a simple WSM, the corresponding v αβ i term that generates the current is simply the Fermi velocity of the Weyl fermions.
To facilitate the calculations, we rewrite the dynamical structure factor tensor of Eq. (4) as where β = (k B T ) −1 , k B is the Boltzmann constant, T is the temperature, jik is the Levi-Civita symbol, and Π im (Q, ω) is the dynamical response function (also known as the magnetic susceptibility tensor or the current-current correlation function) defined as where ω n = 2πn/β (n = 0, ±1, ±2, . . . ) is the bosonic Matsubara frequency. Eq. (7) originates from the imaginary-time current-current correlation function where T τ is the imaginary time-ordering operator, . . . denotes the imaginary time average over the whole canonical ensemble, and the total current operator is given in the Heisenberg picture as J m (Q, τ ) = e Hτ J m (Q)e −Hτ . Further details of the calculation of Eq. (6) can be found in Appendix B.
Having defined the dynamical structure factor, we can write the final form of the double differential cross-section as follows [56] with γ = +1.913 being the dimensionless neutron gyromagnetic ratio and r e = e 2 /m e c 2 , the classical electron radius. An important consequence of Eq. (9) is that the poles in the susceptibility Π im (Q, ω) can be reflected as elementary excitations peaks in the scattering spectra.

III. MAGNETIC NEUTRON RESPONSE IN PROTOTYPICAL WEYL SEMIMETALS
The aim of this section is to compute the dynamical structure factor S jl (Q, ω) for a toy model consisting of a single Weyl node using the formalism of Section II. The low energy effective Hamiltonian for a single massless Weyl node is given by is the vector of Pauli matrices which operates on the pseudospin degrees of freedom, and µ is the bare chemical potential. Note that the eigenstates of H have definite chirality. From Eq. (10), the charge density operator is given by ρ Q = −(e/V ) kαβ ψ † kα ψ k−Qβ and the corresponding current operator is Generally, the dynamical response function will depend on both the magnitude and the direction of the wavevector Q. However, for isotropic Weyl nodes with rotational symmetry, we choose Q = Qẑ without loss of generality. Correspondingly, with this choice of Q direction, only the perpendicular components of the dynamical structure factor tensor S xx (Q, ω) and S yy (Q, ω) will contribute to the differential cross-section in Eq. (9). Indeed, by using Eq. (6), we obtain from which we can see that only the components of the response function that are perpendicular to that of the dynamical structure factor need to be computed. For instance, only Π xx (Q, ω) and Π zz (Q, ω) components contribute to S yy (Q, ω). Using Eq. (3) and the fact that Q is chosen is lie along theẑ axis, Π zz (Qẑ, ω) does not contribute to the overall magnetization. Thus, for a component of the current along thex axis (Fig. 1), the magnetic neutron scattering signal will originate from S yy (Qẑ, ω) which, for its part, only possesses a contribution from Π xx (Qẑ, ω).
We explicitly compute the current-current correlation function Π xx (Qẑ, ω) from the effective low energy Hamiltonian in Eq. (10). The corresponding non-interacting Green's function of a single Weyl node in the Matsubara frequency domain is given by where I is the 2 × 2 identity matrix. The non-interacting current-current correlation function is given as By summing over the Matsubara frequency and performing analytical continuation, we arrive at where the summation is over the sign denoted by s, s ; η → 0; n F (E) = (e βE + 1) −1 is the Fermi-Dirac distribution function, and the prefactor h ss is defined as To proceed, we first present the zero temperature calculation of Eq. (15) and leave the finite temperature result for Section VIII. Only the imaginary component of the response function is presented as it contributes to the measured neutron signal. In this derivation, we focus on the case of ω > 0 without loss of generality as the ω < 0 case can be obtained by taking the complex conjugate (i.e., Π(Q, −ω) = [Π(Q, ω)] * ). The currentcurrent correlation function can be decomposed into two separate pieces as corresponding to an undoped contribution Π UD (Qẑ, ω), where the chemical potential is located at the Weyl node, and a doped contribution Π D (Qẑ, ω), where the chemical potential is situated at a finite value away from the Weyl node. The current-current correlation function can be further separated into interband and intraband contributions as illustrated in Fig. 2 with the inclusion of a negative contribution stemming from extra, yet forbidden, transitions of the interband type [57,58]. Schematic of Weyl cones with chemical potentials µ1, µ2 above (purple) and below (orange) the Fermi level. The first contribution corresponds to interband transitions with µ = 0; the second, to intraband transitions in the upper (purple) and lower (orange) bands; and the third, as a negative contribution from forbidden transitions that arose from interband transitions when setting µ = 0.
For the undoped portion, the contribution is solely from s = −s = ±1 interband transitions and we obtain where Θ(x) denotes the Heaviside step function and the doped contribution is obtained as where the summation is over the sign denoted by s, s and for conciseness, we define auxiliary functions g (±) and Rect(x) as By substituting Eqs. (16) and (17) into Eq. (6), the dynamical structure factor can be written as This is valid for a single Weyl node; for multiple Weyl nodes, the dynamical structure factor is obtained by summing over the independent intra-node, interband contributions of each Weyl node. Eq. (18) contains rich physical meaning as shown in Fig. 3. First, the dependence of S yy on v F suggests that a larger Fermi velocity can contribute larger signal. Second, the term involving Θ(−v F Q + ω) indicates the presence of a large and discontinuous signal, here in the second derivative of S yy located along the line ω = v F Q. Such discontinuity signifies the physical onset for the interband electronic transition within the Weyl node. There are additional discontinuities along ω = v F Q ± 2µ for finite chemical potential. This leads to a hallmark "triplet" signature equally spaced in both energy and wavevector space. It is also worthwhile mentioning that although the relation ω = v F Q resembles the form of linear phonon dispersion, the neutron energy transfer ω and momentum transfer Q are linked by the Fermi velocity v F of the Weyl fermions. This offers a clear signature as to how neutron scattering can probe the value of the Fermi velocity in a nodal system. Since a time-of-flight neutron scattering measurement spans a four-dimensional (Q, ω) space, the neutronbased approach can readily determine the Fermi velocity vector for anisotropic Weyl cones. . Dynamical structure factor at zero temperature. Syy (in arbitrary units) as a contour plot against wavevector Q and frequency ω at zero (left) and finite (right) chemical potential taken to be 25 meV. Line cuts of Syy at a constant frequency are shown below the contour plots along with values of the second derivative with respect to Q. There is a large discontinuity in the second derivative located along ω = vF Q (black dotted line) and additional ones along ω = vF Q ± 2µ (purple dotted lines) at finite chemical potential, leading to an overall "triplet" signature.

IV. RESPONSE FUNCTION OF A NODAL LINE SEMIMETAL
The response function can be generalized for the case of a NLSM by adopting the Hamiltonian H = i d i σ i (i = 1, 2), where σ i are the Pauli matrices and m is the mass, p 0 is a constant, k i are momentum components and we define d(k) = i d 2 i (k). We define the Green's function We are interested in computing the correlation function of Eq. (7), using the expression for the Green's function for a NLSM (Eq. (20)). After defining the resulting currents from the Hamiltonian in second quantization (Eq. (C1) in Appendix C), we compute the polarization operator, analogous to Eq. (14), as where In Appendix C, we show a detailed calculation of the response function by substituting the explicit expression for the Green's function (Eq. (20)) into the polarization operator (Eq. (21)) followed by performing the Matsubara summation and momentum integration. As shown in Eq. (6), since we are interested in the neutron dynamical structure factor, we take the imaginary part of the response function Im [Π xx (Q, ω)]. This leads to four separate contributions corresponding to intraband and interband transitions within the valence and conduction bands from k to k + Q (Eq. (C4)).
Here, for simplicity, we compute only one component in Im [Π xx (Q, ω)] for the case of zero Fermi level and zero temperature. Only one contribution out of the four, associated with the transition between states with energies −d(k) and d(k + Q), survives (i.e., n F (−d(k)) = 1 and n F (d(k + Q)) = 0). We assume that Q = Qẑ such that only the Π xx component is under consideration, as was done previously. After the Matsubara summation, the resulting expression is where the S (±) (k) matrices are defined in Eq. (C3). From Eq. (22), we can remove the δ-function by taking an integral over k z . By setting the argument of the δfunction to zero, we can obtain expressions for k z and subsequently, for d(k) and d(k + Q) at this value of k z . Setting the argument of the delta function to zero and solving it requires the condition that ω 2 − v 2 F Q 2 > 0 thereby determining the limits of integration as shown in Appendix C. We compute the remaining integral over k ⊥ with limits of integration that are defined by conditions imposed by the δ-function in Eq. (C4) to obtain the response function for a NLSM at µ, T = 0 as where we define an auxiliary function ξ as .
The resulting neutron dynamical structure factor S yy (Qẑ, ω) given through the imaginary part of the response function of Eq. (23) also displays non-analytical behavior at ω = v F Q resulting in a similar plot to Fig. 3 (with µ is taken to be zero). Analogous to Section III, the effect of nonzero chemical potential in NLSMs will introduce additional discontinuities at ω = v F Q ± 2µ, which can revealed through a similar procedure of calculations. Ultimately, this captures the relationship between the energy transferred by the neutron ω and the Fermi velocity v F of the electronic quasiparticles within a NLSM. This link can be inferred through neutron scattering experiments to extract the value of v F , much like the case of a simple WSM.

V. USING NEUTRON PROBES TO EXPLORE WEYL NODES
We proceed to calculate the dynamical structure factor for a realistic WSM with multiple pairs of Weyl nodes in the Brillouin zone. Without loss of generality, we assume the Weyl nodes to be located at momentum space locations k W,ζ with Fermi velocity v F,ζ and chemical potential µ ζ , where ζ = 1, 2 labels the subset of Weyl node. The low-energy Hamiltonian for each Weyl node with explicit momentum space position follows from Eq. (10) with the shift ∇ → ∇ − ik W,ζ . Here we assume that the two Weyl nodes in consideration share the same chirality. The corresponding Matsubara Green's function for each This is the same result as Eq. (13) with a shift due to the position of the Weyl node as k ζ = k − k W,ζ and with proper indexing. If intra-Weyl node scattering is solely considered, the calculated dynamical structure factor between two Weyl nodes is actually a sum of two contributions originating from taking into account one single Weyl node. Inter-Weyl node scattering may also occur between two Weyl nodes at distinct positions k W,1 and k W,2 (with chemical potentials µ 1 , µ 2 , respectively). In this circumstance, analogous to Eq. (14), the currentcurrent response function for inter-Weyl node scattering with G ζ (k, iω n , µ ζ ) defined in Eq. (24). Here, as a proof of principle, we do not consider a specific interaction potential that can lead to the inter-node scattering, nor the effect from the finite relaxation time. This assumes that strong inter-Weyl node scattering can be determined from phenomena like the Kohn anomaly [59]. To clarify the experimental signature of inter-Weyl node scattering, we assume that the two Weyl nodes share the same isotropic Fermi velocity, i.e., v F = v F,1 = v F,2 . Similar to the intra-Weyl node scattering of Section III, the current-current response function of Eq. (25) can be split into an undoped (µ 1 = µ 2 = 0) and a doped (µ 1 , µ 2 = 0) contribution as The imaginary part of the undoped contribution is computed as represents the momentum transfer with the consideration of inter-Weyl node scattering andω = ω − µ 1 + µ 2 is the associated energy transfer. Note thatω can be both positive and negative depending on the relative importance between the incoming neutron frequency and the chemical potentials of the Weyl nodes. If we define an auxiliary functionf (with associated auxiliary functions g (±) (x) and Rect(x) defined previously) given bỹ then the imaginary part of the doped contribution is The firstf function in the expression for Im Π int xx,D (Q, ω, µ 1 , µ 2 ) produces a nonzero contribution only when the condition µ 1 , µ 2 > 0 is satisfied, whereas the secondf function is nonzero only for µ 1 , µ 2 < 0. These two terms correspond to the electron and the hole responses, respectively.
The corresponding structure factor S yy (Qẑ, ω) originating from both undoped and doped contributions of the neutron response function for the case of inter-Weyl node scattering, at zero temperature, is computed as Due to the presence of terms involving the Heaviside function in Eq. (26), including those inf , one expects a large non-analytical neutron signal atQ =ω or, equivalently, at v F Q = ± v F (k W,2 − k W,1 ) +ω as a result of inter-Weyl node scattering, as seen in Fig. 4. Since one usually maintains control overω through a priori knowledge of the chemical potential of the Weyl nodes and selection of the incoming neutron energy, the measured neutron signal can serve as a probe for the momentum space location of the distinct Weyl nodes within the material system. Measuring the neutron signal with suitable instruments such as triple-axis and time-of-flight neutron spectroscopies followed by performing derivatives of the signal with respect to Q reveals the non-analytical behavior. Since the latter forms along the line Q ≈ k W,2 −k W,1 , it is readily discernible in the neutron data and provides information about the different nesting conditions between the nodes that stem from the material system, thereby revealing their momentum space locations.

VI. SIGNATURE OF CHIRAL ANOMALY IN NEUTRON SPECTRA
In this section, we focus on an important phenomenon associated with WSMs which is the chiral anomaly [60][61][62][63]. It is observed when a pair of parallel electric and magnetic fields, E and B, induces a non-local pumping from one node possessing positive chirality to another with negative chirality for E · B > 0 (and vice-versa for E · B < 0). The charge transfer is stabilized by an inter-Weyl node scattering mechanism with characteristic relaxation time given by τ . The low-energy Hamiltonian in the vicinity of a Weyl node is given by Eq. (10), but unlike Section V, we reintroduce chirality λ along with the shift in momentum position from the origin of the Weyl nodes by taking ∇ → ∇ − iλk W and µ → µ λ with λ = ±1. The chiral chemical potential µ λ is given by [64,65] where µ is the chemical potential without any consideration of chirality and τ is the internode scattering timescale. Eq. (27) is only valid in the limit of weak magnetic field whereby the discreteness of Landau levels can be ignored. The current-current response function in WSMs in the context of the chiral anomaly has been previously studied in Ref. [65]. It is mentioned that for the chiral anomaly to be observable in experiments, the shift in the chemical potential seen in Eq. (27) should be on the order of the chemical potential µ itself. By taking typical values of v F and µ for WSMs, it is shown that the required value of E · B for observation is attainable in present-day experiments. Here, we follow a similar approach, but apply the results towards the calculation of the dynamical structure factor for neutron scattering experiments. The corresponding Matsubara Green's function in the vicinity of the Weyl node with chirality λ is given by Eq. (13) with k → k λ = k+λk W . The response function for two Weyl nodes with opposite chirality can be obtained by simply summing the contributions from the two Weyl nodes with a modified chemical potential µ → µ λ . The imaginary part of the current-current response function can be found in terms of this modified chemical potential. The corresponding dynamical structure factor at zero temperature is obtained as Similar to the discussion in Section III, at a constant value of wavevector (which is a common experimental measurement for triple-axis neutron spectroscopy), the second derivative in the dynamical structure factor S yy with respect to the frequency ω displays divergence-like behavior at ω = v F Q and at ω = v F Q ± 2µ λ . In this case, the chiral anomaly shifts the chemical potential from µ → µ λ in accordance with Eq. (27) due to an extra term that depends on E · B. Hence upon application of parallel electric and magnetic fields, one expects that the divergence in the second derivative will also be shifted and occur at ω = v F Q ± 2µ λ whereas the divergence at ω = v F Q remains unaltered. This behavior is illustrated in Fig. 5 and serves, by itself, as a novel experimental signature of the chiral anomaly in WSMs. Provided that the shift µ → µ λ is on the order of the original value of the chemical potential, the shift in ω whereupon the divergence in the second derivative occurs can be measured as a function of the angle θ between E and B to unveil information on v F and remarkably, on τ , the latter of which is generally difficult to quantify in experiment. . These values of µ λ correspond to a shift from the original value µ induced by the chiral anomaly as implied by the E · B term in Eq. (27). The E · B-induced shift is taken to be on the same order as µ 3 where µ is set to 25 meV and positive chirality is assumed. The right-most plot shows the second derivative of Syy with respect to the frequency ω to highlight divergent behavior at ω = vF Q and ω = vF Q ± 2µ λ , the latter of which is heavily influenced by the E · B shift.

VII. MASSIVE DIRAC FERMIONS
We investigate how introducing a finite mass to the massless Weyl fermions of Eq. (10) induces a change to the dynamical structure factor. It is convenient to rewrite Eq. (10) into a massive four-band Hamiltonian in terms of the anti-commuting four-component γ matrix [66] whereψ ≡ ψ † γ 0 is the Dirac conjugate spinor field, γ = (γ 1 , γ 2 , γ 3 ), γ 5 = iγ 0 γ 1 γ 2 γ 3 , and (k w ) µ = (k 0 w , k w ) is a 4-vector. In this derivation, we set = v F = 1 for simplicity and employ the Feynman slash notation, wherein / W = γ µ W µ for the W µ 4-vector. We see that ψ/ k w γ 5 ψ term leads to an induced Chern-Simons term of the form (k w ) µ µνλσ F λσ A ν in a WSM [61,67]. The chiral shift k w is present in the free Hamiltonian of WSMs. Additionally, this shift can also be dynamically generated in DSMs in the normal phase with the presence of a magnetic field [68]. Similar to the case of graphene, the generation of the Chern-Simons term implies a topolog-ical nature of the normal state in this material [69]. We compute the current operators of the Hamiltonian in Eq. (29) to obtain The first current operator in Eq. (30) is known as the gauge current with global U(1) gauge invariance and the second, an axial (chiral) current operator. As a result of the chiral magnetic effect, the axial current density J 5 µ is generated in the free theory when a fermion charge density and a magnetic field are present [70]. From the expression of J 5 µ , the chiral shift k w is already induced to the lowest order in perturbation theory even if k 0 w = 0. As a result, a DSM with a nonzero charge density is transformed into a WSM immediately when an external magnetic field is applied to the system.
We proceed to calculate the current-current response function using the relativistic-like model of Eq. (29) and the current operators of Eq. (30). This response function is written in relativistic notation [66] as where iG −1 (k) = / k − / k w γ 5 − m and d 4 k/(2π) 4 denotes an integral over the four-dimensional Euclidean momentum space.
When we rewrite the propagator G −1 (k) in terms of k µ w , we obtain three contributions: the first one is proportional to the vacuum polarization of quantum electrodynamics, while the second and third ones are the contributions that are linear and quadratic in k µ w , respectively. Although nonlinear terms may introduce further interesting phenomena, we restrict our analysis to the term that is linear in k µ w such that The Π 0 im (Q) is the usual vacuum polarization tensor of quantum electrodynamics [71]. The Π 1 im term that is linear in k µ w is given by where iG −1 0 (k) = / k − m is the bare propagator at k µ w = 0 and G 1 (k) = i / k− / kwγ 5 −m / k w γ 5 G 0 (k). We compute the imaginary part of the Π 0 im term in the whole domain of Q 2 following analytical continuation [71] as is the squared norm, and η im is the metric tensor. In Eq. (34), we considered the case corresponding to µ = 0 and the dimensional regularization scheme was implemented [72]. For Q 2 < 4m 2 , the imaginary part of the vacuum polarization function is zero. The currentcurrent response function for the undoped case in Section III, shown in Eq. (16), corresponds to the polarization bubble diagram of quantum electrodynamics in Eq. (34) with m = 0. We calculate the imaginary part of Π 1 im to first order in k µ w as [73] Im We can write the dynamical structure factor at zero temperature using Eqs. (6), (34) and (35) to obtain

VIII. EFFECT OF FINITE TEMPERATURE
We calculate the dynamical structure factor of a Weyl node through the current-current response function with consideration of finite temperature. For intra-Weyl node scattering (Section III), the charge neutral (undoped) part of the response function Π xx,UD (Q, ω) remains unaffected whereas the doped contribution Π xx,D (Q, ω) is modified by finite temperature. From the Matsubara summation of Eq. (13), we obtain the finite temperature expression with auxiliary functions G ± and H defined as In particular, it can be shown that under the T → 0 limit, Eq. (37) reduces to the zero temperature result of Eq. (17). The final form of the dynamical structure factor at finite temperature can be written as Figure 6. Dynamical structure factor at finite temperature. Syy (in arbitrary units) as a contour plot against wavevector Q and frequency ω at different temperatures of 1K, 5K, 100K, and 300K. The chemical potential is taken to be 25 meV. Line cuts of Syy at a constant frequency are shown below these contour plots along with values of the second derivative with respect to Q. There is a large and discontinuity in the second derivative located along ω = vF Q (black dotted line). The discontinuities associated with chemical potential at ω = vF Q ± 2µ (purple dotted lines) are smeared out by the effect of finite temperature with a small remnant signal upwards to 10K (as suggested by the two left-most subplots).
The effect of finite temperature from Eq. (38) is shown in Fig. 6. As mentioned in Section III, at zero temperature, three discontinuities in the second derivative were observable: at ω = v F Q, at ω = v F Q + 2µ and another branch at ω = v F Q − 2µ, the latter of which may contribute to plasmon excitations [74]. At finite temperatures, the branches associated with the chemical potential of the Weyl mode remain visible at low temperatures, but are quickly smeared out as the temperature increases as observed in Fig. 7. In contrast, the contribution from the Weyl fermion collective excitation persists as the discontinuity along the line ω = v F Q lingers even at higher temperatures. This highlights the suitability of these signals to be measured in both time-of-flight and triple-axis neutron spectroscopies within experimental conditions.

IX. DISCUSSION
The previous sections illustrate how the nontrivial neutron response of TNMs is specific to these types of exotic materials and can be used to probe for topological phenomena. Remarkably, calculations of the dynamical structure factor for different models of TNMs unveil how the neutron scattering response displays non-analytical behavior arising from the existence of topological nodes in momentum space in addition to unique electromagnetic responses such as the chiral anomaly in the case of WSMs.
In TNMs, the non-analytical smoothness of the neutron scattering signal with respect to wavevector or energy lies along a line ω = v F Q which persists even at high Figure 7. Temperature dependence of second derivative. Second derivative of Syy (in arbitrary units) with respect with wavevector Q at different temperatures. The chemical potential is taken to be 25 meV. Each line plot corresponding to a temperature is vertically offset by 175 for clarity. The divergence in the second derivative at ω = vF Q ± 2µ (left and right dotted lines) decreases with increase of temperature, whereas the signal at ω = vF Q (central dotted line) prevails up to room temperature.
temperatures, demonstrating topological robustness as well as a pointer of how one can use this feature to explore the topological nodes in momentum space and determine v F . Furthermore, this characteristic non-analytic smoothness may also exist along lines ω = v F Q ± 2µ for finite µ so long as the temperature remains sufficiently low as this signal diminishes with temperature increase. This can be used to experimentally determine parameters associated with the chiral anomaly in WSMs provided that the temperature is sufficiently low.
As a proof of concept, we apply the calculations made in the previous sections for the case of type-I WSM tantalum phosphide (TaP). TaP crystallizes in a body-centered tetragonal lattice with space group I4 1 md (109) and point group C 4ν with a = b = 3.32Å and c = 11.34Å and lacks inversion symmetry [75]. Fig. 8 illustrates the calculation of the second derivative of S yy with respect to wavevector Q, averaged along all directions of Q, within the Brillouin zone of TaP. These calculations reveal a large intensity of this signal in regions that correspond to combinations of nesting wavevectors between the locations of the Weyl nodes (corresponding to divergences at ω = v F Q). Elevated intensities can also be coincidentally seen near the locations of the Weyl nodes themselves due to the crystal symmetry. For simplicity, the calculation does not account for the chemical potential difference between the W1 and W2 subset of nodes in this material which would further reveal elevated intensities at ω = v F Q ± 2µ. Figure 8. Calculation for type-I WSM TaP. Intensity map of the second derivative of the dynamical structure factor Syy with respect to wavevector Q (averaged along all directions in Q-space) in the Brillouin zone of type-I WSM tantalum phosphide (TaP). Elevated intensity is seen at locations corresponding to combinations of nesting wavevectors amongst the different locations of Weyl nodes, the latter illustrated as circles in the plot with indication of subset (blue or magenta) and chirality (open or closed). The calculation shown here does not consider for the chemical potential difference between two Weyl nodes which would lead to higher intensities at ω = vF Q ± 2µ.
By performing inelastic neutron scattering with a triple-axis or a time-of-flight experiment and calculating the second derivative of the measured signal, one can extract the nesting wavevectors between the Weyl nodes by examining the location of these elevated intensity points in momentum space and reconstruct the locations of the nodes themselves. This could effectively provide experimental parameters of the Weyl nodes such as momentum space locations and chemical potential.
The calculations performed in this study have mainly focused on isotropic systems. However, several DSMs and WSMs have been shown to exhibit anisotropic topological points. This anisotropy in the response function can be captured in our calculations of S yy by performing a similar substitution for the wavevector Q → Q as described in Ref. [65] which depends on the anisotropic velocities {v x , v y , v z } of the Dirac or Weyl cone. As the magnetic neutron scattering signal depends on the component of the current taken to be in thex direction in the picture described in Section III, one can orient the sample such that the largest velocity component points in this direction to obtain an increased signal. Conversely, one can also map out the magnitudes of the Fermi velocity components using different orientations of the sample if this information is not known a priori. This study highlights the use of neutron scattering for probing exotic topological phenomena. In fact, the advantageous opportunity of concurrently applying external magnetic and electric fields during the measurement, which is not applicable with the aforementioned techniques and has possibly enhanced resolution capabilities, positions inelastic neutron scattering as a possible enabler of further discoveries in the realm of TNMs. The direct probing of electronic bands with neutron scattering also differs from probing topology with its lattice degrees of freedom [59,76,77].
As the repertoire of TNMs continually enlarges, the requirements of ultrahigh purity and large sample mass for neutron scattering experiments are attainable through gradually-improved synthesis and thus, many materials may be ripe for magnetic field-based neutron scattering studies. These include DSMs, NLSMs, inversionsymmetry breaking WSMs such as the TaAs family and more intriguingly, magnetic WSMs such as Co 3 Sn 2 S 2 [78,79], Co 2 MnGa [80], CeAlGe [81,82] whereby fielddriven topological phase transitions may be uncovered. Preliminary works on inelastic neutron scattering serving as a platform for studying exotic phenomena in TNMs have only recently emerged. These include a theoretical overview on probing bulk excitations in type-I WSMs with unpolarized and polarized neutrons [83] and an experimental study of coupling between Dirac fermions and spin waves [84]. A non-exhaustive compilation of other neutron scattering experiments on topological semimetals (including candidates) up until now can be noted in Refs. [59,77,82,[84][85][86][87][88][89][90][91][92][93]. Although these earlier surveys are encouraging, a great deal of work remains to accentuate the potential usefulness of neutron scattering towards characterizing future TNMs sensitive to topological bands.

X. CONCLUSION
In summary, we provide a theoretical framework showing that TNMs can have large magnetic neutron scattering response even for non-magnetic materials. The calculation of the dynamical structure factor demonstrates that this large response is uncovered through a nonanalytical behavior of the measured spectra, seen as a divergent behavior of the second derivative with respect to wavevector Q or frequency ω. Notably, the high-intensity peaks of the second derivative are sensitive to the location of the Weyl nodes and the chemical potential of the Weyl nodes through the values of Q or ω at which the divergence takes place. The elevated sensitivity of this signal to these criteria may serve as a valuable method to extract information on parameters that are used to explicitly characterize TNMs, especially WSMs, using an experimental technique that has been hitherto underutilized for this purpose.
Notably, this study showcases that one can perform this type of measurement with triple-axis or time-of-flight neutron spectroscopies to determine the reciprocal-space location of the topological nodes, the magnitude of their Fermi velocities and the characteristic scattering time within these types of materials. In addition, upon application of parallel electric and magnetic fields, one can probe for a unique experimental signature of the chiral anomaly, a defining element of WSMs. Our study opens up new opportunities for the use of neutron scattering as an innovative experimental probe of the nontrivial topology of TNMs and complements the current repertoire of experimental techniques that serve in ongoing characterization efforts of these exotic quantum materials. We evaluate the magnetic scattering double differential cross-section term in Eq. (9). Let V represent the interaction potential operator of a neutron with the magnetic field. The double differential cross-section is [56] where |ξkσ = |ξ ⊗ |k ⊗ |σ is the tensor product of the material state |ξ , the neutron state |k , and the neutron spin state |σ . In Eq. (A1), m n is the neutron mass, V is the volume, is the reduced Planck constant, E is the energy, p is the probability of the neutron being in the spin up or spin down state and the subscripts i, f indicate the initial and final states, respectively. The interaction operator is given by where µ n (r) = −γµ N σ N is the magnetic moment operator of the neutron and B(r) is the magnetic field experienced by the neutron. In the definition of µ n (r), γ is the neutron gyromagnetic ratio, µ N is the nuclear magneton and σ N is the vector of neutron Pauli matrices. The magnetic field is where J is the current density operator. We proceed to calculate the matrix element k f σ f |V|k i σ i in Eq. (A1) using Eqs. (A2) and (A3): where Q ≡ k i − k f is the neutron momentum transfer and Q = |Q|. In the derivation above, we defined the Fourier transform of the current operator as To proceed further, we define the neutron magnetization operator in momentum space as where µ B = e /2m e c is the Bohr magneton in Gaussian units. The term in Eq. (A1) can be simplified as where we simplified the prefactor as where r e is the classical electron radius. In Eq. (A5), we label a = 1, 2, 3 for the spatial direction and the correlation is taken with respect to the initial material state.
We performed the summation over the neutron Pauli matrices with an assumption of non-polarized neutrons. Furthermore, we notice that the energy conservation factor in Eq. (A1) can be rewritten using the integral representation of the δ-function as where and M a (Q, t) = e iHt/ M a (Q, 0)e −iHt/ is the neutron magnetization operator in the Heisenberg picture. The neutron magnetization operator M defined in Eq. (A6) is written in a way that is perpendicular to Q. In the general case where localized spins may exist, we can rewrite the correlation in Eq. (A6) using the full neutron magnetization operator as where we suppress the time dependence in the notation. Using the correlation of the neutron magnetization operator defined above, we obtain the corresponding magnetic dynamical structure factor S jl (Q, ω) written in Eq. (4). In a state of thermal equilibrium, the double differential cross-section from Eq. (A1) can be written in terms of S jl (Q, ω) as where we used dE f = dω in the derivation. where j, l label the spatial direction. Taking the sum over initial and final states and using the same integral representation of the δ-function for energy conservation as in Eq. (A6), we obtain (B2) Now we use the fluctuation-dissipation theorem [94] in order to relate the dynamical structure factor S jl (Q, ω) with the dynamical response function Π im (Q, ω) as We obtain Eq. (6) by combining Eqs. (B1)-(B3).

(C3)
In the following, we use the shorthand notation v i,k = v i (k) corresponding to the velocity operator where v x,y = (1/m)(k x,y + Q x,y /2) and v z = v F σ y . For additional conciseness, we will write d i,k = d i (k), d k = d(k) and S (±) k = S (±) (k) when appropriate.
After rewriting in terms of S (±) (k) matrices, we can insert the Green's function into the polarization operator of Eq. (21) and perform the summation over Matsubara frequencies to obtain where n F (x) = (e βx + 1) −1 is the Fermi-Dirac distribution function. The imaginary part of Π im (Q, ω), which appears in the expression for the neutron dynamical structure factor written in Eq. (6), has the form where δ(x) is the delta function. We assume ω > 0 without loss of generality and take Q = Qẑ. Assuming zero Fermi level and zero temperature, the xx-component of the imaginary part of the response function is simplified to Eq. (22) because only one of the four contributions in Eq. (C4) remains nonzero as discussed in the main text.
To eliminate the δ-function, we need to take an integral over k z . In order to do so, we compute the derivative . By setting the argument of the δ-function to zero and solving for k z , we find the expressions .

(C5)
As we require that d k , d k+Q > 0, we can solve these inequalities using Eq. (C5) to obtain the condition v 2 F Q 2 − ω 2 < v 2 F Q 2 (k 2 ⊥ − p 2 0 ) 2 m 2 (ω 2 − v 2 F Q 2 ) which, in turn, implies that ω 2 − v 2 F Q 2 > 0. By enforcing this condition onto Eq. (22) and performing the integral over k z , we can simplify the expression to where the summation over k z,± refers to the two possibilities of sign in the expression for k z in Eq. (C5). The limits of integration of Eq. (C6) for k ⊥ are determined by the condition that the argument in the square roots of Eq. (C5) be strictly positive which is equivalent to