# Coulomb Corrections for Coherent Neutrino Nucleus Scattering

Ryan Plestid  Walter Burke Institute for Theoretical Physics, California
Institute of Technology, Pasadena, CA 91125

###### Abstract

In this work we consider sub-leading $O(Z^{2}\alpha^{3})$ corrections to
coherent elastic neutrino nucleus scattering (CEvNS). These corrections are
not negligible by power counting since nuclei with large coherent cross
sections have sizeable nuclear charges e.g. $Z\alpha\sim 0.4$. We find that
the corrections are much smaller than naive power counting in $Z\alpha$ would
suggest and that Coulomb corrections do not substantially alter predictions
for CEvNS in the Standard Model. We comment on similarities to older
literature on mesonic and muonic atoms.

††preprint: CALT-TH-2023-008

##  I Introduction

Coherent elastic neutrino nucleus scattering Stodolsky (1966); Freedman
(1974); Kopeliovich and Frankfurt (1974) (CEvNS) offers a clean probe of the
weak neutral current. It is a leading tool for constraining neutrinophillic
models of new physics deNiverville  _et al._ (2015); Papoulias  _et al._
(2019); Billard  _et al._ (2018); Altmannshofer  _et al._ (2019); Abdullah
_et al._ (2018); Blanco  _et al._ (2020); Farzan  _et al._ (2018); Miranda
_et al._ (2020), and it furthermore supplies valuable information on nuclear
structure Amanik and McLaughlin (2007); Huang and Chen (2019); Coloma  _et
al._ (2020); Van Dessel  _et al._ (2020); Abdullah  _et al._ (2022). CEvNS was
recently measured for the first time Akimov  _et al._ (2017) and there now
exist dozens of experimental collaborations that will attempt to measure a
CEvNS signal in the near future Abdullah  _et al._ (2022). The COHERENT
collaboration, in particular, has aspirations of sub-percent level precision
(see e.g. Fig. 19 of Akimov  _et al._ (2022)). Similar projections of percent
level measurements have been proposed for reactor neutrinos with low-threshold
detectors Fernandez-Moroni  _et al._ (2021).

Motivated by recent interest in precision measurements of the CEvNS cross
section, the authors of Tomalak  _et al._ (2021) performed a study of
radiative corrections to CEvNS at $O(\alpha)$ (see also earlier work Sehgal
(1985); Botella  _et al._ (1987) along these lines). The calculation was
performed in a low energy effective field theory (EFT) in which any large “RG-
logs” had already been resummed into Wilson coefficients Hill and Tomalak
(2020). Accounting for hadronic, nuclear, and perturbative uncertainties, the
authors of Tomalak  _et al._ (2021) concluded that CEvNS represents a bonafide
precision observable within the SM. In this work I consider Coulomb
corrections to the charged lepton loop that mediates photon exchange with the
nucleus; the diagram of interest is depicted in Fig. 1. This constitutes the
only relevant radiative corrections that was not considered previously; other
corrections exist but are suppressed by higher powers of $\alpha$ with no
coherent enhancements. Formally, Fig. 1 represents a relative
$O(Z^{2}\alpha^{2})$ correction to the results computed in Tomalak  _et al._
(2021) (which are themselves $O(Z\alpha)$). For heavy nuclei this relative
correction could be sizeable. As we will now discuss, however, the correction
is actually much smaller than naive power counting would suggest and can be
safely neglected for practical purposes.

The rest of the paper is organized as follows: Section II discusses how
Coulomb corrections enter in CEvNS and reviews existing results on vacuum
polarization in an electrostatic field. Next in Section III, I compute Coulomb
corrections from lepton loops for CEvNS in the point-like limit making use of
existing work from the muonic atom literature; the effect is found to be
small. In Section IV I discuss finite size corrections and hadronic vacuum
polarization and argue these effects, while interesting theoretically, are not
relevant for a percent-level determination of the cross section. Finally in
Section V I summarize, conclude, and comment on the position of CEvNS as a
precision SM observable.

Figure 1: Graph that contributes to CEvNS considered in this work. The square
denotes a four-Fermi contact operator (neutrinos not drawn) and the double
line the nucleus. In the infinite mass limit the nucleus behaves as a static
potential of charge $Z$. The loop is any dynamical charged particle in the
effective theory e.g. muons, electrons, or pions. In practice only electrons
are relevant phenomenologically.

##  II Coulomb corrections

CEvNS is a particularly simple system from the perspective of radiative
corrections due to the absence of charge particles on the external legs.
Nuclei carry macroscopic charge, however their coupling to real photons is
mass-suppressed111This is most easily seen by using a heavy particle
Lagrangian and working in Coulomb gauge. and their only influence at leading
order in $1/M$ is via the exchange of Coulomb photons. Short-distance
radiative corrections can be incorporated via standard EFT techniques. The
only residual contribution that must then be considered involve loops of light
particles.

In what follows we discuss radiative corrections at $O(Z^{n}\alpha^{n+1})$.
This restriction allows us to express the amplitude as

$\mathcal{M}=\mathcal{M}_{0}+\frac{\alpha}{\pi}\sum_{n}(Z\alpha)^{2n}\mathcal{M}_{1}^{(n)}~{},$
(1)

Only even powers of $Z\alpha$ appear by Furry’s theorem Furry (1937). The
tree-level amplitude $\mathcal{M}_{0}$ is proportional to the weak nuclear
form factor, $F_{W}(Q^{2})$, whereas all $O(\alpha)$ corrections in
$\mathcal{M}_{1}^{(n)}$ depend only on the charge form factor, $F_{\rm
ch}(Q^{2})$.

At one-loop the amplitude is modified in a simple way that depends only on the
QED vacuum polarization Tomalak  _et al._ (2021)

$\mathcal{M}_{1}^{(0)}=\frac{\alpha}{\pi}F_{\rm
ch}(Q^{2})(\delta^{\nu_{\ell}}+\delta^{\rm QCD})~{},$ (2)

where $\delta^{\nu_{\ell}}$ is the sum of vacuum polarizations
$\Pi(Q^{2},m_{\ell})$, weighted by flavor dependent Wilson coefficients. The
correction $\delta^{\rm QCD}$ is likewise related to hadronic vacuum
polarization, and a charge isospin correlator Jegerlehner (1986); Tomalak and
Hill (2020); Tomalak  _et al._ (2021).

The first non-trivial effect from the nuclear Coulomb field appears at three-
loop order as shown in Fig. 1. In what follows we will consider the point-like
limit, $F_{\rm ch}(Q^{2})\rightarrow 1$, for which analytic results are known
Wichmann and Kroll (1956)

$\mathcal{M}_{1}^{(1)}=\frac{Z}{Q_{w}}\mathcal{M}_{0}\sum c_{V}^{(f)}\Pi_{\rm
WK}(Q^{2},m_{f})~{},$ (3)

where $Q_{w}$ is the weak nuclear charge, $c_{V}^{(f)}$ is a short-hand for
the vectorial Wilson coefficients that appear in Eqs. (2.4) and (2.5) of
Tomalak  _et al._ (2021), and the subscript “WK” stands for Wichmann-Kroll
Wichmann and Kroll (1956). The Wichmann-Kroll contribution is typically
studied in coordinate space and can be expressed as a weighted integral over
Yukawa-like contributions Blomqvist (1972). Fourier transforming to momentum
space, a similar representation may be found for the vacuum polarization

