Neutrino spin oscillations in gravitational fields in higher dimensions

S.A. Alavi†,‡ ,T. Fallahi.Serish‡
‡ _Department of physics, Hakim Sabzevari University, P.O. Box 397_
_Sabzevar, Iran_
†<EMAIL_ADDRESS>alaviag.gmail.com_

Neutrino physics is one of the most active fields of research with important
implications for particle physics, cosmology and astrophysics. On the other
hand, motivated by some theories including string theory, formulation of
physical theories in more than four space-time dimensions has been the subject
of increasing attention in recent years. Interaction of neutrinos with
gravitational fields is one of the interesting phenomena which can lead to
transition between different helicity states(spin oscillations). We study
neutrino spin oscillations in Schwarzschild and RN backgrounds in higher-
dimensional gravitational fields. We calculate the transition probability as a
function of time and also study the dependence of the oscillation frequency on
the orbital radius. The results help us to better understand the behavior of
gravity and neutrinos in higher dimensions.

_Keywords_ : Neutrino spin oscillation; gravitational fields; Reissner-
Nordstrom (RN) and Schwarzschild metrics; higher dimensions.

##  1 Introduction

Neutrino physics is one of the most interesting fields of study in various
branches of physics from particle physics to cosmology and astrophysics. Since
1930, when Pauli first hypothesized the existence of a neutral particle later
called neutrino, much research was done on various aspects of this fundamental
particle, leading to the discovery of neutrinos and other developments in
subsequent years. The types of neutrinos include electron neutrinos $\nu_{e}$,
tau neutrinos $\nu_{\tau}$, muon neutrinos $\nu_{\mu}$, and their
antiparticles. Significant differences in the experimental values of neutrino
fluxes produced in reactors, accelerators and solar sources, compared to the
values predicted by Standard Models, led to several hypotheses, among which
the theory of neutrino oscillation was considered more seriously. There are
currently three possibilities for neutrino oscillations: in the first type of
oscillation, different flavors of neutrinos can be converted to each other,
which is important in the physics of solar and atmospheric neutrinos. The
second type is the oscillation or transition between neutrino helicity states
of the same flavor, called neutrino spin oscillations. The third one, which is
a combination of the previous two types, is called spin-flavor oscillations.
The interaction of neutrinos with gravitational fields is one of the
interesting phenomena that can lead to transitions between different helicity
states (spin oscillations). It was realized more than 20 years ago that
external fields drastically change the process of vacuum neutrino
oscillations. The influence of the gravitational interaction on neutrino
oscillations has been studied in many publications, see Refs. 1, 2, 3, 4. On
the other hand, according to some theories such as string theory, the
formulation of physical theories in more than four dimensions (4d) of space-
time received increasing attention in recent years. Since Einstein proposed
his general theory of relativity in 1915, a lot of research has been devoted
to unify general relativity (GR) and electromagnetism as two fundamental
interactions in nature. However, the early proposals date back to the 1920s,
through Kaluza–Klein theory to unify these interactions [5, 6] that was a
classical unified field theory built in five-dimensional space-time. Recently,
motivated by string theory as a requirement for describing a consistent theory
of quantum gravity, extra dimensions have been the subject of much attention.
Besides string theory, there are some other theories proposing the necessity
of extra dimensions.[7, 8] In this work, we study neutrino spin oscillations
in Schwarzschild and Reissner–Nordstrom (RN) backgrounds in higher dimensions.
It is worth mentioning that gravitational measurements in higher dimensions
are studied in Ref. 9.

##  2 Neutrino Spin Oscillations in Gravitational Fields

The relation between the vierbein four-velocity
$u^{a}=(u^{0},u_{1},u_{2},u_{3})$ and four-velocity of the particle in its
geodesic path, $U^{\mu}$ is [10, 11, 12, 13]

${u^{a}}=e^{a}_{\mu}U^{\mu},$ (2.1)

where $e^{a}_{\mu}$ are nonzero vierbein vectors and the components of
$U^{\mu}$in spherical coordinates are as follows:

${U^{a}}=(U^{0},U_{r},U_{\theta},U_{\varphi}),U^{0}=\frac{dt}{d\tau},U_{r}=\frac{dr}{d\tau},U_{\theta}=\frac{d\theta}{d\tau},U_{\varphi}=\frac{d\varphi}{d\tau}.$
(2.2)

In the relevant metric, $U^{\mu}$ is related to the world velocity of the
particle through $\overrightarrow{U}=\gamma\overrightarrow{v}$ where
$\gamma=\frac{dt}{d\tau}$ and $\tau$ is the proper time. To study the spin
oscillation of a particle in a gravitational field we calculate
$G_{ab}=(\overrightarrow{E},\overrightarrow{B})$ which is the analog tensor of
the electromagnetic field and is defined as

$G_{ab}=e_{a\mu;\nu}e^{\mu}_{b}U^{\nu},G_{0i}=E_{i},G_{ij}=-\varepsilon_{ijk}B_{k}.$
(2.3)

$e_{a\mu;\nu}$ are the covariant derivatives of vierbein vectors:

$e_{a\mu;\nu}=\frac{\partial
e_{a\mu}}{\partial_{x^{\nu}}}-\Gamma^{\lambda}_{\mu\nu}e_{a\lambda},$ (2.4)

where $\Gamma^{\lambda}_{\mu\nu}$ is the Christopher symbol. Geodesic equation
of a particle in a gravitational field is[14]

$\frac{d^{2}x^{\mu}}{dq^{2}}+\Gamma^{\mu}_{\sigma\upsilon}\frac{dx^{\sigma}}{dq}\frac{dx^{\upsilon}}{dq}=0,$
(2.5)

where the variable $q$ parameterizes the particle’s world line. Neutrino spin
oscillations frequency is given by the expression
$\overrightarrow{\Omega}=\frac{\overrightarrow{G}}{\gamma}$, where vector
$\overrightarrow{G}$ is defined in the following way[10, 11, 15]:

$\overrightarrow{G}=\frac{1}{2}(\overrightarrow{B}+\frac{1}{1+u^{0}}[\overrightarrow{E}\times\overrightarrow{u}]).$
(2.6)

##  3 Neutrino Spin Oscillations in Schwarzschild Metric in Higher Dimensions

Amongst the various types of black hole solutions of Einstein field equations,
a natural higher-dimensional generalization of the Schwarzschild metric, also
known as the Schwarzschild-Tangherlini metric [8], has been assumed to be
stable, like its four-dimensional counterpart. The spacetime around such a
stationary, neutral, non-rotating and spherically symmetric black hole in
$(d+1)$ dimensions is described by the most general standard form of the
spacetime element:

$ds^{2}=A^{2}(r)dt^{2}-A^{-2}(r)dr^{2}-r^{2}d\Omega^{2}_{d-1},$ (3.1)

where $d\Omega^{2}_{d-1}$ denotes the element of unit $(d-1)$ -sphere with
area $S_{d-1}$ and we have

$S_{d-1}=\frac{2\pi^{\frac{d}{2}}}{\Gamma[\frac{d}{2}]},A(r)=\sqrt{1-\frac{\mu_{0}}{r^{d-2}}}.$
(3.2)

The constant parameter $\mu_{0}$ is related to the mass of the black hole by
the following relation[16]:

$M=\frac{(d-1)S_{d-1}\mu_{0}}{16\pi G_{d+1}},$ (3.3)

where $G_{d+1}=G_{4}L^{d-3}$ is $(d+1)$-dimensional gravitational constant and
L is the size of the extra dimensions, so

$A(r)=\sqrt{1-\frac{8G_{4}ML^{d-3}\pi^{1-\frac{d}{2}}\Gamma[\frac{d}{2}]}{(d-1)r^{d-2}}}.$
(3.4)

The mass $M$ is related to the Schwarzschild radius through
$M=\frac{r_{g}}{2}$. It should be mentioned that if a gravitational radius of
a black hole is much smaller than the characteristic length of the extra
dimensions, then the black hole can be very well described by asymptotically
flat solutions[8, 17, 18, 19](see also references therein). For later
convenience, we set $G_{4}=1$ and define dimensionless variables
$x=\frac{r}{\ell_{p}},\beta=\frac{r_{g}}{\ell_{p}}$ and
$\alpha=\frac{L}{\ell_{p}}$ where $\ell_{p}$ is the Planck length.

Figure 1: (Color online) The $A(r)$ component of the Schwarzschild metric vs.
x for different dimensions $d=3-7$ and with certain values of $\beta=2.4$ and
$\alpha=1.5$.

We have plotted the $A(r)$ component of the higher dimensional Schwarzschild
metric as a function of $x$ for different spatial dimensions in Fig. 1. The
location of the event horizon is determined by the equation $A(r)=0$, so as
seen in Fig. 1, this occurs in smaller distances in higher dimensions which
asserts that for large distances gravity in four-dimensional spacetime is
stronger than higher dimensions. This fact can also be checked by noting that,
in higher dimensions, the $A(r)$ curves tend more rapidly to the $A(r)$ of
flat spacetime.
The vierbein components of the vector $u^{a}$ and $\overrightarrow{B}$ are
defined as follows [10, 11, 13]:

$u^{a}=(\gamma
A,U_{r}A^{-1},U_{\theta}r,U_{\phi}r\sin\theta),\overrightarrow{B}=(U_{\phi}\cos\theta,-U_{\phi}A\sin\theta,U_{\theta}A)$
(3.5)

To derive Eq.(3.5) we used the fact that
$\overrightarrow{U}=\gamma\overrightarrow{v}$. Let us now study a neutrino
orbiting a black hole. For simplicity, we consider only circular orbits with
constant radius $r$ $(U_{r}=\frac{dr}{d\tau}=0)$. It is important to note that
not all the orbits with arbitrary radius are stable. Because of the spherical
symmetry of the gravitational field, we constraint ourselves to the study of
the orbits lying only in the equatorial plane $(\theta=\frac{\pi}{2}$ or
equivalently $v_{\theta}=0$ so $U_{\theta}=0$ ). $\Omega_{1}$ and $\Omega_{3}$
depend on $\cos\theta$ and $v_{\theta}$ so that in this case
$\Omega_{1}=\Omega_{3}=0$ [10, 11]. It should be noted that the vierbein four
velocity now takes the subsequent form:

$u^{a}=(\gamma A,0,0,\gamma v_{\varphi}r).$ (3.6)

We also have the following forms for the gravi-electric and gravi-magnetic
fields:

$\overrightarrow{E}=(E_{0},0,0),\overrightarrow{B}=(0,-v_{\varphi}A\gamma,0),$
(3.7)

where $E_{0}=\frac{1}{A^{2}(r)}(1-\frac{rF(r)}{2A^{2}(r)})$, $v_{\varphi}$ and
$\gamma^{-1}$ are defined a long these lines[14]:

$v_{\varphi}=\frac{d\varphi}{dt}=\sqrt{\frac{F(r)}{2r}},$ (3.8)

$\gamma^{-1}=\frac{d\tau}{dt}=\sqrt{A^{2}(r)-\frac{1}{2}rF(r)},$ (3.9)

where $F(r)=\frac{dA^{2}(r)}{dr}$.

After substituting Eq. (3.4) in Eqs. (3.8) and (3.9), we have

$\gamma^{-1}=\sqrt{1-\frac{4L^{-3+d}\pi^{1-\frac{d}{2}}r^{2-d}\Gamma(\frac{d}{2})r_{g}}{-1+d}+\frac{2(2-d)L^{-3+d}\pi^{1-\frac{d}{2}}r^{2-d}\Gamma[\frac{d}{2}]r_{g}}{-1+d}},$
(3.10)
$v_{\varphi}=\sqrt{-\frac{2(2-d)L^{-3+d}\pi^{1-\frac{d}{2}}r^{-d}\Gamma[\frac{d}{2}]r_{g}}{-1+d}}.$
(3.11)

By substituting Eqs. (3.4),(3.6),(3.7),(3.10) and (3.11) in Eq. (2.6) and then
$\overrightarrow{\Omega}=\frac{\overrightarrow{G}}{\gamma}$, the only non-zero
component of the frequency of neutrino spin oscillations in higher dimensions
in Schwarzschild metric i.e .$\Omega_{2}$ is obtained:

$\Omega_{2}=(+\frac{L^{-3+d}\pi^{1-\frac{d}{2}}r^{-d}\Gamma[\frac{d}{2}]r_{g}}{-1+d}-\frac{dL^{-3+d}\pi^{1-\frac{d}{2}}r^{-d}\Gamma[\frac{d}{2}]r_{g}}{2(-1+d)}$
$-\frac{2dL^{-6+2d}\pi^{2-d}r^{2-2d}\Gamma[\frac{d}{2}]^{2}r_{g}^{2}}{(-1+d)^{2}}+\frac{d^{2}L^{-6+2d}\pi^{2-d}r^{2-2d}\Gamma[\frac{d}{2}]^{2}r_{g}^{2}}{(-1+d)^{2}})^{\frac{1}{2}}=-\frac{v_{\varphi}}{2}\gamma^{-1}$
(3.12)

We have plotted $2|\Omega_{2}|M$ in Schwarzschild metric as a function of $x$
for $d=3$ and $d>3$ in Figs. 2 and 3.

Figure 2: Neutrino spin oscillations frequency for Schwarzschild metric vs.
$x$ for $d=3$ with $\beta=2.4$ and $\alpha=1.5$.

It is seen from Fig. 2 that the neutrino spin oscillations frequency tends to
zero at long distances.

Figure 3: (color online) Neutrino spin oscillations frequencies for
Schwarzschild metric vs. $x$ for various dimensions with $\beta=2.4$ and
$\alpha=1.5$.

To draw Fig. 3, we have set $c=G=\hbar=1$. It is observed that the frequencies
of oscillations decrease and tend to zero when $x\rightarrow\infty$. It is
also seen that the maximum of frequencies (peaks) for different dimensions
have different values. By increasing the dimensions of space, the peaks occur
in smaller $x$. Furthermore, it is noted that for higher dimensions by
increasing $x$, neutrino spin oscillations frequency decreases more rapidly
than (3+1)-dimensional spacetime, which shows that for Large distances
(compared to the size of extra dimensions) gravity becomes weaker in higher
dimensions.

Figure 4: (color online) Neutrino transition Probability for Schwarzschild
metric.

Figure 4 shows the neutrino spin transition probability,
$P(t)=\sin^{2}(\Omega_{2}t)$. The period time of oscillations increases as the
dimension of space increases, which indicates that for short periods of time
the contribution of three dimensional space to the probability of neutrino
spin oscillation is more than higher dimensions.

##  4 Neutrino Spin oscillations in Reissner-Nordstrom (RN) Background in
Higher Dimensions

RN metric describes the gravitational field of a charged non-rotating black
hole, where its temporal component is defined as

$A(r)=\sqrt{1-\frac{r_{g}}{r}+\frac{Q^{2}}{r^{2}}}$ (4.1)

where $Q$ is the charge of the black hole. The temporal component in higher
dimensions takes the following form [16, 18]:

$A(r)=\sqrt{1-\frac{4r_{g}L^{d-3}\pi^{1-\frac{d}{2}}\Gamma[\frac{d}{2}]}{(d-1)r^{d-2}}+\frac{Q^{2}L^{2d-6}}{\frac{(d-1)(d-2)}{2}r^{2d-4}}}$
(4.2)

By inserting Eq. (4.2) in Eqs. (3.8) and (3.9) and doing necessary
calculations, the following expressions are obtained for the neutrino angular
velocity and $\gamma^{-1}$:

$v_{\varphi}=\frac{d\phi}{dt}=(-\frac{2L^{-6+2d}Q^{2}r^{2-2d}}{-1+d}-\frac{2(2-d)L^{-3+d}\pi^{1-\frac{d}{2}}r^{-d}r_{g}\Gamma[\frac{d}{2}]}{-1+d})^{\frac{1}{2}}.$
(4.3)

$\gamma^{-1}=(1+\frac{2L^{-6+2d}Q^{2}r^{4-2d}}{(-2+d)(-1+d)}-\frac{4L^{-3+d}\pi^{1-\frac{d}{2}}r^{2-d}r_{g}\Gamma[\frac{d}{2}]}{-1+d}$
$+\frac{2L^{-6+2d}Q^{2}r^{4-2d}}{-1+d}+\frac{2(2-d)L^{-3+d}\pi^{1-\frac{d}{2}}r^{2-d}r_{g}\Gamma[\frac{d}{2}]}{-1+d})^{\frac{1}{2}}.$
(4.4)

The only nonzero component of the frequency in higher dimensions is obtained
by using Eqs. (4.2),(4.3),(4.4),(3.6)and(3.7), and inserting in relations
(2.6) and $\overrightarrow{\Omega}=\frac{\overrightarrow{G}}{\gamma}$:

$\Omega_{2}=(+\frac{L^{-12+4d}Q^{4}r^{6-4d}}{(-1+d)^{2}}+\frac{L^{-12+4d}Q^{4}r^{6-4d}}{(-1+d)^{2}(-2+d)}+\frac{L^{-6+2d}Q^{2}r^{2-2d}}{2(-1+d)}$
$+\frac{2L^{-9+3d}\pi^{1-\frac{d}{2}}Q^{2}r^{4-3d}r_{g}\Gamma[\frac{d}{2}]}{(-1+d)^{2}}+\frac{2L^{-9+3d}\pi^{1-\frac{d}{2}}Q^{2}r^{4-3d}r_{g}\Gamma[\frac{d}{2}]}{(-1+d)^{2}(-2+d)}$
$-\frac{2dL^{-9+3d}\pi^{1-\frac{d}{2}}Q^{2}r^{4-3d}r_{g}\Gamma[\frac{d}{2}]}{(-1+d)^{2}}-\frac{dL^{-9+3d}\pi^{1-\frac{d}{2}}Q^{2}r^{4-3d}r_{g}\Gamma[\frac{d}{2}]}{(-1+d)^{2}(-2+d)}$
$+\frac{L^{-3+d}\pi^{1-\frac{d}{2}}r^{-d}r_{g}\Gamma[\frac{d}{2}]}{-1+d}-\frac{dL^{-3+d}\pi^{1-\frac{d}{2}}r^{-d}r_{g}\Gamma[\frac{d}{2}]}{2(-1+d)}$
$-\frac{2dL^{-6+2d}\pi^{2-d}r^{2-2d}r_{g}^{2}\Gamma[\frac{d}{2}]^{2}}{(-1+d)^{2}}$
$+\frac{d^{2}L^{-6+2d}\pi^{2-d}r^{2-2d}r_{g}^{2}\Gamma[\frac{d}{2}]^{2}}{(-1+d)^{2}})^{\frac{1}{2}}=-\gamma^{-1}\frac{v_{\varphi}}{2}.$
(4.5)

To proceed further, we define a new variable $\eta=\frac{Q}{\ell_{p}}$. In
Fig. 5, we have plotted $2|\Omega_{2}|M$ in three dimensional RN metric as a
function of $x$ for $\beta=2.4,\eta=0.6$ and $\alpha=1.5$.

Figure 5: Neutrino spin oscillations frequency for RN metric vs. $x$ in three
dimension.

$2|\Omega_{2}|M$ in three as well as higher-dimensional RN metric as a
function of $x$ are presented in Fig. 6.

Figure 6: Neutrino spin oscillations frequency vs. $x$ for RN metric in higher
dimensions.

There are some points that can be inferred from Fig. 6: ( i ) There is a
regular peak for the value of frequency of oscillation in each spacetime
dimension. ( ii ) The maximum values (peaks) reduce by decreasing the
spacetime dimensions and the peaks of the frequencies take place in smaller x
as the number of dimensions increase. ( iii ) As we get away from the horizon,
the values of frequencies go faster to zero by increasing spacetime
dimensions, which means that, in long distances, the dominant contribution
belongs to three dimension, but is short distances (small $x$), the
contributions of higher dimensions become stronger.

Figure 7: Neutrino spin transition Probability for RN metric in Higher
dimensions.  Figure 8: Neutrino spin transition Probability for
$Q=0,0.25M,0.5M$ in three dimension.

In Fig. 7, the neutrino spin transition probability $P(t)$ for RN metric is
plotted in three together with higher dimensions. It is seen that the behavior
is the same as Schwarzschild metric i.e., by increasing the dimension of
space, the time period of oscillations increases. Fig. 8 shows the effects of
charge of the black hole on neutrino spin transition probability $P(t)$ in
three dimensions. It is observed that by increasing the black hole charge the
time period of oscillations decreases.

Figure 9: The effects of the black hole charge on neutrino spin oscillation
frequency for RN metric in different dimensions with $\eta=0.7$.

Figure 9 displays the effects of charge of the black hole on $\Omega_{2}$ in
different dimensions. the following points can be deduced from the figure: (i)
By increasing the dimension of space, the effects of the black hole charge
decrease and diminish to zero. (ii) In each dimension by increasing the charge
of the black hole, the peak occurs in smaller distances and finally the $Q=0$
curves are in agreement with the results obtained for the Schwarzschild
metric.

##  5 Phenomenological Application

Let us present a phenomenological application of our results. We consider a
simple neutrino bipolar system which helps us to understand many of the
qualitative features of collective neutrino oscillations in supernovae. The
system composed of a homogeneous and isotropic gas that initially consists of
mono-energetic $\nu_{e}$ and $\overline{\nu_{e}}$ and is described by the
flavor pendulum [20]. We introduce $\epsilon$ as the fractional excess of
neutrinos over antineutrinos, $\epsilon=\frac{n_{\nu}}{n_{\overline{\nu}}}$.
Now let us consider such systems in the presence of a neutral (Schwarzschild)
black hole with the same initial value of $\epsilon$ in different dimensions:

$\epsilon^{d=3}_{t=0}=\epsilon^{d=4}_{t=0}=.....\epsilon^{d=7}_{t=0}$ (5.1)

It can be checked from Fig. 4 and 7 that the time period of the neutrino spin
oscillations probability increases with increasing dimensions:

$T_{d=3}<T_{d=4}<T_{d=5}<T_{d=6}<T_{d=7},$ (5.2)

so at a later time $t>0$, we have

$\epsilon^{d=3}_{t}<\epsilon^{d=4}_{t}<\epsilon^{d=5}_{t}<\epsilon^{d=6}_{t}<\epsilon^{d=7}_{t}.$
(5.3)

The same result is obtained for the RN metric. This implies that the spin
oscillations frequencies of the flavor pendulum as a function of the neutrino
number density for both metrics decrease in higher dimensions.

##  6 Discussion

In this section, we discuss some more details and clarify some points. We
start with a discussion about the values of neutrino energy in the problem of
spin oscillation around the black holes. There is a wide range of neutrino
energies from one-millionth of an electronvolt to $10^{18}$ electronvolt. This
implies that there is a great deal of neutrinos to explore and getting
interesting information about the processes that formed those neutrinos. Some
of the lowest-energy neutrinos come from the Big Bang (cosmic neutrino
background or relic neutrinos), while the most energetic ones reach us from
extragalactic sources. It is essential to recognize that the cosmic neutrino
background (CNB) is the main known source of non-relativistic neutrinos in the
universe. There have been much interest and efforts to detect these non-
relativistic neutrinos here on earth, see Ref. 21. So in principle, there is a
broad possible range of energy for neutrinos but, it is worth mentioning that,
usually the spin-flip induced by only gravitation is strongly suppressed in
the relativistic regimes and in fact, a magnetic field is needed.[22] A
comment about the phenomenological application is also in order. We have
interpreted the neutrino spin precession as neutrino-antineutrino
oscillations, which is true for Majorana neutrinos. In other words, for
Majorana neutrinos, the neutrino spin flip is equivalent to neutrino-
antineutrino oscillations. For the Dirac case, only the left-handed neutrinos
are available and the right-handed neutrinos’ population is not included in
the CNB or other neutrino energy ranges. Therefore, for the Majorana case, the
neutrino flux contains both left and right-helical neutrinos, and the capture
rate is increased. For Dirac neutrinos, after spin transition from left-handed
to right-handed neutrinos, the latter becomes sterile, and hence does not
interact. Consequently, the related flux should be modified by dimensional
effects, along the lines discussed by Dvornikov.[23, 24]

##  7 Conclusion

In this paper, we have investigated neutrino spin oscillations in
gravitational fields in higher dimensions. We have analyzed the dependence of
the oscillation frequency on the orbital radius. For both Schwarzschild and RN
metrics, the maximum of the frequency increases by increasing the dimensions
of space and the peaks occur in smaller distances compared to the three-
dimensional case. It is also observed that for long distances the main
contribution to the neutrino spin transition probability belongs to three
dimensions. In addition, in both metrics, the period time of oscillation
increases in extra dimensions. It is also shown that the effects of the black
hole charge on the frequency of spin oscillations decrease and tend to zero in
higher dimensions. We have also briefly studied the effects of higher
dimensions on a bipolar neutrino system.