$\begin{split}\Pi_{\rm
WK}(Q^{2},m_{f})=\frac{1}{2}\int_{0}^{\infty}&\differential\zeta\frac{Q^{2}}{Q^{2}+4\zeta^{2}m_{f}^{2}}\frac{1}{\zeta^{4}}\\\
&\quantity[F(\zeta)-\frac{\pi^{2}}{12}\sqrt{\zeta^{2}-1}\Theta(\zeta-1)]~{}.\end{split}$
(4)

The function $F(\zeta)$ has relatively lengthy integral representation which
can be found in the literature (see e.g. pages 99-100 of Blomqvist
(1972)).222It is useful to note that $F(\zeta)\sim O(\zeta^{7})$ for small
$\zeta$ such that Eq. 4 is a well defined convergent integral.

In the limit of a massless fermion, or equivalently when $Q^{2}\gg
m_{\ell}^{2}$, results are known to all orders in $Z\alpha$ Brown  _et al._
(1974a, b). Structure dependent corrections have also been investigated in for
$QR\lesssim 1$ in this limit Brown  _et al._ (1975). In what follows we apply
these well known results to CEvNS and quantify the size of higher order
Coulomb corrections for this process from closed lepton loops.

##  III Results in the point-like limit

We have numerically computed the Wichmann-Kroll vacuum polarization using the
integral representation given in Eq. 4. Corrections to CEvNS turn out to be
small. For the sake of simplicity, we quantify this with the lepton-flavor
asymmetry defined in Tomalak  _et al._ (2021). Defining the total vacuum
polarization via

$\Pi=\Pi_{1}+(Z\alpha)^{2}\Pi_{\rm WK}+\ldots,$ (5)

with $\Pi_{1}$ the leading-in-$\alpha$ contribution, the flavor asymmetry in
CEvNS is then given by

$\begin{split}\mathcal{A}_{\ell,\ell^{\prime}}&=\frac{\differential\sigma_{\nu_{\ell}}-\differential\sigma_{\nu_{\ell^{\prime}}}}{\differential\sigma_{\nu_{\ell}}}\\\
&=\frac{4\alpha}{\pi}\frac{Z}{Q_{w}}\quantity[\Pi(Q^{2},m_{\ell})-\Pi(Q^{2},m_{\ell^{\prime}})]~{}.\end{split}$
(6)

This observable has the advantage of being scheme independent with no residual
scale dependence. The Wichmann-Kroll contribution to $\Pi(Q^{2},m_{\ell})$
shifts this ratio at $O(Z^{2}\alpha^{2})$. We express this dependence as

$\mathcal{A}_{\ell,\ell^{\prime}}=\frac{4\alpha}{\pi}\frac{Z}{Q_{w}}\quantity(\mathcal{A}_{\ell,\ell^{\prime}}^{(0)}+(Z\alpha)^{2}\mathcal{A}_{\ell,\ell^{\prime}}^{(1)}+\ldots)$
(7)

The numerical size of these effects is shown in Table 1 for three different
values of $Q^{2}$. One can clearly see that coefficient of three-loop
$O(Z^{2}\alpha^{3})$ correction is much smaller than that of the vacuum
polarization. Similar conclusions apply to the total cross section and indeed
any other observable.

Table 1: The CEvNS flavor asymmetry for $Q^{2}=(1~{}{\rm MeV})^{2}$, $(10~{}{\rm MeV})^{2}$, and $Q^{2}=(30~{}{\rm MeV})^{2}$ in the point-like limit. The result is essentially $Q$-independent above a few MeV. In all cases the Coulomb corrections are very small (even when $Z\alpha$ is not). The $\mu-\tau$ asymmetries are completely negligible.  $Q$ | Asymmetry | $\mathcal{A}_{\ell,\ell^{\prime}}^{(0)}$ | $\mathcal{A}_{\ell,\ell^{\prime}}^{(1)}$
---|---|---|---
|  $\ell=e$ $\ell^{\prime}=\mu$ | $3.36$ | $0.02$
$1~{}{\rm MeV}~{}$ |  $\ell=e$ $\ell^{\prime}=\tau$ | $5.24$ | $0.02$
|  $\ell=\mu$ $\ell^{\prime}=\tau$ | $1.88$ | —
|  $\ell=e$ $\ell^{\prime}=\mu$ | $2.12$ | $0.03$
$10~{}{\rm MeV}~{}$ |  $\ell=e$ $\ell^{\prime}=\tau$ | $4.00$ | $0.03$
|  $\ell=\mu$ $\ell^{\prime}=\tau$ | $1.88$ | —
|  $\ell=e$ $\ell^{\prime}=\mu$ | $1.40$ | $0.03$
$30~{}{\rm MeV}~{}$ |  $\ell=e$ $\ell^{\prime}=\tau$ | $3.28$ | $0.03$
|  $\ell=\mu$ $\ell^{\prime}=\tau$ | $1.88$ | —

Let us now discuss why the corrections are small for heavy leptons, ($\mu$ and
$\tau$) v.s. light dynamical ones (the electron). All contributions involving
3 or more photon attachments to the lepton loop are UV-convergent,333This is
not obvious by superficial power counting but is guaranteed by gauge
invariance Karplus and Neuman (1951); Brown  _et al._ (1974b); Liang and
Czarnecki (2012). and must decouple as $1/m_{\ell}^{2}$ in the limit that
$m_{\ell}\rightarrow\infty$. This is in contrast to the vacuum polarization
where a residual logarithmic dependence, correcting for the running value of
$\alpha$, persists. This then guarantees that corrections are parametrically
suppressed by $(Z\alpha)^{2}Q^{2}/m_{\ell}^{2}$. For electrons, $Q^{2}\gg
m_{e}^{2}$ and so this argument does not apply, however the short-distance
coefficient of the Wichmann-Kroll potential is small being given by

$\lim_{Q\rightarrow\infty}\Pi_{\rm
WK}(Q/m_{\ell})=\frac{1}{2}\quantity[\frac{\pi^{2}}{6}-\frac{2}{3}\zeta(3)-\frac{7}{9}]\approx
0.03~{}.$ (8)

This suppression has been shown to persist to all orders in $Z\alpha$ for the
hierarchy $m_{e}\ll Q\ll 1/R$ where $R$ is the nuclear radius; we discuss
finite size effects in the next section. A cursory examination of Eq. 4 and
the observation that $0\leq Q^{2}/(Q^{2}+4m_{\ell}^{2}\zeta^{2})\leq 1$
guarantees that $\Pi_{\rm WK}(Q^{2})$ is a monotonically increasing function
of $Q^{2}$; this is shown explicitly in Fig. 2. Therefore, since corrections
are small in the limit that $Q^{2}\gg m_{e}^{2}$ they are small for all values
of $Q^{2}$.

Figure 2: Numerical calculation of the Wichmann-Kroll vacuum polarization as a
function of $Q/m_{\ell}$. One clearly sees the Coulomb corrections are always
much smaller than naive power-counting in $Z\alpha$ would suggest.

##  IV Additional corrections

In the previous section we have presented results for higher order corrections
to the vacuum polarization in the point-like limit. We focused on leptonic
contributions, however in the full CEvNS cross section hadronic loops also
contribute. One may ask how hadronic contributions to the vacuum polarization
and finite size corrections may alter our results. In this section we will
argue these effects are never large, nor are finite size corrections from the
charge distribution of the nucleus.

Hadronic contributions will always be small on dimensional grounds. For the
same reason that heavy leptons decouple, hadrons must also. Since the lightest
hadron is $m_{\pi}$ the effect of hadronic loops will always be suppressed by
the same factor of $Q^{2}/m_{\pi}^{2}$. This ratio is typically smaller than
$0.1$ for CEvNS kinematics. It is possible the Coulomb corrections from pion
loops may be larger than the leptonic contribution, however they will never be
“large” in so far as $(Z\alpha)^{2}Q^{2}/m_{\pi}^{2}\ll 1$.

Our discussion so far has focused on the point-like limit. Finite size effects
have been previously studied in the context of mesonic and muonic atoms where
a coordinate space representation for the potential was derived Brown  _et
al._ (1975), valid to all orders in $\alpha(Z\alpha)^{2n+1}$, for $r\gg R$;
finite size effects were found to be small as one would expect. For neutrino
energies above $30~{}{\rm MeV}$ or so, a sizeable fraction of the CEvNS cross
section comes from regions where $QR\sim 1$ due to the phase-space enhancement
of large-$Q$ regions. In these regions of phase space the point-like
calculations presented above do not strictly speaking apply, but should
provide a reasonable estimate of the radiative corrections. Since they are so
small we do not attempt to compute the three-loop integral or to construct
modified Greens functions for an extended nuclear charge distribution. In
addition to anticipating that the corrections will never be larger than in the
point-like limit, we note that in the regime where $QR\sim 1$ uncertainties
due to the tree-level nuclear form factor dominate the error budget
contributing an uncertainty at the level of a few percent (see Table 2 of
Tomalak  _et al._ (2021)).

Finally, one may wonder if the cancellation in Eq. 8 is accidental such that
higher order terms may be important. The analysis of Brown  _et al._ (1974b)
shows this not to be the case; all higher order corrections are similarly
suppressed by a factor of roughly $1/100$.

##  V Conclusions & Outlook

In this work the leading Coulomb corrections to CEvNS, directly related to the
Wichmann-Kroll potential, have been calculated in the limit of a point-like
nucleus. The corrections are found to be small. The influence of finite size
corrections, mixed hadronic correlators, and other SM effects not explicitly
computed have been estimated and will not affect the CEvNS cross section at a
level relevant for experiments (i.e. they are sub-permille).

In a broader context, this work demonstrates that Coulomb corrections are
under control, a necessary requirement for a precision CEvNS program. This
compliments other recent additions to the literature such as an estimate of
inelastic invisible contributions to the cross section Dutta  _et al._ (2022),
and percent-level determinations of neutrino fluxes Rapp (2019); Tomalak
(2022). Further work on hadronic contributions, via the mixed correlator
$\Pi_{3\gamma}$ discussed in Tomalak  _et al._ (2021), and a reduction of
uncertainties in the neutron and proton distributions inside nuclei offer a
path towards a sub percent prediction of the CEvNS cross section.

## Acknowledgements

I would like to thank Oleksandr Tomalak, Richard Hill, and especially Robert
Szafron for helpful discussions, and both Robert Szafron and Oleksandr Tomalak
for reading this manuscript before completion. This work was instigated at the
Mainz Institute for Theoretical Physics, and I thank them for their
hospitality. This work is funded by the Neutrino Theory Network Program Grant
under Award Number DEAC02-07CHI11359 and the US DOE under Award Number DE-
SC0020250. This work is supported in part by the U.S. Department of Energy,
Office of Science, Office of High Energy Physics, under Award Number DE-
SC0011632 and by the Walter Burke Institute for Theoretical Physics.