## References

  * [1] D. Piriz, M. Roy and J. Wudka,  Phys. Rev. D, 54, 1587 (1996).
  * [2] C. Y. Cardall and G. M. Fuller,  Phys. Rev. D, 55, 7960 (1997).
  * [3] H. J. M. Cuesta and G. Lambiase,  , Astrophys. J. 689, 371 (2008).
  * [4] G. Lambiase, G. Papini, R. Punzi and G. Scarpetta,  Phys. Rev. D, 71, 073011 (2005).
  * [5] T. Kaluza,,  Sitz. Kgl. Preuß. Akad. Wiss. (Berlin), 1921, 966 (1921).
  * [6] O. Klein,  Z. Phys. 37, 895 (1926).
  * [7] M. Quirós,  arXiv:hep-ph/0606153.
  * [8] F. R. Tangherlini,  Il Nuovo Cimento, 27, 636 (1963).
  * [9] D. M. Yekta, S. A. Alavi and M. Karimabadi,  Galaxies, 9, 4 (2021).
  * [10] M. Dvornikov,  Int. J. Mod. Phys. D, 15, 1017 (2006).
  * [11] S. A. Alavi and S. F. Hosseini,  Gravit. Cosmol. 19, 129 (2013).
  * [12] M. Dvornikov,  Phys. Rev. D99, 116021 (2019).
  * [13] S. A. Alavi and S. Nodeh,  Phys. Scripta 90 , 035301 (2015).
  * [14] S. Weinberg,  Gravitation and Cosmology: Principles and Applications of the General Theory of Relativity, (John Wiley and Sons, 1972).
  * [15] G. Lambiase and L. Mastrototaro,  Phys. Rev. D 104, 024021 (2021).
  * [16] R. C. Myers,  Phys. Rev. D 60, 046002 (1999).
  * [17] K. Nouicer,  Class. Quantum Grav. 24, 5917(2007).
  * [18] R. C. Myers and M. J. Perry,  Ann. Phys. 172, 304(1986).
  * [19] R. Emparan and H. S. Reall,  Living Rev. Relativ. 11, 1(2008).
  * [20] H. Duan, G. M. Fuller, J. Carlson and Y. Z. Qian,  Phys. Rev. D, 75, 125005 (2007).
  * [21] A. J. Long, C. Lunardini and E. Sabancilar,  J. Cosmol. Astropart. Phys, 2014, 038(2014).
  * [22] M. Dvornikov,  J. Cosmol. Astropart. Phys.2021, 005 (2021).
  * [23] M. Dvornikov,  , Eur. Phys. J. C 80, 1 (2020).
  * [24] M. Dvornikov,  J. Cosmol. Astropart. Phys. 2013, 015 (2013).


generalized eigenfunctions for Maxwell’s equations using a basis of vector
spherical harmonics.

Acknowledgements: The authors wish to thank Alexander Strohmaier for useful
discussions during the course of this paper. Part of this research was carried
while the author A. W. was at the Erwin Schrödinger Institute Thematic Program
”Spectral Theory and Mathematical Relativity.” A.W. wishes to thank the Erwin
Schrödinger Institute for its hospitality (and coffee!) during this time.

## References

  * [1] H. Bahouri, J.-Y. Chemin and R. Danchin.  Fourier Analysis and Nonlinear Partial Differential Equations.  Grundlehren der mathematischen Wissenschaften, 2011.
  * [2] C. Bardos, G. Lebeau, and J. Rauch.  Sharp sufficient conditions for the observation, control and stabilization of waves from the boundary.  SIAM Journal of Control and Optimization, (305):1024– 1065, 1992.
  * [3] N. V. Borisov, W. Mul̈ler, and R. Schrader.  Relative index theorems and supersymmetric scattering theory.  Communications in mathematical physics 114(3):475–513, 1988.
  * [4] N. Burq and R. Schippa  Strichartz estimates for Maxwell equations on domains with perfectly conducting boundary conditions  preprint arXiv:2304.13368, 2023
  * [5] N. Burq, P. Gerard, and N. Tzvetkov.  Strichartz inequalities and the nonlinear Schrödinger equation on compact manifolds.  Amer. J. Math., 126(3):569–605, 2004.
  * [6] P. E. Conner.  The Green’s and Neumann’s problems for differential forms on Riemannian manifolds.  Proceedings of the National Academy of Sciences, 40(12):1151–1155, 1954.
  * [7] T. Duyckaerts and D. Lafontaine  Scattering for critical radial Neumann waves outside a ball  Revista Matemática Iberoamericana, 2021.
  * [8] G. F. D. Duff and D. C. Spencer  Harmonic Tensors on Riemannian Manifolds with Boundary.  Annals of Mathematics, 56(1):128–156, 1952.
  * [9] S. Dyatlov and M. Zworski  Mathematical Theory of Scattering Resonances.  Graduate Studies in Mathematics, 2019.
  * [10] A. Ericsson, D. Sjöberg, C. Larsson, T. Martin  Scattering from a multilayered sphere- Applications to electromagnetic absorbers on double curved surfaces  Technical Report, Lund University, Sweden. Vol 7249, 2017.
  * [11] M. P. Gaffney.  The heat equation method of Milgram and Rosenbloom for open Riemannian mani- folds.  Annals of Mathematics 458–466, 1954.
  * [12] J. Galkowski, P. Marchang, J. Wang, and M. Zworski  The scattering phase seen at last  https://arxiv.org/abs/2210.09908
  * [13] C. Guillarmou and D. Sher.  Low energy resolvent for the Hodge Laplacian: Applications to Riesz transform, Sobolev estimates and analytic torsion.  International Mathematics Research Notices, 1–75, 2015.
  * [14] C. Guillarmou, A. Hassell  The resolvent at low energy and Riesz transform for Schrodinger operators on asymptotically conic manifolds, Part I.  Math Annalen. 341(4), 859-896. 2018.
  * [15] C. Guillarmou, A. Hassell,  The resolvent at low energy and Riesz transform for Schrodinger operators on asymptotically conic manifolds, Part II  Ann. Inst. Fourier 59(2), 1553-1610, 2009
  * [16] C. Guillarmou, A. Hassell,  Uniform Sobolev estimates for non-trapping metrics,  Journal of Inst. Math. Jussieu, 13(3), 599-632, 2014.
  * [17] A. Kirsch and F. Hettlich.  The Mathematical theory of time-harmonic Maxwell’s equations  Springer International. Switzerland 2015.
  * [18] T. Harge and G. Lebeau.  Diffraction par un convexe.  Invent. Math., 118:161–196, 1994.
  * [19] P. Hintz,  A Sharp Version of Price’s Law for Wave Decay on Asymptotically Flat Spacetimes,  Commun. Math. Phys., 389, 491-542, 2022.
  * [20] L. Hörmander, _The Analysis of Linear Partial Differential Operators I: Distribution Theory and Fourier Analysis_ Springer-Verlag. Copyright 1963. Reprint of the 1990 edition, 2003.
  * [21] O. Ivanovici.  On the Schrödinger equation outside strictly convex obstacles.  Anal. PDE, 3(3):261–293, 2010.
  * [22] O.Ivanovici and G.Lebeau  Dispersion estimates for the wave and the Schrödinger equations outside a ball and counterexamples  preprint arXiv:2012.08366, 2021.
  * [23] O.Ivanovici and G.Lebeau  Dispersion for the wave and the Schrödinger equations outside strictly convex obstacles and counter examples  C.R. Acad. Sci. Paris, Ser. I. 355(7) 774–779, 2017.
  * [24] P. Lax and R. Phillips  Scattering theory  Academic Press. New York and London, 1967
  * [25] R. B. Melrose.  Singularities and energy decay in acoustical scattering.  _Duke Mathematical Journal_ , 46(1):43–59, 1979.
  * [26] R. B. Melrose  Geometric Scattering Theory. Cambridge, UK: Cambridge University Press, 1995.
  * [27] R. B. Melrose and M. E. Taylor.  Near peak scattering and the corrected Kirchhoff approx- imation for a convex obstacle.  _Adv. in Math._ , 55(3):242–315, 1985.
  * [28] R. B. Melrose and M. E. Taylor.  The radiation pattern of a diffracted wave near the shadow boundary.  _Comm. Partial Differential Equations_ , 11(6):599–672, 1986.
  * [29] R. B. Melrose and M.E. Taylor.  Boundary problems for the wave equations with grazing and gliding rays, 1987.
  * [30] G. Mie  Beiträge zur Optik trüber Medien, speziell kolloidaler Metallösungen  Annalen der Physik. 330 (3): 377–445, 1908
  * [31] K. Morgan and J. Wunsch,  Generalized Price’s law on fractional-order asymptotically flat stationary spacetimes.  Mathematical Research letters.
  * [32] W. Mueller and A. Strohmaier,  Scattering at low energies on manifolds with cylindrical ends and stable systoles.  Geom. Funct. Anal. 20(3):741–778. 2010
  * [33] F. W. J. Olver,  Asymptotics and special functions  Academic Press, Inc. San Diego, 1974.
  * [34] F. W. J. Olver, D. W. Lozier, R. F. Boisvert, and C. W. Clark.  Nist digital library of mathematical functions.  Online companion to: http://dlmf. nist. gov, 2010.
  * [35] V. Petkov and M. Zworski.  Semi-classical estimates on the scattering determinant.  In Annales Henri Poincare 2, pages 675–711. Springer, 2001.
  * [36] L. Parnovski  Scattering matrix for manifolds with conical ends.  J. London Math. Soc. 61(2):555–567. 2000.
  * [37] J. Shatah and M. Struwe  Geometric wave equations.  Courant Lecture Notes in Mathematics, (1998).
  * [38] C. D. Sogge and S. Zelditch  Riemannian manifolds with maximal eigenfunction growth.  Duke Math. J., 114(3):387–437, 2002.
  * [39] H. Smith and C. Sogge.  On the critical semilinear wave equation outside convex obstacles  Journal of the American Mathematical Society 8(4) 879-916,1995
  * [40] J.A. Stratton  Electromagnetic Theory  New York: McGraw-Hill, 1941
  * [41] A. Strohmaier and A. Waters.  Geometric and obstacle scattering at low energy.  Comm. Partial. Diff. Eq., 45(11): 1451-1511, 2020.
  * [42] A. Strohmaier  The classical and quantum photon field for non-compact manifolds with boundary and in possibly inhomogeneous media.  Comm. Math. Phys., 1–21, 2021.
  * [43] A. Strohmaier and A. Waters. _The Birman-Krein formula for differential forms and electromagnetic scattering_. 17, paper no. 103166, 2022.
  * [44] M. Zworski.  High frequency scattering by a convex obstacle.  Duke Math. J., 61(2):545–634, 1990.


# Upmixing via style transfer: a variational autoencoder for disentangling
spatial images and musical content

###### Abstract

In the stereo-to-multichannel upmixing problem for music, one of the main
tasks is to set the directionality of the instrument sources in the
multichannel rendering results. In this paper, we propose a modified
variational autoencoder model that learns a latent space to describe the
spatial images in multichannel music. We seek to disentangle the spatial
images and music content, so the learned latent variables are invariant to the
music. At test time, we use the latent variables to control the panning of
sources. We propose two upmixing use cases: transferring the spatial images
from one song to another and blind panning based on the generative model. We
report objective and subjective evaluation results to empirically show that
our model captures spatial images separately from music content and achieves
transfer-based interactive panning.

Index Terms—  Stereo-to-multichannel upmixing, variational autoencoders,
panning, information disentanglement

© 2022 IEEE. Personal use of this material is permitted. Permission from IEEE
must be obtained for all other uses, in any current or future media, including
reprinting/republishing this material for advertising or promotional purposes,
creating new collective works, for resale or redistribution to servers or
lists, or reuse of any copyrighted component of this work in other works.

##  1 Introduction

While home theater and surround sound have become more common and affordable,
most music is still in stereo. Approaches to play stereo content on
multichannel speakers can be categorized into two kinds. The first uses matrix
[1] to map the left and right content channels into each of the surround
speakers. These operate in different sound modes, e.g. direct, stereo, and all
channel stereo mode, without sources decomposition. The second approach relies
on the separation of sound components from the stereo mixture, e.g., into
music content and other ambient sounds [2, 3] or into individual musical
instruments [4] that are then placed in the desired panning locations. There
are different techniques proposed to the primary-ambient separation, including
frequency-domain inter-channel coherence index computation [5, 6, 7], time or
frequency-domain adaptive-filtering [8, 9], PCA-based primary-ambient
decomposition [10, 11], and deep neural network (DNN)-based methods [12, 13,
14]. Perceived azimuth direction has also proven to be beneficial for direct
and ambient separation [15]. Jeon et al. proposed an upmixing algorithm based
on robust source separation with a post-scaling algorithm to compensate for
the interference [4]. Another DNN-based upmixing model, introduced by Park
et.al [16], bypasses any explicit decomposition and generates the center and
surround channels from stereo channels directly based on non-linear
transformations learned in the DNN.

In this paper, we develop a generative model for upmixing. The goal is to
implicitly separate the instrumental components contained in the stereo signal
and rearrange them into the five-channel setting that renders the desired
perceptual source directions. This setting of the problem leads to the fact
that there can be many legitimate upmixed versions that correspond to a stereo
signal, thus making optimization ill-defined on this one-to-many mapping
function .

To address this issue, we firstly make a basic assumption, that the spatial
images of a multi-channel music signal can be represented independently of the
musical content. This means that two different 5-channel music signals with
the same instrumentation and panning will map to an identical spatial
representation even though the content—melody, harmony, etc—differs. Then, we
hope to extract the source-specific spatial representation from the
multichannel audio and utilize the learned representation to guide the test-
time upmixing.To this end, we propose a modified variational autoencoder (VAE)
model [17] and train it in a supervised manner. Similar to the ordinary VAE
training, we use a 5-channel signal as both input and target. However, we also
provide a stereo version of the 5-channel signal as input to the decoder and
train the VAE bottleneck layer to exclusively capture the spatial images of
the sources rather than an entangled representation of music components and
spatial information.

We empirically show that the learned latent representation reflects the
spatial images of the multichannel input, i.e., it is correlated to the
panning strategy used to render the 5-channel input, and is invariant to the
music content. Hence, the decoder of the trained VAE can generate the upmixed
version of the input stereo with the guidance provided in the form of the
latent variables. Moreover, we also show that panning information represented
in the latent space can be transferred in-between different songs, enabling
spatial style transfer.

##  2 Upmixing: a probabilistic formulation

We firstly describe a general probabilistic upmixing model where spatial
information and the stereo music signal are still entangled. Then, we showcase
how the assumption that spatial information and stereo audio are independent
alters the proposed formulation. We further introduce an alternative
deterministic downmixing process, to achieve the desired disentanglement.

We define the upmixing process probabilistically, as the likelihood of
observing the 5-channel sample $\bar{\bm{X}}^{(5)}$ given the spatial
information $\bar{\bm{h}}$ and stereo audio $\bar{\bm{X}}^{(2)}$:
$p_{\theta}(\bm{X}^{(5)}|\bm{X}^{(2)},\bm{h})$, where $\theta$ stands for the
model parameters and the bar notation is the sample of a random variable,
e.g., $\bar{\bm{X}}^{(5)}\sim\bm{X}^{(5)}$. We use superscript to denote the
number of channels. However, we cannot compute the likelihood due to the
intractable marginal distribution $p(\bm{X}^{(5)})$ and posterior distribution
$p(\bm{X}^{(2)},\bm{h}|\bm{X}^{(5)})$. VAE employs variational inference to
resolve this issue by proposing an encoder function
$q_{\phi}(\bm{X}^{(2)},\bm{h}|\bm{X}^{(5)})$ that best approximates the true
posterior. Coupled with concept of evidence lower bound (ELBO), the
optimization is eventually defined as follows:

$\begin{split}\mathcal{L}(\theta,\phi;\bm{X}^{(5)})=&\mathbb{E}_{q_{\phi}(\bm{X}^{(2)},\bm{h}|\bm{X}^{(5)})}[\log
p_{\theta}(\bm{X}^{(5)}|\bm{X}^{(2)},\bm{h})]\\\
-&D_{KL}\big{(}q_{\phi}(\bm{X}^{(2)},\bm{h}|\bm{X}^{(5)})||p(\bm{X}^{(2)},\bm{h})\big{)},\end{split}$
(1)

where $p(\bm{X}^{(2)},\bm{h})$ denotes the prior distribution of the latent
variables, thus forming a regularization term. Meanwhile, the main
reconstruction objective is to maximize the expectation of the log-likelihood
w.r.t. the approximated posterior distribution.

The desired disentanglement of stereo music content and spatial information
assumes that $\bm{X}^{(2)}$ and $\bm{h}$ are independent from each other.
Therefore, the joint posterior distribution can be factorized, i.e.,
$p(\bm{X}^{(2)},\bm{h}|\bm{X}^{(5)})=p(\bm{X}^{(2)}|\bm{X}^{(5)})p(\bm{h}|\bm{X}^{(5)})$.
To further constrain that $\bm{X}^{(2)}$ only reflects stereo music content
while $\bm{h}$ represents 5-channel spatial information, we employ a hard
regularization trick by replacing the random variable $\bm{X}^{(2)}$ with a
simple downmixed version $\bar{\bm{X}}^{(2)}$. Then, the probabilistic model
$p(\bm{X}^{(2)}|\bm{X}^{(5)})$ is simplified into a deterministic downmixing
process, $\bar{\bm{X}}^{(2)}=\mathcal{G}(\bar{\bm{X}}^{(5)})$, where
$\mathcal{G}$ denotes the pre-defined downmix function. We treat
$p(\bm{X}^{(2)}|\bm{X}^{(5)})$ as a constant and rewrite eq. (1) as follows:

$\begin{split}\mathcal{L}(\theta,\phi;\bm{X}^{(5)},\bm{X}^{(2)})\propto&\mathbb{E}_{q_{\phi}(\bm{h}|\bm{X}^{(5)})}[\log
p_{\theta}(\bm{X}^{(5)}|\bar{\bm{X}}^{(2)},\bm{h})]\\\
-&D_{KL}\big{(}q_{\phi}(\bm{h}|\bm{X}^{(5)})||p(\bm{h})\big{)}.\end{split}$
(2)

The loss function represents our proposed disentanglement method. As a
modified VAE, the spatial information is learned via $\bm{h}$ while the
decoder takes in extra stereo input to generate the multi-channel output
$\bm{X}^{(5)}$. The second term regularizes each latent variable by a standard
normal distribution, i.e., $p(\bm{h})=\mathcal{N}(\bm{0},\bm{I})$, where
$\bm{I}$ is a $\\!J$ dimensional identical matrix.

Fig. 1: Model architecture

##  3 Model Description

###  3.1 DNN architecture

Fig. 1 summarizes the overall model architecture. We employ densely connected
convolutional blocks introduced in DenseNet [18] as our building blocks in
both the encoder and the decoder. In both modules, we alternate 5 dense blocks
with 4 transition layers. A dense block contains 5 convolutional layers, each
of which produces $C=20$ output channels. We set the stride to be 2 for the
transition layers in the encoder to reduce the data dimension. During
training, the encoder takes as input a 5-channel audio in the form of stacked
magnitude spectrograms, $\bar{\bm{X}}_{r}^{(5)}\in\mathbb{R}_{+}^{5\times
F\times T}$, where $F$ and $T$ are the frequency subbands and time-domain
frames. Here, we introduce a subscript $r$ to denote the particular panning
strategy for creating a 5-channel input signal. The encoder outputs two
$J$-dimensional vectors $\bm{\mu}_{r},\bm{\sigma}_{r}\in\mathbb{R}^{J}$, which
are the parameters of the $J$-dimensional multivariate latent normal
distribution. Latent variables are then sampled using the reparameterization
technique [17], i.e.,
$\bar{\bm{h}}_{r}=\bm{\mu}_{r}+\bm{\epsilon}\odot\bm{\sigma}_{r}$ where
$\bm{\epsilon}\sim\mathcal{N}(\bm{0},\bm{I})$.

The stacked stereo spectrogram, $\bar{\bm{X}}^{(2)}\in\mathbb{R}^{2\times
F\times T}$, is a unique feature of the proposed model. It is input to the
decoder to provide the music signal for upmixing. During training, we first
prepare a 5-channel signal $\bar{\bm{X}}_{r}^{(5)}$ under a panning condition
$r$, from which the encoder learns the latent vector ${\bm{h}}_{r}$.
Meanwhile, the same music stems are rearranged into another 5-channel version
$\bar{\bm{X}}_{a}^{(5)}$ using a different panning condition $a$, whose stereo
downmix $\bar{\bm{X}}_{a}^{(2)}$ is fed to the decoder. We rearrange the
source locations in $\bar{\bm{X}}_{a}^{(2)}$ rather than reusing
$\bar{\bm{X}}_{r}^{(2)}$, to ensure that the left- and right panning in the
stereo is unrelated to the left- and right panning in the 5-channel signal.

We repeat $\bar{{\bm{h}}}_{r}$ and stack to $\bar{\bm{X}}_{a}^{(2)}$ at the
concatenation step to condition the spatial information on the stereo input.
The transition layer that receives the stereo $\bar{\bm{X}}_{a}^{(2)}$ does
not change the tensor size, which is $2\\!\times\\!F\\!\times\\!T$. After
repeating the values of $J$-dimensional embedding vector $\bar{{\bm{h}}}_{r}$
at all $F\\!\times\\!T$ bins, the input to the decoder is of size
$(2\\!+\\!J)\\!\times\\!F\\!\times\\!T$.

###  3.2 The Test-Time Generation Processes

For test-time generation, the decoder operates as the upmix generator, which
takes stereo signal $\bm{X}_{a}^{(2)}$ and the spatial information
${\bm{h}}_{r}$ as a seed. Our model can accommodate two use-case scenarios:

  * •

Style transfer-based upmixing: We posit that panning style transfer from one
music piece to another is possible. Here we define the panning style as the
set of apparent source directions for each instrument. In this case, the user
provides two input signals: $\bar{\bm{Y}}_{s}^{(5)}$ as the source 5-channel
audio with the panning style $s$ and $\bar{\bm{Z}}^{(2)}$ as the stereo music
to upmix. Note that during the test time we do not specify the panning method
of $\bar{\bm{Z}}^{(2)}$ as it doesn’t influence the process. The actual music
in $\bar{\bm{Y}}_{s}^{(5)}$ and $\bar{\bm{Z}}^{(2)}$ differs and the aim of
the generator is to create a 5-channel audio $\hat{\bm{Z}}_{s}^{(5)}$ which is
upmixed from $\bar{\bm{Z}}^{(2)}$ and has the same panning method as
$\bar{\bm{Y}}_{s}^{(5)}$. The proposed style transfer is conducted by sharing
the latent variables $\bm{h}_{s}$ in style extraction (i.e. encoding) and
upmixing (i.e. decoding) processes. The encoder takes in
$\bar{\bm{Y}}_{s}^{(5)}$ and learns $\bm{h}_{s}$ which encodes the source’s
spatial panning. The generator takes in stereo seed signal
$\bar{\bm{Z}}^{(2)}$ and $\bar{\bm{h}}_{s}$, and synthesizes the multi-channel
version of the target $\hat{\bm{Z}}_{s}^{(5)}$. This indirect way enables a
user interface where the user handles spatial control by providing an example
song.

  * •

Blind upmixing: The blind upmixing method literally “generates” spatial
panning information from the learned latent variables $\bm{h}$ and applies it
to the seed stereo signal $\bar{\bm{Z}}^{(2)}$. In our VAE, the generation
process begins with a random sample $\bar{\bm{h}}$ only, hence the spatial
image of the generated 5-channel audio sample $\hat{\bm{Z}}^{(5)}$ is out of
the user’s control. In the blind upmixing scenario, the sampling process is
completely random without involving any intuition about the latent space.

###  3.3 The Baseline

Although blind upmixing has been studied extensively in the literature, our
approach to upmixing via style transfer is novel. Hence, we are unable to use
any existing model as a baseline. Instead, we build a simple baseline model by
spreading each channel of $\bar{\bm{Z}}^{(2)}$ to the front and rear channels
of the same side in the 5-channel output. This is the same as the “all channel
stereo” upmixing mode commonly used by home receivers. It is a straightforward
way to perform upmixing and can perfectly preserve the sound quality of the
stereo input because the process does not generate any artifacts. By comparing
with the baseline model, we mainly seek to prove the functionality of our
model on style-transfer upmixing.