## References

  * Stodolsky (1966) Leo Stodolsky, “Application of Nuclear Coherence Properties to Elementary-Particle Reactions,” Phys. Rev. 144, 1145–1153 (1966).
  * Freedman (1974) Daniel Z. Freedman, “Coherent Neutrino Nucleus Scattering as a Probe of the Weak Neutral Current,” Phys. Rev. D 9, 1389–1392 (1974).
  * Kopeliovich and Frankfurt (1974) V.B. Kopeliovich and L.L. Frankfurt, “Isotopic and chiral structure of neutral current,” JETP Lett. 19, 145–147 (1974).
  * deNiverville  _et al._ (2015) Patrick deNiverville, Maxim Pospelov, and Adam Ritz, “Light new physics in coherent neutrino-nucleus scattering experiments,” Phys. Rev. D 92, 095005 (2015), arXiv:1505.07805 [hep-ph] .
  * Papoulias  _et al._ (2019) D. K. Papoulias, T. S. Kosmas, and Y. Kuno, “Recent probes of standard and non-standard neutrino physics with nuclei,” Front. in Phys. 7, 191 (2019), arXiv:1911.00916 [hep-ph] .
  * Billard  _et al._ (2018) Julien Billard, Joseph Johnston, and Bradley J. Kavanagh, “Prospects for exploring New Physics in Coherent Elastic Neutrino-Nucleus Scattering,” JCAP 11, 016 (2018), arXiv:1805.01798 [hep-ph] .
  * Altmannshofer  _et al._ (2019) Wolfgang Altmannshofer, Michele Tammaro, and Jure Zupan, “Non-standard neutrino interactions and low energy experiments,” JHEP 09, 083 (2019), [Erratum: JHEP 11, 113 (2021)], arXiv:1812.02778 [hep-ph] .
  * Abdullah  _et al._ (2018) Mohammad Abdullah, James B. Dent, Bhaskar Dutta, Gordon L. Kane, Shu Liao, and Louis E. Strigari, “Coherent elastic neutrino nucleus scattering as a probe of a Z’ through kinetic and mass mixing effects,” Phys. Rev. D 98, 015005 (2018), arXiv:1803.01224 [hep-ph] .
  * Blanco  _et al._ (2020) Carlos Blanco, Dan Hooper, and Pedro Machado, “Constraining Sterile Neutrino Interpretations of the LSND and MiniBooNE Anomalies with Coherent Neutrino Scattering Experiments,” Phys. Rev. D 101, 075051 (2020), arXiv:1901.08094 [hep-ph] .
  * Farzan  _et al._ (2018) Yasaman Farzan, Manfred Lindner, Werner Rodejohann, and Xun-Jie Xu, “Probing neutrino coupling to a light scalar with coherent neutrino scattering,” JHEP 05, 066 (2018), arXiv:1802.05171 [hep-ph] .
  * Miranda  _et al._ (2020) O. G. Miranda, D. K. Papoulias, O. Sanders, M. Tórtola, and J. W. F. Valle, “Future CEvNS experiments as probes of lepton unitarity and light-sterile neutrinos,” Phys. Rev. D 102, 113014 (2020), arXiv:2008.02759 [hep-ph] .
  * Amanik and McLaughlin (2007) Philip S. Amanik and Gail C. McLaughlin, “Neutron Form Factor from Neutrino-Nucleus Coherent Elastic Scattering,” (2007), arXiv:0707.4191 [hep-ph] .
  * Huang and Chen (2019) Xu-Run Huang and Lie-Wen Chen, “Neutron Skin in CsI and Low-Energy Effective Weak Mixing Angle from COHERENT Data,” Phys. Rev. D 100, 071301 (2019), arXiv:1902.07625 [hep-ph] .
  * Coloma  _et al._ (2020) Pilar Coloma, Ivan Esteban, M. C. Gonzalez-Garcia, and Javier Menendez, “Determining the nuclear neutron distribution from Coherent Elastic neutrino-Nucleus Scattering: current results and future prospects,” JHEP 08, 030 (2020), arXiv:2006.08624 [hep-ph] .
  * Van Dessel  _et al._ (2020) N. Van Dessel, V. Pandey, H. Ray, and N. Jachowicz, “Nuclear Structure Physics in Coherent Elastic Neutrino-Nucleus Scattering,” (2020), arXiv:2007.03658 [nucl-th] .
  * Abdullah  _et al._ (2022) M. Abdullah _et al._ , “Coherent elastic neutrino-nucleus scattering: Terrestrial and astrophysical applications,” (2022), arXiv:2203.07361 [hep-ph] .
  * Akimov  _et al._ (2017) D. Akimov _et al._ (COHERENT), “Observation of Coherent Elastic Neutrino-Nucleus Scattering,” Science 357, 1123–1126 (2017), arXiv:1708.01294 [nucl-ex] .
  * Akimov  _et al._ (2022) D. Akimov _et al._ , “The COHERENT Experimental Program,” in  _Snowmass 2021_ (2022) arXiv:2204.04575 [hep-ex] .
  * Fernandez-Moroni  _et al._ (2021) Guillermo Fernandez-Moroni, Pedro A. N. Machado, Ivan Martinez-Soler, Yuber F. Perez-Gonzalez, Dario Rodrigues, and Salvador Rosauro-Alcaraz, “The physics potential of a reactor neutrino experiment with Skipper CCDs: Measuring the weak mixing angle,” JHEP 03, 186 (2021), arXiv:2009.10741 [hep-ph] .
  * Tomalak  _et al._ (2021) Oleksandr Tomalak, Pedro Machado, Vishvas Pandey, and Ryan Plestid, “Flavor-dependent radiative corrections in coherent elastic neutrino-nucleus scattering,” JHEP 02, 097 (2021), arXiv:2011.05960 [hep-ph] .
  * Sehgal (1985) L. M. Sehgal, “Differences in the Coherent Interactions of $\nu_{e}$, $\nu_{\mu}$ and $\nu_{\tau}$,” Phys. Lett. B 162, 370–372 (1985).
  * Botella  _et al._ (1987) F. J. Botella, C. S. Lim, and W. J. Marciano, “Radiative Corrections to Neutrino Indices of Refraction,” Phys. Rev. D 35, 896 (1987).
  * Hill and Tomalak (2020) Richard J. Hill and Oleksandr Tomalak, “On the effective theory of neutrino-electron and neutrino-quark interactions,” Phys. Lett. B 805, 135466 (2020), arXiv:1911.01493 [hep-ph] .
  * Furry (1937) W. H. Furry, “A Symmetry Theorem in the Positron Theory,” Phys. Rev. 51, 125–129 (1937).
  * Jegerlehner (1986) F. Jegerlehner, “Hadronic Contributions to Electroweak Parameter Shifts: A Detailed Analysis,” Z. Phys. C 32, 195 (1986).
  * Tomalak and Hill (2020) Oleksandr Tomalak and Richard J. Hill, “Theory of elastic neutrino-electron scattering,” Phys. Rev. D 101, 033006 (2020), arXiv:1907.03379 [hep-ph] .
  * Wichmann and Kroll (1956) Eyvind H. Wichmann and Norman M. Kroll, “Vacuum Polarization in a Strong Coulomb Field,” Phys. Rev. 101, 843–859 (1956).
  * Blomqvist (1972) J. Blomqvist, “Vacuum polarization in exotic atoms,” Nucl. Phys. B 48, 95–103 (1972).
  * Brown  _et al._ (1974a) Lowell S. Brown, Robert N. Cahn, and Larry D. McLerran, “Analytic Calculation to All Orders in Z-alpha of Nuclear Size Effects in Vacuum Polarization,” Phys. Rev. Lett. 32, 562 (1974a).
  * Brown  _et al._ (1974b) Lowell S. Brown, Robert N. Cahn, and Larry D. McLerran, “Analytical Calculation to All Orders in Z alpha of Vacuum Polarization at Short Distances,” Phys. Rev. Lett. 33, 1591 (1974b).
  * Brown  _et al._ (1975) Lowell S. Brown, Robert N. Cahn, and Larry D. McLerran, “Vacuum Polarization in a Strong Coulomb Field. 3. Nuclear Size Effects,” Phys. Rev. D 12, 609 (1975).
  * Karplus and Neuman (1951) Robert Karplus and Maurice Neuman, “The scattering of light by light,” Phys. Rev. 83, 776–784 (1951).
  * Liang and Czarnecki (2012) Yi Liang and Andrzej Czarnecki, “Photon-photon scattering: A Tutorial,” Can. J. Phys. 90, 11–26 (2012), arXiv:1111.6126 [hep-ph] .
  * Dutta  _et al._ (2022) Bhaskar Dutta, Wei-Chih Huang, Jayden L. Newstead, and Vishvas Pandey, “Inelastic nuclear scattering from neutrinos and dark matter,” Phys. Rev. D 106, 113006 (2022), arXiv:2206.08590 [hep-ph] .
  * Rapp (2019) Rebecca Rapp (COHERENT), “SNS Neutrino Flux,” Workshop on Fundamental Physics at the Second Target Station (2019), Link.
  * Tomalak (2022) Oleksandr Tomalak, “Radiative (anti)neutrino energy spectra from muon, pion, and kaon decays,” Phys. Lett. B 829, 137108 (2022), arXiv:2112.12395 [hep-ph] .


# Large dust fractions can prevent the propagation of soundwaves

Timothée David–Cléris1, Guillaume Laibe1,2\-
1Univ Lyon, Univ Lyon1, Ens de Lyon, CNRS, Centre de Recherche Astrophysique
de Lyon UMR5574, F-69230, Saint-Genis,-Laval, France.
2Institut Universitaire de France
timothee.david–cleris<EMAIL_ADDRESS>
###### Abstract

Dust plays a central role in several astrophysical processes. Hence the need
of dust/gas numerical solutions, and analytical problems to benchmark them. In
the seminal dustywave problem, we discover a regime where sound waves can not
propagate through the mixture above a large critical dust fraction. We
characterise this regime analytically, making it of use for testing accuracy
of numerical solvers at large dust fractions.

###### keywords:

(ISM:) dust, extinction — methods: analytical — protoplanetary discs

††pagerange: Large dust fractions can prevent the propagation of
soundwaves–5††pubyear: 2020

##  1 Introduction

Quantitative study of dust is of prime importance in astrophysics. Numerical
simulations are used for determining the 3D evolution of dust/gas systems.
Providing accurate tests to benchmark these numerical codes is therefore
critical to ensure reliability of the results. Astrophysical dust is usually
modelled by a pressureless continuum that exchanges momentum with the gas
through a drag force (Saffman, 1962; Baines et al., 1965; Clair et al., 1970;
Marble, 1970). Several analytical problems involving advection, waves, shocks,
settling or dust/gas instabilities have been used to benchmark dust/gas codes
(e.g. Benítez-Llambay et al. 2019; Stoyanovskaya et al. 2020 and references
therein for recent discussions).