##  4 Experimental Setup

###  4.1 Vector-Based Amplitude Panning for Data Preparation

To train and evaluate the model, it is necessary for us to know the panning
method that is used to create the 5-channel audio, so as to examine whether
the information captured in the latent space is correlated to the correct
spatial images that we want to learn. Also, when having stem sources, we can
expose the model to various artificial upmixing configurations where source
locations change freely. However, such ground-truth spatial maps are not
readily available for most 5-channel music signals we have access to.

To that end, we choose to build our own 5-channel dataset from individual
musical instruments via the vector-based amplitude panning (VBAP) method [19].
VBAP provides an efficient equation for virtual sound source positioning,
enabling control of an unlimited number of loudspeakers in an arbitrary two-
or three-dimensional placement around the listener. In this paper, we employ a
two-dimensional rendering space with 5 speakers, following the ITU’s 5.1
channel standard [20] without the subwoofer, which is not considered in our
upmix algorithm. For each instrument, we specify a virtual source direction.
Then, we pan each instrument independently using the two adjacent speakers to
the desired coming direction, based on the vector base formulation [19].

###  4.2 Datasets

We build a synthesized 5-channel dataset using the MUSDB18 dataset [21].
MUSDB18 provides pop songs in stereo format and four separated instrumental
stems: vocals, drums, bass, and other. We split data set into training,
validation, and testing subsets, which amount to about 5, 0.5, and 1.5 hours,
respectively. The 5-channel versions are created using the VBAP method
described in Sec. 4.1 by controlling the stem tracks. Stereo input signals
$\bar{\bm{X}}^{(2)}$ are rendered by applying passive downmixing on the
synthesized 5-channel, i.e.,
$\bar{\bm{X}}^{(2)}=\mathcal{G}(\bar{\bm{X}}^{(5)})$, where
$\bar{\bm{X}}^{(2)}_{[1,:,:]}=\bar{\bm{X}}^{(5)}_{[1,:,:]}+\bar{\bm{X}}^{(5)}_{[2,:,:]}+\bar{\bm{X}}^{(5)}_{[3,:,:]}/2$
and
$\bar{\bm{X}}^{(2)}_{[2,:,:]}=\bar{\bm{X}}^{(5)}_{[4,:,:]}+\bar{\bm{X}}^{(5)}_{[5,:,:]}+\bar{\bm{X}}^{(5)}_{[3,:,:]}/2$.
The order of the five channels in $\bar{\bm{X}}^{(5)}$ are : front left, real
left, center, front right, and real right.

We mainly run our experiments on the synthesized MUSDB18 dataset. However,
when validating the blind upmixing model, we employ another internal real-
world dataset of professionally mixed surround music with approximately 17
hours for training, 2.5 hours for validation, and 4 hours for testing, in
order to test our model’s generalization ability on real music data.

###  4.3 Training Setup

The directions of sources are randomly sampled from the entire $360^{\circ}$
circle. The models are trained on 2.2 second-long segments, which have at
least one instrument not silent. Each segment is processed by short-time
Fourier transform (STFT) on Hann-windowed frames of 1024 samples with a 75%
overlap, resulting in a spectrogram of size $513\times 384$. We apply the
phase of the stereo input to recover 5-channel magnitude spectrograms. The
left channel’s phase in stereo is used for both front and rear left channels
in the 5-channel audio, and similarly for the right channels; the mean of the
left and right channels is used for the center channel recovery. We use the
Adam optimizer with an initial learning rate of 0.005 [22].

###  4.4 Objective Evaluation Methods

We evaluate our style transfer-based upmixing model against the baseline using
various objective metrics as follows:

  * •

Overall reconstruction quality: Scale-dependent source-to-distortion ratio
(SD-SDR) is employed to report the overall reconstruction quality by comparing
the time-domain reconstruction $\hat{\bm{z}}_{s}^{(5)}$ and the ground-truth
${\bm{z}_{s}^{(5)}}$ as SD-SDR is proven to better reflect the scale
reconstruction compared to ordinary SDR [23].

  * •

Spatial reconstruction quality: To validate the spatial reconstruction
quality, we propose a new metric using the Wasserstein distance on
interchannel level differences (WILD). Given the 5-channel magnitude
spectrogram $\bm{Z}\in\mathbb{R}_{+}^{5\times F\times T}$, we compute an ILD
matrix for each pair of channels $(i,j)$, e.g.,
$\bm{M}^{(i,j)}_{[f,t]}=20\log_{10}\big{(}\|\bm{Z}_{[i,f,t]}\|/\bm{\|Z}_{[j,f,t]}\|\big{)}$.
The five channels result in ten such ILD matrices in total, i.e.,
$\bm{M}\in\mathbb{R}^{10\times F\times T}$. The Wasserstein distance computes
on the histograms of target and estimation ILD matrices, i.e.,
$\sum_{k=1}^{10}\text{WILD}\big{(}\text{hist}(\bm{M}_{[k,:,:]}),\text{hist}(\hat{\bm{M}}_{[k,:,:]})\big{)}$.
Note that the Wasserstein distance improves the robustness of the comparison
when the two histograms are too different from each other, where KL-divergence
fails to quantify the dissimilarity.

  * •

Average angle difference: We use the least square algorithm to approximately
decompose each of the 5 channels into a linear combination of four
instrumental sources. Based on the relative amount of each instrument spread
in the 5 channels, we can estimate the virtual location of the instrument. We
report the difference between the ground-truth and estimated direction of the
upmixed instruments in terms of angles in degree.

###  4.5 Subjective Evaluations

We conduct an ABX test to subjectively evaluate the proposed spatial style
transfer algorithm. The ground-truth 5-channel version $\bm{Z}_{s}^{(5)}$ is
provided to the participants as the reference. Participants are asked to
choose the one that is more similar to $\bm{Z}_{s}^{(5)}$ between the style-
transferred reconstruction $\bm{\hat{Z}}_{s}^{(5)}$ and the baseline upmix
$\bm{\tilde{Z}}^{(5)}$, in terms of the incoming direction of the four sources
separately, as well as the overall spatial image. The test contains 5 trials
and is conducted in a professional surround-sound listening room.

##  5 Experimental Results

###  5.1 Analysis on the Learnt Latent Space

(a)

(b)

(c)

(d)

Fig. 3: Visualization of latent variable activations.

To validate the amount of information disentanglement, we analyze the
statistical properties of the learned latent variables $\bm{h}$, and whether
the latent variables $\bm{h}$ learn spatial information independently from
music content. In Fig. 3(a), we show the activation of the mean of the latent
distribution from two different panning configurations of the same song. The
two uncorrelated graphs indicate that the $J=50$ dimensions learned in
$\bm{h}$ differ from each other although the two versions share the exact same
stem tracks. On the contrary, Fig. 3(b) shows highly correlated activation of
two completely different songs, because their corresponding latent variables
succeed to learn the same representation based on their same panning
locations. These results suggest that the learned latent variables are
influenced more by the spatial arrangement of the stem tracks and not by music
content, achieving the desired disentanglement.

Another visualization of the latent variables further strengthens our claim.
We prepare a set of data samples, which are created from the combination of
five different songs, and five different spatial configurations. If the latent
space reflects the spatial information more than the music content, there must
be a latent structure that reflects the spatial characteristics more. In Fig.
3(c) and 3(d), we show a dimension-reduced latent space via principal
component analysis. In 3(c), the coloring is based on the songs, while in 3(d)
it is based on the spatial configurations. We can see that the segments from
the same spatial configuration form a cluster in 3(d), while segments that
belong to the same song are scattered everywhere. The average standard
deviation for each color group in 3(c) is relatively higher than the ones in
3(d), i.e., $0.060>0.038$. This observation again indicates that the encoder
extracts music-invariant spatial features successfully.

###  5.2 Style Transfer-Based Upmixing Results

We conduct the style transfer-based upmixing and compare the model’s
performance against the baseline model. To this end, we prepare 540 test
examples that are 10 seconds long and are rendered in one of the 20
arbitrarily defined panning configurations.

Table 1: Objective performance of the style transfer method.

---
| SD-SDR (dB) |  $\Delta$ Angle (∘) | WILD
Style-transfer | 8.71$\pm$1.50 | 19.5$\pm$15 | 58.02$\pm$13.78
Baseline |  4.20$\pm$0.71 |  49.75$\pm$19.75 |  70.84$\pm$10.03

  * •

Objective evaluations: Table 1 summarizes the objective evaluation results.
Our proposed style transfer-based upmixing method shows an obvious advantage
over the baseline upmixing model with a better reconstruction score, less
difference in angle, and WILD. It is noteworthy that the style-transferred
results bear higher SD-SDR than the fixed-panning-based baseline, even in the
presence of the machine learning model’s algorithmic artifact. The SD-SDR
score indicates that our model can properly fulfill the mission of arbitrarily
positioning instruments during 5-channel generation without significant harm
to the sound quality.

  * •

Subjective evaluations: Fig. 4 reports the percentage of subjects preferring
the proposed style transfer-based upmixing over the baseline, regarding the
direction similarity. Ten listening experts participated in our subjective
test. More than 80% of participants believe the proposed model is more similar
to the reference, in terms of the overall perceptual panning. When it comes to
the individual instruments, the reconstruction of vocals is the most favored,
which might account for the overall performance due to its dominance in pop
music. Our model gains preference from the listeners on drums and other,
although with higher variance on the latter one, while the advantage doesn’t
pass to bass. We believe that the different preference on different
instruments results from the fact that the decoder implicitly performs
instrument separation and relocation for estimated sources: when an instrument
is difficult to separate from the mixture, its relative scale in each channel
and the resulting perceived direction in the 5-channel may be negatively
affected, which also explains the phenomenon that our preference scores align
to the diverse instrument-specific separation performances on the MUSDB18
dataset in the literature [24].

###  5.3 Blind upmixing and other experiments

To test our model’s potential as a generator, we apply blind upmixing to
MUSDB18 signals and to real-world music signals from our internal dataset.
Because all the models are initially trained on MUSDB18, to let the model
better accommodate real-world music signals, we fine-tune the model with real-
world 5-channel music before blindly generating 5-channel examples. Due to the
fact that the blind upmixing process is absolutely random, we cannot validate
its quality and the surrounding effect via a quantitative measure. However, we
find that the upmixing results are convincing based on informal listening by
our peers. We provide some blind upmixing samples, both from MUSDB18 and the
real-world dataset at https://saige.sice.indiana.edu/research-
projects/generative-upmixing.

We evaluated different latent dimensions ranging from $J=20$ to $J=60$, while
there is no strong enough trend to recommend a single number. We eventually
chose $J=50$ for all the experiments.

Fig. 4: Subjective tests on the rendered multichannel samples.

##  6 Conclusion

In this paper, we proposed to tackle the stereo upmixing problem in a
generative way, where spatial images and music content in the 5-channel signal
are disentangled to allow style transfer-based upmixing. We formulated the
problem into a modified VAE model and train the latent space to capture music-
invariant spatial information, e.g., panning locations. Our experiments showed
that the learned latent variables successfully capture spatial information
separately from the music contents. Both the objective and subjective
evaluations demonstrate that our style transfer-based upmixing model achieved
interactive upmixing.

## References

  * [1] M. A. Gerzon,  “Optimal reproduction matrices for multispeaker stereo,”  in Audio Engineering Society Convention 91. Audio Engineering Society, 1991.
  * [2] A. Walther and C. Faller,  “Direct-ambient decomposition and upmix of surround signals,”  in 2011 IEEE Workshop on Applications of Signal Processing to Audio and Acoustics (WASPAA), 2011, pp. 277–280.
  * [3] C. Faller,  “Multiple-loudspeaker playback of stereo signals,”  Journal of the Audio Engineering Society, vol. 54, no. 11, pp. 1051–1064, 2006.
  * [4] S.-W. Jeon, Y.-C. Park, S.-P. Lee, and D.-H.Youn,  “Robust representation of spatial sound in stereo-to-multichannel upmix,”  in Audio Engineering Society Convention 128. Audio Engineering Society, 2010.
  * [5] C. Avendano and J. M. Jot,  “Ambience extraction and synthesis from stereo signals for multi-channel audio up-mix,”  in 2002 IEEE International Conference on Acoustics, Speech, and Signal Processing, 2002, vol. 2, pp. II–1957–II–1960.
  * [6] C. Avendano and J. M. Jot,  “Frequency domain techniques for stereo to multichannel upmix,”  Journal of the Audio Engineering Society, June 2002.
  * [7] J. Merimaa, M. M. Goodwin, and J.-M. Jot,  “Correlation-based ambience extraction from stereo recordings,”  Journal of the Audio Engineering Society, october 2007.
  * [8] Roy Irwan and Ronald M Aarts,  “Two-to-five channel sound processing,”  Journal of the Audio Engineering Society, vol. 50, no. 11, pp. 914–926, 2002.
  * [9] J. Usher and J. Benesty,  “Enhancement of spatial sound quality: A new reverberation-extraction audio upmixer,”  IEEE Transactions on Audio, Speech, and Language Processing, vol. 15, no. 7, pp. 2141–2150, 2007.
  * [10] Y.-H. Baek, S.-W. Jeon, Y.-C Park, and S. Lee,  “Efficient primary-ambient decomposition algorithm for audio upmix,”  in Audio Engineering Society Convention 133. Audio Engineering Society, 2012.
  * [11] K. M. Ibrahim and M. Allam,  “Primary-ambient extraction in audio signals using adaptive weighting and principal component analysis,”  in Proc. Sound and Music Computing Conference (SMC), 2016, pp. 227–232.
  * [12] C.Uhle and C.Paul,  “A supervised learning approach to ambience extraction from mono recordings for blind upmixing,”  in Proc. Int. Conf. Digital Audio Effects (DAFx), 2008.
  * [13] K. M. Ibrahim and M. Allam,  “Primary-ambient source separation for upmixing to surround sound systems,”  in 2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP). IEEE, 2018, pp. 431–435.
  * [14] J. Choi and J. H. Chang,  “Exploiting deep neural networks for two-to-five channel surround decoder,”  Journal of the Audio Engineering Society, vol. 68, no. 12, pp. 938–949, 2021.
  * [15] S. Kraft and U. Zölzer,  “Stereo signal separation and upmixing by mid-side decomposition in the frequency-domain,”  in 18th International Conference on Digital Audio Effects (DAFx), 2015.
  * [16] S.Y. Park, C. J. Chun, and H. K. Kim,  “Subband-based upmixing of stereo to 5.1-channel audio signals using deep neural networks,”  in 2016 International Conference on Information and Communication Technology Convergence (ICTC), 2016.
  * [17] D. P. Kingma and M. Welling,  “Auto-encoding variational bayes,”  in Proceedings of the International Conference on Learning Representations (ICLR), 2014.
  * [18] G. Huang, Z. Liu, L. Van Der Maaten, and K. Q. Weinberger,  “Densely connected convolutional networks,”  in Proceedings of the IEEE International Conference on Computer Vision and Pattern Recognition (CVPR), 2017, pp. 4700–4708.
  * [19] V. Pulkki,  “Virtual sound source positioning using vector base amplitude panning,”  Journal of the audio engineering society, vol. 45, no. 6, pp. 456–466, 1997.
  * [20] “Multichannel stereophonic sound system with and without accompanying picture,”  Standard, International Telecommunication Union, Aug. 2012.
  * [21] Z. Rafii, A. Liutkus, F.-R. Stöter, S. I. Mimilakis, and R. Bittner,  “MUSDB18-HQ - an uncompressed version of MUSDB18,” Aug. 2019.
  * [22] D.P. Kingma and J. Ba,  “Adam: A method for stochastic optimization,”  in Proceedings of the International Conference on Learning Representations (ICLR), 2015.
  * [23] J. Le Roux, S. Wisdom, H. Erdogan, and J. R. Hershey,  “SDR – half-baked or well done?,”  in IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2019.
  * [24] D. Stoller, S. Ewert, and S. Dixon,  “Wave-u-net: A multi-scale neural network for end-to-end audio source separation,”  in Proceedings of the International Society for Music Information Retrieval Conference (ISMIR), 2018.


#  Thermal conductivity of strongly coupled Yukawa fluids

Sergey A. Khrapak<EMAIL_ADDRESS>Joint Institute for High
Temperatures, Russian Academy of Sciences, 125412 Moscow, Russia

###### Abstract

A vibrational model of heat conduction in liquids with soft pairwise
interactions is applied to estimate the thermal conductivity coefficient of
strongly coupled Yukawa fluids. A reasonable agreement with the available data
from numerical simulations is observed. The results can be useful in the
context of strongly coupled plasma and complex (dusty) plasma fluids, when
Yukawa (or screened Coulomb) interaction potential is applicable.

A simple vibrational model of heat transfer in simple liquids with soft
pairwise interatomic interactions has been put forward recently Khrapak
(2021a). In this model the thermal conductivity coefficient is expressed via
an effective frequency of atomic vibrations around a local temporary
equilibrium position, mean interatomic separation and specific heat capacity.
The vibrational frequency is determined by averaging over the liquid
collective mode excitation spectrum. The model has been applied to quantify
heat transfer in a dense Lennard-Jones liquid and a strongly coupled one-
component plasma (OCP). A remarkable agreement with the available numerical
results has been documented Khrapak (2021a). The purpose of this Brief
Communication is to apply the model to strongly coupled Yukawa fluids and to
examine its accuracy in this case.

In the vibrational model of thermal conduction, a liquid is approximated by a
layered structure with layers perpendicular to the temperature gradient and
separated by a distance $\Delta=n^{-1/3}$, where $n$ is the particle (or
atomic) density. The average interparticle separation within such quasi-layers
is also $\Delta$. The particle positions in the liquid’s quasi-layers are not
regular and are not fixed. The particles can migrate within a given layer as
well as between different layers. However, the time scale of these migrations
is relatively long in dense liquids, much longer than the characteristic time
of latticle-like vibrations of a particle around its temporary equilibrium
position Frenkel (1955). This resembles the applicability condition of the
quasi-localized charge approximation (QLCA), which is very successively used
to describe collective modes in strongly coupled plasma-related systems Golden
and Kalman (2000).

The average difference in energy between the atoms of adjacent layers is
$\Delta(dU/dx)=\Delta c_{\rm v}(dT/dx)$, where $U$ is the internal energy and
$c_{\rm v}=dU/dT$ is the specific heat at constant volume. Although the
specific heat at constant pressure, $c_{\rm p}$, would be apparently more
appropriate for such an estimate, the difference between $c_{\rm p}$ and
$c_{\rm v}$ is vanishingly small in strongly coupled Yukawa fluids. At the
same time $c_{\rm v}$ is easier to evaluate (see below) so it is appropriate
to use $c_{\rm v}$ for practical purposes. The energy exchange between
neighbouring quasi-layers occurs at a characteristic frequency $\nu$, which is
set equal to the average vibrational frequency of a particle,
$\nu=\langle\omega\rangle/2\pi$. The energy flux between two layers per unit
area is ${\mathcal{E}}=-\nu\Delta c_{\rm v}(dT/dx)/\Delta^{2}$. This should be
compared with the Fourier law, ${\mathcal{E}}=-\lambda(dT/dx)$, resulting in

$\lambda=c_{\rm v}\frac{\langle\omega\rangle}{2\pi\Delta}.$ (1)

This is a simple general expression without free parameters, which we are
going to examine in detail using strongly coupled Yukawa fluids.

Since the actual frequency distribution can be quite complex in liquids, and
can vary from one type of liquid to another, some simplifying assumptions have
to be employed in evaluating $\langle\omega\rangle$. Before we do that for
Yukawa fluids, let us point out two special cases of Eq. (1). In the simplest
Einstein approximation all atoms vibrate with the same (Einstein) frequency
$\Omega_{\rm E}$ so that $\langle\omega\rangle=\Omega_{\rm E}$. We obtain the
expression proposed by Horrocks and McLaughlin (to within some difference in
numerical coefficients) Horrocks and McLaughlin (1960)

$\lambda=c_{\rm v}\frac{\Omega_{\rm E}}{2\pi\Delta}.$ (2)

As a more involved approximation, assume that a dense liquid supports one
longitudinal and two transverse modes, which exhibit acoustic dispersion
relations in the long-wavelength limit. The Debye-like averaging procedure
yields Khrapak (2021a)

$\lambda\simeq\frac{1}{4}\left(\frac{3}{4\pi}\right)^{1/3}c_{\rm
v}\frac{c_{l}+2c_{t}}{\Delta^{2}},$ (3)

where $c_{l}$ and $c_{t}$ are the longitudinal and transverse sound
velocities, respectively. This is similar (again to within numerical
coefficients) to a formula proposed by Cahill and Pohl to describe thermal
conductivity coefficient of amorphous inorganic solids Cahill and Pohl (1989);
Cahill  _et al._ (1992).

The Yukawa systems considered here represent a collection of point-like
charged particles immersed into a neutralizing medium (usually conventional
electron-ion plasma), which provides screening. The pairwise Yukawa repulsive
interaction potential (also known as screened Coulomb or Debye-Hückel
potential) is

$\phi(r)=(Q^{2}/r)\exp(-\kappa r/a),$ (4)

where $Q$ is the particle charge and $\kappa$ is the dimensionless screening
parameter, which is the ratio of the Wigner-Seitz radius $a=(4\pi n/3)^{-1/3}$
to the plasma screening length. Yukawa potential is widely used as a
reasonable first approximation for actual interactions in three-dimensional
isotropic complex plasmas and colloidal suspensions Tsytovich (1997); Fortov
_et al._ (2004, 2005); Khrapak  _et al._ (2008); Khrapak and Morfill (2009);
Chaudhuri  _et al._ (2011); Ivlev  _et al._ (2012); Lampe and Joyce (2015).

The dynamics and thermodynamics of Yukawa systems are conventionally
characterized by the screening parameter $\kappa$ and the coupling parameter
$\Gamma=Q^{2}/aT$, where $T$ is the system temperature (in energy units). The
coupling parameter characterizes the ratio between the potential energy of
interparticle interaction and the kinetic energy. In strongly coupled Yukawa
fluids the condition $\Gamma\gg 1$ should be satisfied. For even higher
coupling Yukawa fluids crystallize, forming either the body-centered-cubic
(bcc) or face-centered-cubic (fcc) lattices. Detailed phase diagrams of Yukawa
systems are available in the literature Robbins  _et al._ (1988); Hamaguchi
_et al._ (1996, 1997); Vaulina and Khrapak (2000); Vaulina  _et al._ (2002).
The screening parameter $\kappa$ determines the softness of the interparticle
repulsion. It varies from the extremely soft and long-ranged Coulomb potential
at $\kappa\rightarrow 0$ (OCP limit) to the hard-sphere-like interaction limit
at $\kappa\rightarrow\infty$. In the context of complex plasmas and colloidal
suspensions the relatively “soft” regime, $\kappa\sim{\mathcal{O}}(1)$, is of
particular interest. Most of previous investigations have focused on this
regime and we follow this tradition.

Transport properties (in particular, self-diffusion, shear viscosity, and
thermal conductivity coefficients) of Yukawa fluids have been relatively well
investigated. Thermal conduction in three-dimensional Yukawa fluids has been
studied and discussed in Refs. Salin and Caillol (2002, 2003); Faussurier and
Murillo (2003); Donkó and Hartmann (2004); Vaulina  _et al._ (2010); Shahzad
and He (2012); Ott  _et al._ (2015); Scheiner and Baalrud (2019); Kählert
(2020). For a detailed comparison we chose the data from non-equilibrium
molecular dynamics simulations reported by Donko and Hartmann Donkó and
Hartmann (2004). This is the most complete set of data in terms of covered
$\Gamma$ and $\kappa$ parameters from those available presently (note that
there exists some deviations between the results from different investigations
as shown for instance in Fig. 5 of Ref. Kählert (2020)). Throughout this paper
we employ Rosenfeld’s normalization Rosenfeld (1999) for the thermal
conductivity coefficient $\lambda_{\rm R}=\lambda n^{-2/3}/v_{\rm T}$ (where
$v_{\rm T}=\sqrt{T/m}$ is the thermal velocity), which is different from that
used originally, $\lambda^{\prime}=\lambda/(n\omega_{p}a^{2})$, where
$\omega_{\rm p}=\sqrt{4\pi Q^{2}n/m}$ is the plasma frequency scale. Note that
since we express temperatures in energy units, the Boltzmann constant does not
show up in the expressions (effectively $k_{\rm B}=1$).