The dustywave problem consists of the propagation in 1D of a sound wave in
such a mixture (Ahuja 1973; Gumerov et al. 1988; Laibe & Price 2011, 2016 –
see Sect. 2). dustywave is one of the most widely used benchmark, since it
associates dust/gas drag and gas compressibility, both in Lagrangian (e.g.
Laibe & Price 2012, 2014a; Lorén-Aguilar & Bate 2014; Booth et al. 2015; Price
& Laibe 2015; Stoyanovskaya et al. 2018; Mentiplay et al. 2020) or in Eulerian
methods (e.g. Porth et al. 2014; Yang & Johansen 2016; Hubber et al. 2018;
McKinnon et al. 2018; Riols & Lesur 2018; Lebreuilly et al. 2019; Moseley et
al. 2019). In attempting to benchmark a numerical code against this test, we
figured out the existence of a regime at large dust-to-gas ratios where waves
cannot propagate. After having recalled briefly the main properties the
dustywave problem in Sect. 2, we show the existence of this regime, derive
analytic values for the corresponding boundaries, provide physical
explanations and numerical tests in Sect. 3.

##  2 Dispersion relation

The equations of evolution for a 1D astrophysical dusty mixture are

$\displaystyle\partial_{t}\rho_{\rm g}+{v}_{\rm g}\partial_{x}\rho_{\rm g}$
$\displaystyle=-\rho_{\rm g}\,\partial_{x}{v}_{\rm g},$ (1)
$\displaystyle\partial_{t}\rho_{\rm d}+{v}_{\rm d}\partial_{x}\rho_{\rm d}$
$\displaystyle=-\rho_{\rm d}\,\partial_{x}{v}_{\rm d},$ (2)
$\displaystyle\partial_{t}{v}_{\rm g}+{v}_{\rm g}\partial_{x}{v}_{\rm g}$
$\displaystyle=+\frac{K}{\rho_{\rm g}}\left({v}_{\rm d}-{v}_{\rm
g}\right)-\frac{\partial_{x}P}{\rho_{\rm g}},$ (3)
$\displaystyle\partial_{t}{v}_{\rm d}+{v}_{\rm d}\partial_{x}{v}_{\rm d}$
$\displaystyle=-\frac{K}{\rho_{\rm d}}\left({v}_{\rm d}-{v}_{\rm g}\right),$
(4)

where $\rm g$ and $\rm d$ stand for gas and dust respectively (e.g. Garaud et
al. 2004) and $K$ denotes the drag coefficient. Assuming isothermal gas
$P=c_{\rm s}^{2}\rho_{\rm g}$, we expand linearly Eqs. 1 – 4 under the generic
form $a=a_{0}+\delta a$, with $v_{\mathrm{d}0}=v_{\mathrm{g}0}=0$. One obtains

$\displaystyle\partial_{t}\delta\rho_{\rm g}$
$\displaystyle=-\rho_{\mathrm{g},0}\,\partial_{x}\delta v_{\rm g},$ (5)
$\displaystyle\partial_{t}\delta\rho_{\rm d}$
$\displaystyle=-\rho_{\mathrm{d},0}\,\partial_{x}\delta v_{\rm d},$ (6)
$\displaystyle\partial_{t}\delta v_{\rm g}$
$\displaystyle=+\frac{K}{\rho_{\mathrm{g},0}}\left(\delta v_{\rm d}-\delta
v_{\rm g}\right)-c_{\rm s}^{2}\frac{\partial_{x}\delta\rho_{\rm
g}}{\rho_{\mathrm{g},0}},$ (7) $\displaystyle\partial_{t}\delta v_{\rm d}$
$\displaystyle=-\frac{K}{\rho_{\mathrm{d},0}}\left(\delta v_{\rm d}-\delta
v_{\rm g}\right).$ (8)

We decompose the perturbation on Fourier space under the form $\delta
a=\tilde{a}e^{i(kx-\omega t)}$ for each perturbed field, giving the condition

$\begin{vmatrix}-i\omega&0&\frac{ic_{\rm
s}k\rho_{\mathrm{g},0}}{\rho_{\mathrm{g},0}+\rho_{\mathrm{d},0}}&0\\\
0&-i\omega&0&\frac{ic_{\rm
s}k\rho_{\mathrm{d},0}}{\rho_{\mathrm{g},0}+\rho_{\mathrm{d},0}}\\\
\frac{ic_{\rm
s}k\left(\rho_{\mathrm{g},0}+\rho_{\mathrm{d},0}\right)}{\rho_{\mathrm{g},0}}&0&-i\omega+\frac{1}{t_{\mathrm{g}}}&-\frac{1}{t_{\mathrm{g}}}\\\
0&0&-\frac{1}{t_{\mathrm{d}}}&-i\omega+\frac{1}{t_{\mathrm{d}}}\end{vmatrix}=0.$
(9)

Figure 1: Imaginary parts of the roots $\omega_{\pm}$ and $\omega_{0}$ for
different values of the dust fraction $\epsilon$ revealing a bifurcation at
$\epsilon\geq\epsilon_{\rm c}=8/9$.

We obtain the following dispersion relation

$\omega^{4}+\frac{i}{t_{\rm s}}\omega^{3}-c_{\rm
s}^{2}k^{2}\omega^{2}-\frac{i}{t_{\rm s}}c_{\rm
s}^{2}k^{2}(1-\epsilon)\omega=0,$ (10)

where the barycentric stopping time is $t_{\rm
s}\equiv\frac{\rho_{\mathrm{g},0}\rho_{\mathrm{d},0}}{K\left(\rho_{\mathrm{g},0}+\rho_{\mathrm{d},0}\right)}$
and
$\epsilon\equiv\frac{\rho_{\mathrm{d},0}}{\rho_{\mathrm{g},0}+\rho_{\mathrm{d},0}}$
is the total dust fraction. Rescaling time and space by $t_{\rm s}$ and
$c_{\rm s}t_{\rm s}$ respectively gives in the dimensionless form

$\omega^{4}+i\omega^{3}-k^{2}\omega^{2}-ik^{2}(1-\epsilon)\omega=0,$ (11)

where we preserved the notations $\omega$ and $k$ for further readability. We
disregard the solution $\omega_{\rm null}=0$ on the null space. On the column
space, Eq. 11 reduces to

$\omega^{3}+i\omega^{2}-\omega k^{2}-ik^{2}(1-\epsilon)=0,$ (12)

which can alternatively be written under the convenient form.

$\omega^{2}-k^{2}+\frac{i}{\omega}\left(\omega^{2}-k^{2}(1-\epsilon)\right)=0.$
(13)

The change of variable $\omega=iy$ gives a cubic with real positive
coefficients

$y^{3}+y^{2}+yk^{2}+k^{2}(1-\epsilon)=0.$ (14)

When Eq. 14 admits two complex conjugated roots and one real root, the latter
is negative since $k^{2}(1-\epsilon)>0$. When Eq. 14 admits three real roots,
Descarte’s rule of sign shows that they are all negative. Since
$\Im\left(\omega\right)=\Re\left(y\right)$, all modes of the dustywave problem
are damped. This result can alternatively be found from the argument principle
(Debras et al., 2020). Let split $\omega$ in its real and imaginary part by
setting $\omega\equiv r+is$. One obtains

$\displaystyle 0$ $\displaystyle=r\left[r^{2}-(3s^{2}+2s+k^{2})\right],$ (15)
$\displaystyle 0$
$\displaystyle=s^{3}+s^{2}+s(k^{2}-3r^{2})+(k^{2}(1-\epsilon)-r^{2}).$ (16)

Eq. 15 shows that the three expected modes decompose as follow:

$\displaystyle r$ $\displaystyle=0,$ (17) $\displaystyle 0$
$\displaystyle=s^{3}+s^{2}+sk^{2}+k^{2}(1-\epsilon),$ (18)

and

$\displaystyle r^{2}$ $\displaystyle=3s^{2}+2s+k^{2},$ (19) $\displaystyle 0$
$\displaystyle=s^{3}+s^{2}+s\frac{1}{4}\left(k^{2}+1\right)+k^{2}\frac{\epsilon}{8}.$
(20)

We note that the symmetry $r\to-r$ implies three solutions
$\omega_{\pm}=\pm|r|+is$ and $\omega_{0}=is$. Eqs. 17 – 18 give solutions that
are always purely damped. Eqs. 19 – 20 give solutions that are contra-
propagative and damped.

##  3 Absence of sound propagation

###  3.1 Analysis

Figure 2: Real part of the root $\omega_{+}$ for different values of the dust
fraction $\epsilon$ (solid lines) superimposed to the asymptotic regimes
(dashed lines). A forbidden band develops for $\epsilon\geq\epsilon_{\rm
c}=8/9$.

Fig. 2 shows that sound waves behave in a specific manner at large dust
fractions. Above a critical value $\epsilon_{\rm c}$, an interval of values
for $k$ where waves can not propagate develops (a so-called forbidden band).
We relate the existence of the bifurcation at $\epsilon_{\rm c}$ to the fact
that above this value, all three solutions come solely from Eq. 18. From the
discriminant of the cubic Eq. 12, one gets the condition
$\epsilon\geq\epsilon_{\rm c}=8/9$ and $k_{\rm min}\leq k\leq k_{\rm max}$
with

$\begin{cases}k_{\min}&\equiv\displaystyle\sqrt{\left(1-\frac{9\epsilon}{8}\right)\left(3\epsilon-\frac{4}{3}+\sqrt{\epsilon}\sqrt{9\epsilon-8}\right)+\frac{1}{3}},\\\
k_{\max}&\equiv\displaystyle\sqrt{\left(1-\frac{9\epsilon}{8}\right)\left(3\epsilon-\frac{4}{3}-\sqrt{\epsilon}\sqrt{9\epsilon-8}\right)+\frac{1}{3}},\\\
\end{cases}$ (21)

The centre of the band $k_{\rm c}$ can be estimated from the relation

$k_{\rm c}^{2}\equiv\sqrt{k^{2}_{\min}k^{2}_{\max}}=\sqrt{1-\epsilon},$ (22)

which indicates that the band is centred around $k_{\rm
c}\sim\left(1-\epsilon_{\rm c}\right)^{1/4}=1/\sqrt{3}\sim 0.57$, except for
values of $\epsilon$ extremely close to unity (see Sect. 3.3 for the physical
explanation). In that regime, this condition gives $r^{2}<0$ in Eq.18. This
ensures three and only three complex roots for the dispersion relation as
expected. Another quick way to find the critical value $\epsilon_{\rm c}=8/9$
consists of solving for $k$ as function of $s$ in Eq. 18 to get
$k(s)^{2}=-\frac{2s+8s^{2}+8s^{3}}{2s+\epsilon}$, which gives by enforcing
$r^{2}<0$ in Eq. 18, $-2s^{2}+(3\epsilon-4)s-(1-\epsilon)>0$. Positivity is
ensured for positive discriminant, i.e.

$\Delta=\epsilon(-8+9\epsilon)\geq 0\Rightarrow\epsilon\geq\frac{8}{9}.$ (23)

At the critical value $\epsilon=\epsilon_{\rm c}$, $k=1/\sqrt{3}$ and the
dispersion relation Eq. 12 factorises according to
$\left(\omega+i/3\right)^{3}=0$. Real and imaginary parts of
$\mathrm{e}^{i\omega t}$ for the different modes at
$\epsilon=0.95>\epsilon_{\rm c}$ are shown on Fig. 3. As expected, at $k=0.5$,
no mode propagates. This plot can be compared to the case $\epsilon=0.1$,
where the modes $\omega_{\pm}$ propagate for any value of $k$ (Fig. 5,
Appendix C).

Figure 3: Real (solid line) and imaginary (dotted line) parts of
$\mathrm{e}^{i\omega t}$ for the modes $\omega_{+}$, $\omega_{0}$ and
$\omega_{-}$ ($\epsilon=0.95$). The middle case $k=0.5$ corresponds to the
regime where no mode propagates. Note the different timescales used for the
plot for readability. The coefficients of the eigenvectors corresponding to
the eigenvalues are provided on top of the plots and can be found in Table.
1–2. The mode $\omega_{\rm null}$ corresponding to the null space of Eqs. 5 –
8 is provided for sanity check.

###  3.2 Eigenvectors

Fig. 4 shows the modulus (top) and the argument (bottom) of the eigenvectors
corresponding to the modes $\omega_{+}$, $\omega_{0}$, $\omega_{-}$ and
$\omega_{\rm null}$ for values of $k$ centred around the bifurcation. The
amplitudes of the modes $\omega_{\rm+}$ and $\omega_{-}$ are similar inside
the forbidden band, although they are distinct outside. Differential phases of
the mode $\omega_{0}$ do not depend on $k$, outside or inside the band.
Differential phases of the modes $\omega_{\pm}$ are also constant inside the
band. This ensures the required differential velocity from which the modes are
damped. The mode $\omega_{\rm null}$ corresponds to a steady static
perturbation on the dust density only, i.e. $\delta\rho_{\rm g}=0$, $\delta
v_{\rm g}=\delta v_{\rm d}=0$.

Figure 4: Modulus (top) and phase (bottom) of the eigenvectors at the
bifurcation for the modes $\omega_{+}$, $\omega_{0}$ and $\omega_{-}$ and
$\omega_{\rm null}$ (from left to right, $\epsilon=0.95$). Values of the
related eigen system can be found in Tables. 1 – 2.

###  3.3 Physical interpretation

In the limit $\epsilon\to 0$, the dispersion relation Eq. 12 reduces to

$\omega^{2}+i\omega-k^{2}=0.$ (24)

Eq. 24 admits the solutions $\tilde{\omega}=-i/2\pm\sqrt{k^{2}-1/4}$, showing
the existence of $k_{\max}=1/2$ below which waved are damped (in this limit,
$k_{\rm min}=0$). Physically, Eq. 24 is also the dispersion relation
associated to the reduced system