Figure 1: (Color online) Reduced thermal conductivity coefficient
$\lambda_{\rm R}$ versus the coupling parameter $\Gamma$. Citcles, triangles,
and rhombs correspond to molecular dynamics results from Ref. Donkó and
Hartmann (2004) for $\kappa=1$, $\kappa=2$, and $\kappa=3$, respectively.
Crosses are the OCP numerical results from Ref. Scheiner and Baalrud (2019).
The solid curves are the theoretical calculations using Eq. (1). The dotted
curves (shown for $\kappa=2$ and $\kappa=3$) are plotted using Eq. (3). The
results of calculation using Eq. (2) are not shown for clarity, but are
discussed in the text instead.

The numerical data for the dependence of $\lambda_{\rm R}$ on $\Gamma$ are
plotted in Fig. 1. The solid symbols correspond to the data from a non-
equilibrium simulation in Ref. Donkó and Hartmann (2004) for strongly coupled
Yukawa fluids with $\kappa=1.0$, $2.0$, and $3.0$. The crosses correspond to
the OCP results from Ref. Scheiner and Baalrud (2019), added for completeness.
The minima in $\lambda_{\rm R}$ separate the regions of gaseous-like and
fluid-like regimes of thermal conductivity. We will consider only the fluid-
like (i.e. strongly coupled) regime to the right from the minima.

Figure 2: (Color online) Dispersion relations of strongly coupled Yukawa
fluids with $\kappa=1.0$ (a) and $\kappa=3.0$ (b). The wave frequency $\omega$
is expressed in units of the plasma frequency $\omega_{\rm p}$, the reduced
wave-number is $q=ka$. The solid curves correspond to the simplified QLCA
dispersion relations of Eq. (6) for the longitudinal mode and Eq. (7) for the
transverse mode. Dashes lines are the corresponding acoustic asymptotes. They
terminate at $q_{\rm max}=(9\pi/2)^{1/3}$, the maximum cutoff wavelength
involved in the averaging procedure. Note different scales of the vertical
axis in (a) and (b).

To get the average vibrational frequency, we perform the conventional
averaging over collective modes that are supported in a strongly coupled
(dense) fluid. Taking into account one longitudinal and two transverse modes
the average frequency becomes Khrapak (2021a)

$\langle\omega\rangle=\frac{2}{9\pi}\int_{0}^{\rm q_{\rm
max}}q^{2}dq\left[\omega_{l}(q)+2\omega_{t}(q)\right],$ (5)

where $q=ka$ is the reduced wave number and $\omega_{l,t}(q)$ are the
dispersion relations of the longitudinal and transverse modes. The cutoff wave
number $q_{\rm max}$ appears due to the normalization condition, which should
ensure that each collective mode contains $n$ oscillations. Mathematically
this implies that

$\frac{4\pi}{3}\left(\frac{k_{\rm
max}}{2\pi}\right)^{3}=n;\quad\Rightarrow\quad q_{\rm
max}=\left(\frac{9\pi}{2}\right)^{1/3}\simeq 2.418.$

This cutoff also ensures that averaging of a quantity that does not depend on
$q$ will not change its actual value. The longitudinal and transverse
collective mode frequencies can be well approximated within the QLCA approach,
which relates the dispersion relations to the pairwise interaction potential
and structural characteristics, expressed in the form of the radial
distribution function (RDF) $g(r)$ Golden and Kalman (2000). Additionally, we
make use of a simple excluded cavity model for the fluid RDF Khrapak  _et al._
(2016). This results in simple and elegant expressions for the dispersion
relations of the longitudinal and transverse modes Khrapak  _et al._ (2016):

$\displaystyle\omega_{l}^{2}=\omega_{\rm
p}^{2}e^{-R\kappa}\left[\left(1+R\kappa\right)\left(\frac{1}{3}-\frac{2\cos
Rq}{R^{2}q^{2}}+\frac{2\sin Rq}{R^{3}q^{3}}\right)\right.$ (6)
$\displaystyle\left.-\frac{\kappa^{2}}{\kappa^{2}+q^{2}}\left(\cos
Rq+\frac{\kappa}{q}\sin Rq\right)\right],$

and

$\omega_{t}^{2}=\omega_{\rm
p}^{2}e^{-R\kappa}\left(1+R\kappa\right)\left(\frac{1}{3}+\frac{\cos
Rq}{R^{2}q^{2}}-\frac{\sin Rq}{R^{3}q^{3}}\right).$ (7)

Here the excluded cavity [the region where $g(r)=0$] radius $R$ is expressed
in units of $a$ and is related to the screening parameter via Khrapak (2017);
Khrapak and Khrapak (2018)

$R(\kappa)\simeq
1+\frac{1}{\kappa}\ln\left[\frac{3\cosh(\kappa)}{\kappa^{2}}-\frac{3\sinh(\kappa)}{\kappa^{3}}\right].$
(8)

These fully analytical dispersion relations of strongly coupled Yukawa fluids
demonstrate a very high accuracy in the long-wavelength regime Khrapak  _et
al._ (2016); Khrapak (2017); Khrapak and Khrapak (2018). The long-wavelength
portion of these dispersion relations for $\kappa=1$ and $\kappa=3$ are
plotted in Fig. 2. Note that although the existence of a so-called “$q$-gap”
(zero frequency non-propagating branch at long wavelengths Murillo (2000);
Goree  _et al._ (2012); Trachenko and Brazhkin (2015); Khrapak  _et al._
(2019); Kryuchkov  _et al._ (2019)) in the transverse mode is not reproduced,
this does not affect our main results significantly, because the contribution
from the long wavelengths (low frequencies) to the thermal conductivity is
relatively small.

We can now compare simulations results with theoretical approximations of Eqs.
(1), (2), and (3).

The integral in Eq. (5) has been evaluated numerically for different screening
parameters and the resulting values for $\langle\omega\rangle/\omega_{\rm p}$
are presented in Table 1. In the same table the values of $\Omega_{\rm
E}/\omega_{\rm p}$ needed to evaluate Eq. (2) are also given. Here we use the
fact that the ratio $\Omega_{\rm E}/\omega_{\rm p}$ is virtually independent
of $\Gamma$ in the strongly coupled regime (be it fluid or crystal) Khrapak
and Klumov (2018) and use the fcc data tabulated by Ohta and Hamaguchi Ohta
and Hamaguchi (2000). Note that the dependence on the coupling parameter is
absorbed in the plasma frequency. Also quoted in the Table 1 are the values
$\sqrt{\langle\omega^{2}\rangle/\omega_{\rm p}^{2}}$ evaluated inserting
$\omega_{l}^{2}(q)+2\omega_{t}^{2}(q)$ under the integral in Eq. (5). Very
close agreement with the values $\Omega_{\rm E}/\omega_{\rm p}$ demonstrates
that all the approximations made to deal with dispersion relations of strongly
coupled Yukawa fluids are appropriate.

The longitudinal and transverse sound velocities can be evaluated from

$c_{l,t}=\lim_{q\rightarrow 0}\frac{\omega_{l,t}(q)a}{q}.$ (9)

Simple expressions are available Khrapak and Khrapak (2018), and these were
used to evaluate $c_{l}$ and $c_{t}$ to be substituted in Eq. (3).

Table 1:  The average of vibrational frequencies $\langle\omega\rangle$ and $\sqrt{\langle\omega^{2}\rangle}$ of strongly coupled Yukawa fluids as well as the Einstein frequency $\Omega_{\rm E}$ of the fcc Yukawa solid for screening parameters up to $\kappa=5$. All frequencies are reduced using the plasma particle frequency. $\kappa$ | 1.0 | 2.0 | 3.0 | 4.0 | 5.0
---|---|---|---|---|---
$\langle\omega\rangle/\omega_{\rm p}$ | 0.420 | 0.279 | 0.164 | 0.088 | 0.044
$\sqrt{\langle\omega^{2}\rangle}/\omega_{\rm p}$ | 0.469 | 0.305 | 0.176 | 0.094 | 0.0469
$\Omega_{\rm E}/\omega_{\rm p}$ | 0.472 | 0.307 | 0.176 | – | 0.0469

The last remaining step is to evaluate the specific heat $c_{\rm v}$.
Thermodynamics of Yukawa fluids has been relatively well investigated
Hamaguchi  _et al._ (1997); Faussurier and Murillo (2003); Tolias  _et al._
(2014); Khrapak and Thomas (2015); Khrapak (2015); Khrapak  _et al._ (2015);
Tolias  _et al._ (2015); Tolias and Castello (2019). Various practical
expressions for thermodynamic functions have been proposed in the literature.
Among these, the freezing temperature scaling of the thermal component of the
internal energy proposed by Rosenfeld and Tarazona Rosenfeld and Tarazona
(1998); Rosenfeld (2000) combines relative good accuracy with simplicity and
relatively wide applicability Ingebrigtsen  _et al._ (2013). In this
apprximation the thermal correction to the internal energy of Yukawa fluids
(the latter is dominated by the static contribution for sufficiently soft
potentials, like Yukawa considered here) scales as $\propto T^{3/5}$ and hence
the reduced correction is $\propto T^{-2/5}\propto\Gamma^{2/5}$. In the range
of $\kappa$ considered the reduced thermal correction per particle can be well
approximated as $u_{\rm th}\simeq 3.1(\Gamma/\Gamma_{\rm fr})^{2/5}$ Khrapak
_et al._ (2015), where $\Gamma_{\rm fr}$ is the value of the coupling
parameter at freezing, which is very close to that at melting, $\Gamma_{\rm
m}$ (the coexistence region is very narrow, so that usually there is no need
to distinguish between $\Gamma_{\rm fr}$ and $\Gamma_{\rm m}$ for Yukawa
systems). The corresponding specific heat at constant volume is therefore

$c_{\rm v}\simeq 1.5+1.86\left(\frac{\Gamma}{\Gamma_{\rm fr}}\right)^{2/5}.$
(10)

The values of $\Gamma_{\rm fr}$ for various $\kappa$ have been obtained and
tabulated in Ref. Hamaguchi  _et al._ (1997) (relatively accurate fits are
also available Vaulina and Khrapak (2000); Vaulina  _et al._ (2002).) We now
have all the necessary information to evaluate each of the equations (1), (2),
and (3).

A comparison between the vibrational model of heat conduction and the results
from numerical experiments is shown in Fig. 1. The solid curves are plotted
using the general expression (1), supplemented by the average vibrational
frequencies from Table 1 and the specific heat from Eq. (10). There is a
reasonable agreement with numerical results in particular in the weak
screening regime with smaller $\kappa$. In the OCP case we have performed
averaging with the help of Eqs. (6) and (7) with $\kappa=0$ Khrapak  _et al._
(2016), but used the accurate OCP equation of state to evaluate $c_{\rm v}$
(for further details see Ref. Khrapak (2021a)). The agreement is excellent.
Note that in the near-OCP regime with $\kappa\lesssim 1$ the thermal
conductivity coefficient exhibits a very weak dependence on $\kappa$ (thus OCP
result can serve as a relevant approximation in this regime).

The dotted curves have been calculated using Eq. (3). This approximation leads
to some overprediction of the thermal conductivity coefficient. The origin of
this overprediction is immediately clear from Fig. 2. The acoustic asymptotes
are overestimating the actual frequencies, in particular for the longitudinal
collective mode. The situation gets particularly worse at weak screening. Here
the longitudinal sound velocity diverges as $\propto\kappa^{-1}$ Khrapak
(2019) (in the OCP limit the longitudinal mode exhibits the plasmon dispersion
so that formally $c_{l}\rightarrow\infty$). This produces grossly inaccurate
results at weak screening and we do not show the results of calculation using
Eq. (3) for $\kappa=1$. For $\kappa=2$ this approximation somewhat
overestimate the numerical data and for $\kappa=3$ it is located rather close
to the numerical results. Thus, the averaging over a Debye-like acoustic
spectrum is inappropriate in the regime of soft long-ranged potential, but
becomes more appropriate when the potential becomes steeper (see Fig. 2 for
illustration). It is not surprising, therefore, that this approximation works
rather well in the case of Lennard-Jones liquids with a rlatively steep
interaction potential Khrapak (2021a).