$\displaystyle\partial_{t}\delta\rho_{\rm g}$
$\displaystyle=-\rho_{\mathrm{g},0}\partial_{x}\delta v_{\rm g},$ (25)
$\displaystyle\partial_{t}\delta v_{\rm g}$
$\displaystyle=-\frac{K}{\rho_{\mathrm{g},0}}\delta v_{\rm g}-c_{\rm
s}^{2}\frac{\partial_{x}\delta\rho_{\rm g}}{\rho_{\mathrm{g},0}},$ (26)

where the gas dissipates its energy in a passive dust phase through back-
reaction. When $\left(1-\epsilon\right)$ is small but finite, large values of
$k^{2}/\omega$ can satisfy the dispersion relation, which reduces to
$i\omega^{2}-\omega k^{2}-ik^{2}(1-\epsilon)=0$. This corresponds to
propagation through the mixture and explains the finite value of $k_{\rm min}$
(e.g. Laibe & Price 2014b and App. A).

For the dustywave problem with multiple dust species, the regime identified in
this study exists when the stopping times of the different species are close
enough (we verified this fact numerically). When this condition is not
fulfilled, gas/mixture mode can propagate through one specie even if it can
not through the other one, and no forbidden band is expected. This example
provides a further situation where a physical effect observed in a mixture
with a single grain size does not occur when several sizes are considered
(e.g. Krapp et al. 2019).

##  4 conclusion

We identified and characterised the existence of a regime in the single specie
dustywave problem where waves can not propagate, neither as a gas mode or a
mixture mode. This regime develops above $\epsilon_{\rm c}=8/9$, when the
dust-to-gas ratio is sufficiently large for gas to dissipate its energy in an
independent dust phase via back-reaction. Numerical solvers can be checked in
this regime to verify their accuracy at large dust fractions.

## Acknowledgements

GL acknowledges funding from the ERC CoG project PODCAST No 864965. This
project has received funding from the European Union’s Horizon 2020 research
and innovation programme under the Marie Skłodowska-Curie grant agreement No
823823. This project was partly supported by the IDEXLyon project (contract
nANR-16-IDEX-0005) under the auspices University of Lyon. We acknowledge
financial support from the national programs (PNP, PNPS, PCMI) of CNRS/INSU,
CEA, and CNES, France. We used Mathematica (Wolfram Research, 14). We thank
the anonymous referee for a thorough and insightful report and for suggesting
the title.

## Data availability

All relevant data are given in the article.

##  Appendix A Asymptotic behaviours

When $k\gg 1$ (the physical wavelength $\lambda\equiv 2\pi c_{\rm s}t_{\rm
s}/k$ is smaller than the physical stopping length $c_{\rm s}t_{\rm s}$), Eqs.
17 – 20 provide

$\begin{cases}\omega_{+}&=\left[+\sqrt{k^{2}+3\epsilon^{2}/4-\epsilon}+\mathcal{O}\left(k^{-2}\right)\right]-i\left[\epsilon/2\hskip
13.7pt+\mathcal{O}\left(k^{-2}\right)\right],\\\
\omega_{0}&=\hphantom{\left[+\sqrt{k^{2}+3\epsilon^{2}/4-\epsilon}+\mathcal{O}\left(k^{-2}\right)\right]}-i\left[(1-\epsilon)+\mathcal{O}\left(k^{-2}\right)\right],\\\
\omega_{-}&=\left[-\sqrt{k^{2}+3\epsilon^{2}/4-\epsilon}+\mathcal{O}\left(k^{-2}\right)\right]-i\left[\epsilon/2\hskip
13.7pt+\mathcal{O}\left(k^{-2}\right)\right],\end{cases}$ (27)

or in a even more simplified form,

$\begin{cases}\omega_{+}&=\left[+k+\mathcal{O}\left(1\right)\right]-i\left[\epsilon/2\hskip
13.7pt+\mathcal{O}\left(k^{-2}\right)\right],\\\
\omega_{0}&=\hphantom{\left[+k+\mathcal{O}\left(1\right)\right]}-i\left[(1-\epsilon)+\mathcal{O}\left(k^{-2}\right)\right],\\\
\omega_{-}&=\left[-k+\mathcal{O}\left(1\right)\right]-i\left[\epsilon/2\hskip
13.7pt+\mathcal{O}\left(k^{-2}\right)\right].\end{cases}$ (28)

The evolution of the plane wave is correctly described by expanding both the
real and the imaginary parts of $\omega$ to their respective leading orders.
In this gas regime, the propagation is supported by the pressure of the gas
and damped by dust back-reaction. Since $\Re\left(w_{\pm}\right)=\pm
k+\mathcal{O}\left(1\right)$, the typical physical oscillation time is
$\sim\lambda/c_{\rm s}\ll t_{\rm s}$. After a typical time
$\sim(1-\epsilon)^{-1}t_{\rm s}$, the initial dust velocity adjusts onto the
one of the gas (mode $\omega_{0}$). Meanwhile, the gas undergoes several
oscillations that are supported by its own pressure (terms $\pm k$, modes
$\omega_{\pm}$) that are progressively damped by dust back-reaction (terms
$-i\epsilon/2$, the factor 2 accounting for dissipation by both modes). In the
limit $k\ll 1$, corresponding to $\lambda\gg c_{\rm s}t_{\rm s}$, one obtains

$\begin{cases}\omega_{+}&=+\sqrt{1-\epsilon}k-i\displaystyle\epsilon
k^{2}/2\hskip 21.1pt+\mathcal{O}\left(k^{3}\right),\\\
\omega_{0}&=\hphantom{+\sqrt{1-\epsilon}k}+i(-1+\epsilon
k^{2})+\mathcal{O}\left(k^{3}\right),\\\
\omega_{-}&=-\sqrt{1-\epsilon}k-i\displaystyle\epsilon k^{2}/2\hskip
21.1pt+\mathcal{O}\left(k^{3}\right).\end{cases}$ (29)

In this mixture regime, the propagation is supported by both gas and dust
simultaneously, and damped by an effective diffusion. The typical oscillation
time for a perturbation satisfies $\lambda/c_{\rm s}\gg t_{\rm s}$. After a
typical time $t_{\rm s}$, the gas and dust velocities have relaxed towards the
barycentric velocity of the mixture (mode $\omega_{0}$, factor $-i$). Since
the stopping time is much shorter than the oscillation time, the drag
maintains the two phases well-coupled and the system tends to oscillate at the
sound speed of the mixture $c_{\rm s}\sqrt{1-\epsilon}$, which accounts for
the inertia of the dust (Laibe & Price, 2012). Damping comes from the
$-i\displaystyle\epsilon k^{2}/2$ term, which originates from the effective
diffusion of the terminal velocity approximation (Laibe & Price, 2014b).

##  Appendix B Parameters for numerical tests

We provide parameters for numerical test – before ($k=0.1$), at ($k=0.5$) and
after ($k=2.$) the bifurcation at $\epsilon=0.95$. The eigenvalues
$\omega_{+}$, $\omega_{0}$ and $\omega_{-}$ are given in Table. 1. The
corresponding values for the eigenvectors are given in Table. 2 Evolution of
the densities and velocities are shown on Fig. 3.

$\omega$ | $k=0.1$ | $k=0.5$ | $k=2.$
---|---|---|---
$\omega_{+}$ | $0.0219514\,-0.00479354i$ | $-0.639778i$ | $1.93044\,-0.474696i$
$\omega_{0}$ | $-0.990413i$ | $-0.0665243i$ | $-0.0506079i$
$\omega_{-}$ | $-0.0219514\,-0.00479354i$ | $-0.293697i$ | $-1.93044\,-0.474696i$
Table 1: Eigenvalues for numerical tests of the dustywave regime ($\epsilon=0.95$). ${}^{\top}{X}=\left(\frac{\delta\rho_{\rm g}}{\rho_{0}},\frac{\delta\rho_{\rm d}}{\rho_{0}},\frac{\delta v_{\rm g}}{c_{\rm s}},\frac{\delta v_{\rm d}}{c_{\rm s}}\right)$ | $k=0.1$ | $k=0.5$ | $k=2.$
---|---|---|---
$X_{+}$ | $\left(\begin{array}[]{c}0.0117633\,+0.0486907i\\\ -0.20202+0.925123i\\\ 0.098324\,+0.202488i\\\ 0.223959i\end{array}\right)$ | $\left(\begin{array}[]{c}0.0388307\\\ -0.0625475\\\ -0.993722i\\\ 0.0842453i\end{array}\right)$ | $\left(\begin{array}[]{c}0.0501937\,+0.00123333i\\\ -0.00576226+0.0234334i\\\ 0.974816\,-0.214459i\\\ 0.0252485i\end{array}\right)$
$X_{0}$ | $\left(\begin{array}[]{c}0.00504115\\\ -0.00509254\\\ -0.998564i\\\ +0.0530918i\end{array}\right)$ | $\left(\begin{array}[]{c}0.0172053\\\ -0.98915\\\ -0.0457829i\\\ +0.138532i\end{array}\right)$ | $\left(\begin{array}[]{c}0.000639653\\\ -0.999645\\\ -0.000323715i\\\ +0.0266263i\end{array}\right)$
$X_{-}$ | $\left(\begin{array}[]{c}0.0117633\,-0.0486907i\\\ -0.20202-0.925123i\\\ -0.098324+0.202488i\\\ +0.223959i\end{array}\right)$ | $\left(\begin{array}[]{c}0.0790522\\\ -0.308168\\\ -0.928696i\\\ +0.190543i\end{array}\right)$ | $\left(\begin{array}[]{c}0.0501937\,-0.00123333i\\\ -0.00576226-0.0234334i\\\ -0.974816-0.214459i\\\ +0.0252485i\end{array}\right)$
$X_{\rm null}$ | $\left(\begin{array}[]{c}0\\\ i\\\ 0\\\ 0\end{array}\right)$ | $\left(\begin{array}[]{c}0\\\ i\\\ 0\\\ 0\end{array}\right)$ | $\left(\begin{array}[]{c}0\\\ i\\\ 0\\\ 0\end{array}\right)$
Table 2: Eigenvectors for numerical tests of the dustywave regime
($\epsilon=0.95$).

##  Appendix C Evolution at low $\epsilon$

Evolution of the perturbations are given for $\epsilon=0.10<\epsilon_{\rm c}$,
for a purpose of comparison with Fig. 3. The regime where no mode propagates
is not observed as expected.

Figure 5: Real (solid line) and imaginary (dotted line) parts of
$\mathrm{e}^{i\omega t}$ for the modes $\omega_{+}$, $\omega_{0}$ and
$\omega_{-}$ and $\omega_{\rm null}$ ($\epsilon=0.10$). No regime is observed
in this case and at least 2 modes propagate through the mixture.

## References

  * Ahuja (1973) Ahuja A. S., 1973, Journal of Applied Physics, 44, 4863
  * Baines et al. (1965) Baines M. J., Williams I. P., Asebiomo A. S., 1965, MNRAS, 130, 63
  * Benítez-Llambay et al. (2019) Benítez-Llambay P., Krapp L., Pessah M. E., 2019, ApJS, 241, 25
  * Booth et al. (2015) Booth R. A., Sijacki D., Clarke C. J., 2015, MNRAS, 452, 3932
  * Clair et al. (1970) Clair B. P. L., Hamielec A. E., Pruppacher H. R., 1970, Journal of the Atmospheric Sciences, 27, 308
  * Debras et al. (2020) Debras F., Mayne N., Baraffe I., Jaupart E., Mourier P., Laibe G., Goffrey T., Thuburn J., 2020, A&A, 633, A2
  * Garaud et al. (2004) Garaud P., Barrière-Fouchet L., Lin D. N. C., 2004, ApJ, 603, 292
  * Gumerov et al. (1988) Gumerov N. A., Ivandaev A. I., Nigmatulin R. I., 1988, Journal of Fluid Mechanics, 193, 53
  * Hubber et al. (2018) Hubber D. A., Rosotti G. P., Booth R. A., 2018, MNRAS, 473, 1603
  * Krapp et al. (2019) Krapp L., Benítez-Llambay P., Gressel O., Pessah M. E., 2019, ApJ, 878, L30
  * Laibe & Price (2011) Laibe G., Price D. J., 2011, MNRAS, 418, 1491
  * Laibe & Price (2012) Laibe G., Price D. J., 2012, MNRAS, 420, 2345
  * Laibe & Price (2014a) Laibe G., Price D. J., 2014a, MNRAS, 440, 2147
  * Laibe & Price (2014b) Laibe G., Price D. J., 2014b, MNRAS, 440, 2147
  * Laibe & Price (2016) Laibe G., Price D. J., 2016, DUSTYWAVE: Linear waves in gas and dust (ascl:1602.004)
  * Lebreuilly et al. (2019) Lebreuilly U., Commerçon B., Laibe G., 2019, A&A, 626, A96
  * Lorén-Aguilar & Bate (2014) Lorén-Aguilar P., Bate M. R., 2014, MNRAS, 443, 927
  * Marble (1970) Marble F. E., 1970, Annual Review of Fluid Mechanics, 2, 397
  * McKinnon et al. (2018) McKinnon R., Vogelsberger M., Torrey P., Marinacci F., Kannan R., 2018, MNRAS, 478, 2851
  * Mentiplay et al. (2020) Mentiplay D., Price D. J., Pinte C., Laibe G., 2020, MNRAS, 499, 3806
  * Moseley et al. (2019) Moseley E. R., Squire J., Hopkins P. F., 2019, MNRAS, 489, 325
  * Porth et al. (2014) Porth O., Xia C., Hendrix T., Moschou S. P., Keppens R., 2014, ApJS, 214, 4
  * Price & Laibe (2015) Price D. J., Laibe G., 2015, MNRAS, 451, 813
  * Riols & Lesur (2018) Riols A., Lesur G., 2018, A&A, 617, A117
  * Saffman (1962) Saffman P. G., 1962, Journal of Fluid Mechanics, 13, 120
  * Stoyanovskaya et al. (2018) Stoyanovskaya O. P., Glushko T. A., Snytnikov N. V., Snytnikov V. N., 2018, Astronomy and Computing, 25, 25
  * Stoyanovskaya et al. (2020) Stoyanovskaya O. P., Okladnikov F. A., Vorobyov E. I., Pavlyuchenkov Y. N., Akimkin V. V., 2020, Astronomy Reports, 64, 107
  * Wolfram Research (14) Wolfram Research I., 14, Mathematica, Version 14.0, https://www.wolfram.com/mathematica
  * Yang & Johansen (2016) Yang C.-C., Johansen A., 2016, ApJS, 224, 39