The calculation using Eq. (2) is not shown in Fig. 2 for clarity, but actually
this is a very good approximation for the regime considered. From the Table 1
we see that the values $\Omega_{\rm E}$ are somewhat higher that those of
$\langle\omega\rangle$. This implies that for finite $\kappa$ this
approximation would be closer to the numerical results than the general
expression (1). In the OCP limit it only slightly overestimates the numerical
results, as has been demonstrated previously Khrapak (2021a).

The main conclusion can be formulated as follows. The vibrational model of
heat transfer works reasonably well for strongly coupled Yukawa fluids. This
model relates the coefficient of thermal conductivity to the specific heat,
mean inter-atomic separation and an average vibrational frequency of atoms
around their temporary equilibrium positions. The general formula, which uses
averaging over collective modes of dense liquids, and its special case, which
assumes that all atoms vibrate with the same Einstein frequency, are
reasonably accurate in the entire regime investigated. The Debye-like
averaging, which uses long-wavelength acoustic asymptotes, is not appropriate
near the OCP limit, but becomes appropriate as the steepness of the
interaction potential (regulated by the screening parameter $\kappa$)
increases.

The model can be straightforwardly generalized to heat transfer in two-
dimensional (2D) Yukawa layers and this has been recently discussed in Ref.
Khrapak (2021b). In the context of 2D complex (dusty) plasma fluids and solids
in ground based laboratories, the regime with $\kappa\lesssim 1$ is of
particular practical interest. Based on the results of the present analysis it
is expected that estimations for the (2D) OCP limit should be relevant in this
regime and the Einstein model should be appropriate as well.

Data sharing is not applicable to this article as no new data were created.

## References

  * Khrapak (2021a) S. A. Khrapak, “Vibrational model of thermal conduction for fluids with soft interactions,” Phys. Rev. E 103, 013207 (2021a).
  * Frenkel (1955) Y. Frenkel,  _Kinetic theory of liquids_ (Dover, New York, NY, 1955).
  * Golden and Kalman (2000) K. I. Golden and G. J. Kalman, “Quasilocalized charge approximation in strongly coupled plasma physics,” Phys. Plasmas 7, 14–32 (2000).
  * Horrocks and McLaughlin (1960) J. K. Horrocks and E. McLaughlin, “Thermal conductivity of simple molecules in the condensed state,” Trans. Faraday Soc. 56, 206 (1960).
  * Cahill and Pohl (1989) D. G. Cahill and R.O. Pohl, “Heat flow and lattice vibrations in glasses,” Solid State Commun. 70, 927–930 (1989).
  * Cahill  _et al._ (1992) D. G. Cahill, S. K. Watson, and R. O. Pohl, “Lower limit to the thermal conductivity of disordered crystals,” Phys. Rev. B 46, 6131–6140 (1992).
  * Tsytovich (1997) V. Tsytovich, “Dust plasma crystals, drops, and clouds.” Phys.-Usp. 40, 53–94 (1997).
  * Fortov  _et al._ (2004) V. E. Fortov, A. G. Khrapak, S. A. Khrapak, V. I. Molotkov, and O. F. Petrov, “Dusty plasmas,” Phys.-Usp. 47, 447 – 492 (2004).
  * Fortov  _et al._ (2005) V. E. Fortov, A. V. Ivlev, S. A. Khrapak, A. G. Khrapak, and G. E. Morfill, “Complex (dusty) plasmas: Current status, open issues, perspectives,” Phys. Rep. 421, 1–103 (2005).
  * Khrapak  _et al._ (2008) S. A. Khrapak, B. A. Klumov, and G. E. Morfill, “Electric potential around an absorbing body in plasmas: Effect of ion-neutral collisions,” Phys. Rev. Lett. 100, 225003 (2008).
  * Khrapak and Morfill (2009) S. Khrapak and G. Morfill, “Basic processes in complex (dusty) plasmas: Charging, interactions, and ion drag force,” Contrib. Plasma Phys. 49, 148–168 (2009).
  * Chaudhuri  _et al._ (2011) M. Chaudhuri, A. V. Ivlev, S. A. Khrapak, H. M. Thomas, and G. E. Morfill, “Complex plasma—the plasma state of soft matter,” Soft Matter 7, 1287–1298 (2011).
  * Ivlev  _et al._ (2012) A. Ivlev, H. Löwen, G. Morfill, and C. P. Royall,  _Complex Plasmas and Colloidal Dispersions: Particle-Resolved Studies of Classical Liquids and Solids_ (World Scientific, 2012).
  * Lampe and Joyce (2015) M. Lampe and G. Joyce, “Grain-grain interaction in stationary dusty plasma,” Phys. Plasmas 22, 023704 (2015).
  * Robbins  _et al._ (1988) M. O. Robbins, K. Kremer, and G. S. Grest, “Phase diagram and dynamics of yukawa systems,” J. Chem. Phys. 88, 3286–3312 (1988).
  * Hamaguchi  _et al._ (1996) S. Hamaguchi, R. T. Farouki, and D. H. E. Dubin, “Phase diagram of Yukawa systems near the one-component-plasma limit revisited,” J. Chem. Phys. 105, 7641–7647 (1996).
  * Hamaguchi  _et al._ (1997) S. Hamaguchi, R. T. Farouki, and D. H. E. Dubin, “Triple point of yukawa systems,” Phys. Rev. E 56, 4671–4682 (1997).
  * Vaulina and Khrapak (2000) O. S. Vaulina and S. A. Khrapak, “Scaling law for the fluid-solid phase transition in yukawa systems (dusty plasmas),” JETP 90, 287–289 (2000).
  * Vaulina  _et al._ (2002) O. Vaulina, S. Khrapak, and G. Morfill, “Universal scaling in complex (dusty) plasmas,” Phys. Rev. E 66, 016404 (2002).
  * Salin and Caillol (2002) G. Salin and J.-M. Caillol, “Transport coefficients of the yukawa one-component plasma,” Phys. Rev. Lett. 88, 065002 (2002).
  * Salin and Caillol (2003) G. Salin and J.-M. Caillol, “Equilibrium molecular dynamics simulations of the transport coefficients of the yukawa one component plasma,” Phys. Plasmas 10, 1220–1230 (2003).
  * Faussurier and Murillo (2003) G. Faussurier and M. S. Murillo, “Gibbs-bogolyubov inequality and transport properties for strongly coupled yukawa fluids,” Phys. Rev. E 67, 046404 (2003).
  * Donkó and Hartmann (2004) Z. Donkó and P. Hartmann, “Thermal conductivity of strongly coupled yukawa liquids,” Phys. Rev. E 69, 016405 (2004).
  * Vaulina  _et al._ (2010) O. S. Vaulina, X. G. Koss, Yu. V. Khrustalyov, O. F. Petrov, and V. E. Fortov, “Thermodynamic and transport properties of nonideal systems with isotropic pair potentials,” Phys. Rev. E 82, 056411 (2010).
  * Shahzad and He (2012) A. Shahzad and M.-G. He, “Thermal conductivity calculation of complex (dusty) plasmas,” Phys. Plasmas 19, 083707 (2012).
  * Ott  _et al._ (2015) T. Ott, M. Bonitz, and Z. Donkó, “Effect of correlations on heat transport in a magnetized strongly coupled plasma,” Phys. Rev. E 92, 063105 (2015).
  * Scheiner and Baalrud (2019) B. Scheiner and S. D. Baalrud, “Testing thermal conductivity models with equilibrium molecular dynamics simulations of the one-component plasma,” Phys. Rev. E 100, 043206 (2019).
  * Kählert (2020) H. Kählert, “Thermodynamic and transport coefficients from the dynamic structure factor of yukawa liquids,” Phys. Rev. Research 2, 033287 (2020).
  * Rosenfeld (1999) Y. Rosenfeld, “A quasi-universal scaling law for atomic transport in simple fluids,” J. Phys.: Condens. Matter 11, 5415–5427 (1999).
  * Khrapak  _et al._ (2016) S. A. Khrapak, B. Klumov, L. Couedel, and H. M. Thomas, “On the long-waves dispersion in yukawa systems,” Phys. Plasmas 23, 023702 (2016).
  * Khrapak (2017) S. A. Khrapak, “Practical dispersion relations for strongly coupled plasma fluids,” AIP Adv. 7, 125026 (2017).
  * Khrapak and Khrapak (2018) S. Khrapak and A. Khrapak, “Simple dispersion relations for coulomb and yukawa fluids,” IEEE Trans. Plasma Sci. 46, 737–742 (2018).
  * Murillo (2000) M. S. Murillo, “Critical wave vectors for transverse modes in strongly coupled dusty plasmas,” Phys. Rev. Lett. 85, 2514–2517 (2000).
  * Goree  _et al._ (2012) J. Goree, Z. Donkó, and P. Hartmann, “Cutoff wave number for shear waves and maxwell relaxation time in yukawa liquids,” Phys. Rev. E 85, 066401 (2012).
  * Trachenko and Brazhkin (2015) K. Trachenko and V. V. Brazhkin, “Collective modes and thermodynamics of the liquid state,” Rep. Progr. Phys. 79, 016502 (2015).
  * Khrapak  _et al._ (2019) S. A. Khrapak, A. G. Khrapak, N. P. Kryuchkov, and S. O. Yurchenko, “Onset of transverse (shear) waves in strongly-coupled yukawa fluids,” J. Chem. Phys. 150, 104503 (2019).
  * Kryuchkov  _et al._ (2019) N. P. Kryuchkov, L. A. Mistryukova, V. V. Brazhkin, and S. O. Yurchenko, “Excitation spectra in fluids: How to analyze them properly,” Sci. Rep. 9, 10483 (2019).
  * Khrapak and Klumov (2018) S. Khrapak and B. Klumov, “High-frequency elastic moduli of two-dimensional yukawa fluids and solids,” Phys. Plasmas 25, 033706 (2018).
  * Ohta and Hamaguchi (2000) H. Ohta and S. Hamaguchi, “Molecular dynamics evaluation of self-diffusion in yukawa systems,” Phys. Plasmas 7, 4506–4514 (2000).
  * Tolias  _et al._ (2014) P. Tolias, S. Ratynskaia, and U. de Angelis, “Soft mean spherical approximation for dusty plasma liquids: One-component yukawa systems with plasma shielding,” Phys. Rev. E 90, 053101 (2014).
  * Khrapak and Thomas (2015) S. A. Khrapak and H. M. Thomas, “Practical expressions for the internal energy and pressure of yukawa fluids,” Phys. Rev. E 91, 023108 (2015).
  * Khrapak (2015) S. A. Khrapak, “Thermodynamics of yukawa systems and sound velocity in dusty plasmas,” Plasma Phys. Control. Fusion 58, 014022 (2015).
  * Khrapak  _et al._ (2015) S. A. Khrapak, N. P. Kryuchkov, S. O. Yurchenko, and H. M. Thomas, “Practical thermodynamics of yukawa systems at strong coupling,” J. Chem. Phys. 142, 194903 (2015).
  * Tolias  _et al._ (2015) P. Tolias, S. Ratynskaia, and U. de Angelis, “Soft mean spherical approximation for dusty plasma liquids: Level of accuracy and analytic expressions,” Phys. Plasmas 22, 083703 (2015).
  * Tolias and Castello (2019) P. Tolias and F. Lucco Castello, “Isomorph-based empirically modified hypernetted-chain approach for strongly coupled yukawa one-component plasmas,” Phys. of Plasmas 26, 043703 (2019).
  * Rosenfeld and Tarazona (1998) Y. Rosenfeld and P. Tarazona, “Density functional theory and the asymptotic high density expansion of the free energy of classical solids and fluids,” Mol. Phys. 95, 141–150 (1998).
  * Rosenfeld (2000) Y. Rosenfeld, “Excess-entropy and freezing-temperature scalings for transport coefficients: Self-diffusion in yukawa systems,” Phys. Rev. E 62, 7524–7527 (2000).
  * Ingebrigtsen  _et al._ (2013) T. S. Ingebrigtsen, A. A. Veldhorst, T. B. Schrøder, and J. C. Dyre, “Communication: The rosenfeld-tarazona expression for liquids’ specific heat: A numerical investigation of eighteen systems,” J. Chem. Phys. 139, 171101 (2013).
  * Khrapak (2019) S. A. Khrapak, “Unified description of sound velocities in strongly coupled yukawa systems of different spatial dimensionality,” Phys. Plasmas 26, 103703 (2019).
  * Khrapak (2021b) S. A. Khrapak, “Thermal conduction in two-dimensional complex plasma layers,” Phys. Plasmas 28, 010704 (2021b).