# Tunable multiwindow magnomechanically induced transparency, Fano resonances,
and slow to fast light conversion

Kamran Ullah<EMAIL_ADDRESS>M. Tahir Naseem<EMAIL_ADDRESS>Özgür E. Müstecaplıoğlu<EMAIL_ADDRESS>Department of Physics, Koç
University, Sarıyer, İstanbul, 34450, Turkey

###### Abstract

We investigate the absorption and transmission properties of a weak probe
field under the influence of a strong control field in a cavity
magnomechanical system. The system consists of two ferromagnetic-material
yttrium iron garnet (YIG) spheres coupled to a single cavity mode. In addition
to two magnon-induced transparencies (MITs) that arise due to magnon-photon
interactions, we observe a magnomechanically induced transparency (MMIT) due
to the presence of nonlinear magnon-phonon interaction. We discuss the
emergence of Fano resonances and explain the splitting of a single Fano
profile to double and triple Fano profiles due to additional couplings in the
proposed system. Moreover, by considering a two-YIG system, the group delay of
the probe field can be enhanced by one order of magnitude as compared with a
single-YIG magnomechanical system. Furthermore, we show that the group delay
depends on the tunability of the coupling strength of the first YIG with
respect to the coupling frequency of the second YIG, and vice versa. This
helps to achieve larger group delays for weak magnon-photon coupling
strengths.

Magnon induced transparency; magnomechanical induced transparency; Fano
resonances; subluminal and superluminal effects.

##  I introduction

Storing information in different frequency modes of light has attracted much
attention due to its critical role in high-speed, long-distance quantum
communication applications Afzelius  _et al._ (2009); Sinclair  _et al._
(2014); Jobez  _et al._ (2016). The spectral distinction of optical signals
eliminates their unintentional coupling to the stationary information or
memory nodes in a communication network. For that aim, multiple transparency
window Electromagnetically Induced Transparency (EIT) schemes have been
considered for multiband quantum memory implementations mainly in the medium
of three-level cold atoms. Experimental demonstrations of three EIT windows
have been reported McGloin  _et al._ (2000), and extended to seven windows
using external fields Ying  _et al._ (2014). Observation of nine EIT windows
has been experimentally demonstrated quite recently, using an external
magnetic field in a vapor cell of Rubidium atoms Bhushan  _et al._ (2019). A
practical question is if such results can be achieved at higher temperatures,
for example, for a room temperature multiband quantum memory.

In recent years, remarkable developments have been achieved to strongly couple
spin ensembles to cavity photons, leading to the emerging field of cavity
spintronics. Quanta of spin waves, magnons, are highly robust against
temperature Zhang  _et al._ (2015); Tabuchi  _et al._ (2014); Goryachev  _et
al._ (2014); Huebl  _et al._ (2013); Zhang  _et al._ (2014), and hence
significant magnon-photon hybridization and magnetically induced transparency
(MIT) have been successfully demonstrated even at room temperature Zhang  _et
al._ (2014). Tunable slow light and its conversion to fast light based upon
room temperature MIT has been theoretically shown recently Liu  _et al._
(2019). Besides, at strong magnon-photon interaction, a wide tunability of
slow light via applied magnetic field has been shown in Kong  _et al._ (2019).
These results demonstrate the promising value of these systems for practical
quantum memories Liu  _et al._ (2019). Here we explore how to split such a MIT
window into multiple bands for a room temperature multimode quantum memory.
Our idea is to exploit the coupling of magnons to thermal vibrations, which is
known to yield magnomechanically induced transparency (MMIT) Zhang  _et al._
(2016), in combination with multiple spin ensembles to achieve multiple bands
in MIT. We also discuss the emergence of Fano resonance in the output spectrum
and explore the suitable system parameters for its observation. Fano resonance
was first reported in the atomic systems Fano (1961), and it emerges due to
the quantum interference of different transition amplitudes which give minima
in the absorption profile. In later years, it has been discussed in different
physical systems, such as photonic crystal Rybin  _et al._ (2009), coupled
microresonators Xiao  _et al._ (2010), optomechanical system Kaur  _et al._
(2016). Recently, Fano-like asymmetric shapes have been experimentally
reported in a hybrid cavity magnomechanical system Zhang  _et al._ (2016).

Our model consists of two ferromagnetic insulators, specifically yttrium iron
garnets (YIGs), hosting long-lived magnons at room temperature, placed inside
a three-dimensional (3D) microwave cavity; we remark that another equivalent
embodiment of our model could be to place the YIGs on top of a superconducting
co-planar waveguide, which can have further practical significance being an
on-chip device Hisatomi  _et al._ (2016). Specific benefits of YIG as the host
of spin ensemble over other systems, such as paramagnetic spin ensembles in
nitrogen-vacancy centers is due to its high spin density of $2.1\times
10^{22}$ $\mu$B $cm^{-3}$ ($\mu$B is the Bohr magneton) and high room
temperature spin polarization below the Curie temperature (559 K). In addition
to multimode quantum memories, our results can be directly advantageous for
readily integrated microwave circuit applications at room temperature such as
multimode quantum transducers coupling different systems at different
frequencies Kim  _et al._ (2016), tunable frequency quantum sensors Huang  _et
al._ (2014) or fast light enhanced gyroscopes Kittel (1948). In addition to
the magnetic dipole interaction between the cavity field and the spin
ensemble, we take into account coupling between the magnons and the quanta of
YIG lattice vibrations, phonons, arising due to the magnetostrictive force
Zhang  _et al._ (2016). We only consider the Kittel mode Li  _et al._ (2018)
of the ferromagnetic resonance modes of the magnons. Such three-body quantum
systems can be of fundamental significance to examine macroscopic quantum
phenomena towards thermodynamic limit and quantum to classical transitions Li
and Zhu (2019).

In our model, tunable slow and fast light emerges as a natural consequence of
tunable splitting of MIT window. Slow-light propagation at room temperature
has been investigated recently in a cavity-magnon system and the group delays
are found to be in the $\sim\mu$s range Liu  _et al._ (2019). In a single YIG
magnomechanical system with strong magnon-photon coupling strength, slow-light
has achieved with a maximum group delay of $<0.8$ ms Kong  _et al._ (2019). In
this paper, we discuss the slow and fast light in a two YIGs magnomechanical
system. Further, we exlain the group delay depends on the tunability of the
magnon-photon coupling of the first YIG (YIG1) with respect to magnon-photon
coupling of the second YIG (YIG2). This not only helps to achieve larger group
delays at weak magnon-photon coupling, but also increase the group delay of
the transmitted probe field by one order of magnitude, which is not possible
with a single YIG system Kong  _et al._ (2019).

The rest of the paper is organized as follows: We describe the model system in
Sec. II and present dynamical equations with steady-state solutions. The
results and discussions for MMIT are presented in the Sec. III. We discuss the
emergence and tunability of the multiple Fano resonances in Sec. IV. Next, in
Sec. V, we present the transmission of the probe field and discuss the group
delays for slow and fast light propagation. Finally, in Sec. VI, we present
the conclusion of our work.

Figure 1: (color online) A schematic illustration of a hybrid cavity
magnomechanical system. It consists of two ferromagnetic yttrium iron garnet
(YIG) spheres placed inside a microwave cavity. A Bias magnetic field is
applied in the $z$ direction on each sphere, which excites the magnon modes,
and these modes are strongly coupled with the cavity field. The bias magnetic
fields activate the magnetostrictive (magnon-phonon) interaction in both YIGs.
The single-magnon magnomechanical coupling strength is very weak Zhang  _et
al._ (2016), and it depends on the spheres diameters and external bias field
directions. Either by considering a larger YIG1 sphere or adjusting the
direction of the bias field on it, the magnomechanical coupling of this sphere
can be ignored. Here, we assume the direction of the bias field on YIG1 such
that the single-magnon magnomechanical interaction becomes very weak and can
be ignored Zhang  _et al._ (2016). However, the magnomechanical interaction of
YIG2 is enhanced by directly driving its magnon mode via a microwave drive (y
direction). This microwave drive plays the role of a control field in our
model. Cavity, phonon, magnon modes are labeled as $a$, $b$, and $m_{i}$
($i=1,2$), respectively.

##  II System Hamiltonian and Theory

We consider a hybrid cavity magnomechanical system that consists of two YIG
spheres placed inside a microwave cavity, as shown in Fig. 1. A uniform bias
magnetic field (z-direction) is applied on each sphere, which excites the
magnon modes and these modes are coupled with the cavity field via magnetic
dipole interaction. The excitation of the magnon modes inside the spheres
leads to the variation magnetization that results in the deformation of their
lattice structures. The magnetostrictive force causes vibrations of the YIGs
which establishes magnon-phonon interaction in these spheres. The single-
magnon magnomechanical coupling strength is very weak Zhang  _et al._ (2016),
and it depends on the spheres diameters and external bias field directions.
Either by considering a larger YIG1 sphere or adjusting the direction of the
bias magnetic field on it, the magnomechanical coupling of this sphere can be
ignored Li and Zhu (2019). Here, we assume the direction of the bias field on
YIG1 such that the single-magnon magnomechanical interaction becomes very weak
and can be ignored Zhang  _et al._ (2016). However, the magnomechanical
interaction of YIG2 is enhanced by directly driving its magnon mode via a
external microwave drive. This microwave drive plays the role of a control
field in our model. In addition, the cavity is driven by a weak probe field.
In this work, we consider high quality YIG spheres, each has a 250 $\mu$m
diameter, and composed of ferric ions Fe+3 of density $\rho=4.22\times
10^{27}m^{-3}$. This causes a total spin $S=5/2\rho V_{m}=7.07\times 10^{14}$,
where $V_{m}$ is the volume of the YIG and $S$ is the collective spin operator
which satisfy the algebra;
$[S_{\alpha},S_{\beta}]=i\varepsilon^{\alpha\beta\gamma}S_{\gamma}$. The
Hamiltonian of the system reads Li and Zhu (2019)

$\begin{split}H/\hbar&=\omega_{a}\hat{a}^{\dagger}\hat{a}+\omega_{b}\hat{b}^{\dagger}\hat{b}+\sum_{j=1}^{2}[\omega_{j}{\hat{m}}_{j}^{{\dagger}}{\hat{m}}_{j}+g_{j}({\hat{m}_{j}}^{{\dagger}}\hat{a}+m_{j}\hat{a}^{{\dagger}})]\\\
&+g_{mb}\hat{m}^{\dagger}_{2}\hat{m}_{2}(\hat{b}+\hat{b}^{\dagger})+i(\Omega_{d}\hat{m}^{\dagger}_{2}e^{-i\omega_{d}t}-\Omega_{d}^{\star}\hat{m}_{2}e^{i\omega_{d}t})\\\
&+i(\hat{a}^{\dagger}\varepsilon_{p}e^{-i\omega_{p}t}-\hat{a}\varepsilon_{p}^{\star}e^{i\omega_{p}t})\end{split}$
(1)

where $a^{\dagger}(a)$ and $b^{\dagger}(b)$ are the creation (annihilation)
operators of the cavity and phonon modes, respectively. The resonance
frequencies of the cavity, phonon and magnon modes are denoted by
$\omega_{a}$, $\omega_{b}$ and $\omega_{j}$, respectively. Moreover, $m_{j}$
is the bosonic operator of the Kittle mode of frequency $\omega_{j}$ and its
coupling strength with the cavity mode is given by $g_{j}$. The frequency
$\omega_{j}$ of the magnon mode $m_{j}$ can be determined by using
gyromagnetic ratio $\gamma_{j}$ and external bias magnetic field $H_{j}$ i.e.,
$\omega_{j}=\gamma_{j}H_{j}$ with $\gamma_{j}/2\pi=28$ GHz. The Rabi frequency
$\Omega_{d}=\sqrt{5}/4\gamma\sqrt{N}B_{0}$ Li  _et al._ (2018), represents the
coupling strength of the drive field of amplitude $B_{0}$ and frequency
$\omega_{d}$. Furthermore, in Eq. (1), $\omega_{p}$ is the probe field
frequency having amplitude $\varepsilon_{p}$ which can be expressed as;
$\varepsilon_{p}=\sqrt{2P_{p}\kappa_{a}/\hbar\omega_{p}}$.

Note that in Eq. (1), we have ignored the non-linear term
$K\hat{m}^{\dagger}_{j}\hat{m}^{\dagger}_{j}\hat{m}_{j}\hat{m}_{j}$ that may
arise due to strongly driven magnon mode Wang  _et al._ (2018, 2016). To
ignore this nonlinear term, we must have $K\left\lvert\langle
m_{2}\rangle\right\rvert^{3}\ll\Omega$, and for the system parameters we
consider in this work, this condition always satisfies. The Hamiltonian in Eq.
(1) is written after applying the rotating-wave approximation in which fast
oscillating terms
$g_{j}(\hat{a}\hat{m}_{j}+\hat{a}^{\dagger}\hat{m}^{\dagger})$ are dropped.
This is valid for $\omega_{a},\omega_{j}\gg g_{j},\kappa_{a},\kappa_{mj}$
which is the case to be considered in the present work. Where $\kappa_{a}$ and
$\kappa_{mj}$ are the decay rates of the cavity and magnon modes,
respectively. In the frame rotating at the drive frequency $\omega_{d}$, the
Hamiltonian of the system is given by

$\begin{split}H/\hbar=&\Delta_{a}\hat{a}^{\dagger}\hat{a}+\omega_{b}\hat{b}^{\dagger}\hat{b}+\sum_{j=1}^{2}[\Delta_{mj}{\hat{m}}_{j}^{{\dagger}}{\hat{m}}_{j}+g_{j}({\hat{m}}_{j}^{{\dagger}}\hat{a}+\\\
&m_{j}\hat{a}^{{\dagger}})]+g_{mb}\hat{m}^{\dagger}_{2}\hat{m}_{2}(\hat{b}+\hat{b}^{\dagger})+i(\Omega_{d}\hat{m}^{\dagger}_{2}-\Omega_{d}^{\star}\hat{m}_{2})+\\\
&i(\hat{a}^{\dagger}\varepsilon_{p}e^{-i\delta
t}-\hat{a}\varepsilon^{\star}_{p}e^{i\delta t}),\end{split}$ (2)

here, $\Delta_{a}=\omega_{a}-\omega_{d}$, $\Delta_{mj}=\omega_{j}-\omega_{d}$,
and $\delta=\omega_{p}-\omega_{d}$. The quantum Heisenberg-Langevin equations
based on the Hamiltonian in Eq. (2) can be written as

$\begin{split}&\dot{\hat{a}}=-i\Delta_{a}\hat{a}-i\sum_{j=1}^{2}g_{j}\hat{m}_{j}-\kappa_{a}\hat{a}+\varepsilon_{p}e^{-i\delta
t}+\sqrt{2\kappa_{a}}\hat{a}^{in}(t),\\\
&\dot{\hat{b}}=-i\omega_{b}\hat{b}-ig_{mb}{\hat{m}}^{\dagger}_{2}\hat{m}_{2}-\kappa_{b}\hat{b}+\sqrt{2\kappa_{b}}\hat{b}^{in}(t),\\\
&\dot{\hat{m}}_{1}=-i\Delta_{m1}\hat{m}_{1}-ig_{1}\hat{a}-\kappa_{m1}m_{1}+\sqrt{2\kappa_{m1}}\hat{m}_{1}^{in}(t),\\\
&\dot{\hat{m}}_{2}=-i\Delta_{m2}\hat{m}_{2}-ig_{2}\hat{a}-\kappa_{m2}m_{2}-ig_{mb}\hat{m}_{2}(\hat{b}+\hat{b}^{\dagger})\\\
&+\Omega_{d}+\sqrt{2\kappa_{m2}}\hat{m}_{2}^{in}(t).\end{split}$ (3)

Where $\kappa_{b}$ is the dissipation rate of the phonon mode, and
$\hat{b}^{in}(t)$, $\hat{m}_{j}^{in}(t)$ and $\hat{a}^{in}(t)$ are the vacuum
input noise operators which have zero mean values i.e.,
$\langle\hat{q}_{\text{in}}\rangle=0$ Xiong and Wu (2018); Aspelmeyer  _et
al._ (2014), and ($q=a,m,b$). The magnon mode $m_{2}$ is strongly driven by a
microwave drive that causes a large steady-state amplitude $\left\lvert\langle
m_{2s}\rangle\right\rvert\gg 1$ of magnon mode, and due to beam splitter
interaction, this leads to the large steady-state amplitude of the cavity mode
$\left\lvert\langle a_{s}\rangle\right\rvert\gg 1$. Consequently, we can
linearize the quantum Langevin equations around the steady-state values and
take only the first-order terms in the fluctuating operator:
$\langle\hat{O}\rangle={O}_{s}+\hat{O}_{-}e^{-i\delta t}+\hat{O}_{+}e^{i\delta
t}$ Huang and Agarwal (2011), here $\hat{O}=a,b,m_{j}$. First, we consider the
zero-order solution, namely, steady-state solutions which are given by

$\begin{split}&a_{s}=-i\sum_{1,2}\frac{g_{j}m_{js}}{\kappa_{a}+i\Delta_{a}},\\\
&b_{s}=\frac{-ig_{mb}\mid{m_{2s}}\mid^{2}}{\kappa_{b}+i\omega_{b}},\\\
&m_{1s}=\frac{-ig_{1}a_{s}}{\kappa_{m1}+i\Delta_{m1}},m_{2s}=\frac{\Omega_{d}-ig_{2}a_{s}}{\kappa_{m2}+i\tilde{\Delta}_{m2}},\\\
&\tilde{\Delta}_{m2}=\Delta_{m2}+g_{mb}(b_{s}+b_{s}^{\star}).\end{split}$ (4)

We assume that the coupling of the external microwave drive on magnon mode
$m_{2}$ is much stronger than the amplitude $\epsilon_{p}$ of the probe field.
Under this assumption, the linearized quantum Langevin equations can be solved
by considering the first-order perturbed solutions and ignoring all higher
order terms of $\epsilon_{p}$. The solution for the cavity mode is given by

$a_{-}=\varepsilon_{p}\left[A^{\prime}+C_{1}^{\prime}+\frac{g_{2}^{2}}{\beta^{\prime}}+\frac{\alpha^{\star}\alpha^{\prime}}{\beta^{\star}\beta^{\prime}+A^{\star}-C_{1}^{\star}+\frac{g_{2}^{2}}{\beta^{\star}}}\right]^{-1},$
(5)

where

$\begin{split}&A=\kappa_{a}+i(\Delta_{a}+\delta),B=\frac{G_{mb}^{2}\omega_{b}}{\omega_{b}^{2}-\delta^{2}+i\delta\kappa_{b}},\\\
&C_{1}=\frac{g_{1}^{2}}{\kappa_{m1}+i(\Delta_{m1}+\delta)},C_{2}=\frac{g_{2}^{2}}{\kappa_{m2}+i(\tilde{\Delta}_{m2}+\delta)},\\\
&A^{\prime}=\kappa_{a}+i(\Delta_{a}-\delta),B^{\prime}=\frac{G_{mb}^{2}\omega_{b}}{\omega_{b}^{2}-\delta^{2}-i\delta\kappa_{b}},\\\
&C_{1}^{\prime}=\frac{g_{1}^{2}}{\kappa_{m1}+i(\Delta_{m1}-\delta)},C_{2}^{\prime}=\frac{g_{2}^{2}}{\kappa_{m2}+i(\tilde{\Delta}_{m2}-\delta)},\\\
&\alpha=\frac{g_{2}^{2}B}{C_{2}+iB},\alpha^{\prime}=\frac{g_{2}^{2}B^{\prime}}{C_{2}^{\prime}+iB^{\prime}},\\\
&\beta=C_{2}-i\frac{C_{2}^{\prime^{\star}}B}{C_{2}^{\prime^{\star}}+i{B}},\beta^{\prime}=C_{2}^{\prime}-i\frac{C_{2}^{\star}B^{\prime}}{C_{2}^{\star}+i{B^{\prime}}}.\end{split}$

Here $G_{mb}=i\sqrt{2}g_{mb}m_{2s}$ is the effective magnon-phonon coupling.
We use the input-output relation for the cavity field
$\varepsilon_{out}={\varepsilon}_{in}-2\kappa_{a}\langle a\rangle$ Gardiner
_et al._ (2004), and the amplitude of the output field can be written as

Figure 2: (Color online). Absorption Re$[\varepsilon_{out}]$ profiles are
shown against the normalized probe field detuning $\delta/\omega_{b}$. (a)
$g_{1}=g_{mb}=0$, $g_{2}/2\pi=1.2$ MHz and (b) $g_{1}=0$, $g_{2}/2\pi=1.2$
MHz, $G_{mb}/2\pi=2.0$ MHz (c) $g_{1}/2\pi=g_{2}/2\pi=1.2$ MHz,
$G_{mb}/2\pi=2$ MHz and (d) $g_{1}/2\pi=g_{2}/2\pi=1.2$ MHz, $G_{mb}/2\pi=3.5$
MHz. The other parameters are given in Sec. III.

$\varepsilon_{out}=\frac{2\kappa_{a}a_{-}}{\varepsilon_{p}}.$ (6)

The real and imaginary parts of $\varepsilon_{out}$ account for in-phase
(absorption) and out of phase (dispersion) output field quadratures at probe
frequency.

##  III MMIT windows profile

For the numerical calculation, we use parameters from a recent experiment on a
hybrid magnomechanical system Zhang  _et al._ (2016), unless stated
differently. Frequency of the cavity field $\omega_{a}/2\pi=10$ GHz,
$\omega_{b}/2\pi=10$ MHz, $\kappa_{b}/2\pi=100$ Hz, $\omega_{1,2}/2\pi=10$
GHz, $\kappa_{a}/2\pi=2.1$ MHz, $\kappa_{m1}/2\pi=\kappa_{m2}/2\pi=0.1$ MHz,
$g_{1}/2\pi=g_{2}/2\pi=1.5$ MHz, $G_{mb}/2\pi=3.5$ MHz,
$\Delta_{a}=\omega_{b}$, $\Delta_{mj}=\omega_{b}$, $\omega_{d}/2\pi=10$ GHz.

We first illustrate the physics behind the multiband transparency by
systematically investigating the role of different couplings in the model.
Fig. 2 displays the response of the probe field in the absorption spectrum of
the output field for different coupling strengths. In Fig. 2(a), we assume the
magnon-phonon coupling ($g_{mb}$) and magnon mode $m_{1}$ coupling ($g_{1}$)
with the cavity are absent. Therefore, only magnon mode $m_{2}$ is coupled
with the cavity. Under these considerations, we observe a magnon induced
transparency (MIT) in which a typical Lorentzian peak of the output spectrum
of the simple cavity splits into two peaks with a single dip, as shown in Fig.
2(a). The width of this transparency window can be controlled via microwave
driving field power and the magnon-photon coupling $g_{2}$. On increasing the
coupling strength $g_{2}$ the width of the window increases, and vice versa.

Figure 3: (Color online) Dispersion Im$[\varepsilon_{out}]$ profiles are shown
against the normalized probe detuning $\delta/\omega_{b}$. (a)
$g_{1}=g_{mb}=0$ and $g_{2}/2\pi=1.2$ MHz and (b) $g_{1}=0$, $g_{2}/2\pi=1.2$
MHz, $G_{mb}/2\pi=2.0$ MHz (c), (d) $g_{1}/2\pi=g_{2}/2\pi=1.2$ MHz, and (c)
$G_{mb}/2\pi=2$ MHz and (d) $G_{mb}/2\pi=3.5$ MHz. The other parameters are
given in Sec. III.

We observe two transparency windows in the absorption as we switch on the
magnon-phonon coupling ($g_{mb}$) and keeping $g_{1}=0$. Due to the non-zero
magnetostrictive interaction, single MIT window in Fig. 2(a) splits into
double window shown in Fig. 2(c). The right transparency window in Fig. 2(c)
is associated with magnon-phonon interaction, and this is so called
magnomechanically induced transparency (MMIT) Zhang  _et al._ (2016) window.
We can observe double MIT by removing magnon-phonon coupling $g_{mb}$, and
considering non-zero couplings between the magnon modes and the cavity field.

Finally, if we consider all three couplings simultaneously non-zero, then the
transparency window splits into three windows consist of four peaks and three
dips, this is shown in Fig. 2(c). In this case, one window is associated with
the magnomechanical interaction, and the rest of the two are induced by
magnon-photon couplings. The width and peaks separation of these windows
increases and broadens, respectively, at higher values of magnon-phonon
coupling $G_{mb}$, which can be seen in Fig. 2(d). Moreover, we have a
symmetric multi-window transparency profile where the splitting of the peaks
occurs at side-mode frequencies $\omega_{p}=\omega_{b}\pm\omega_{d}$.

In Figs. 3(a-d), we plot the dispersion spectrum of the output field versus
normalized frequency of the probe field. The single MIT dispersion spectrum in
the absence of YIG1 and magnon-phonon coupling $g_{mb}$ is shown in Fig. 3(a).
The dispersion spectra for the case of $g_{1}=0$, $g_{2}\neq 0$ and
$g_{mb}\neq 0$ is plotted in the Fig. 3(b). In the presence of all three
couplings, the dispersion spectrum of the output field is given in the Figs.
3(c-d). It is clear from Figs. 3(c-d), by the increase in the effective
magnon-phonon coupling $G_{mb}$, the transparency windows become wider. We
like to point out that the magnomechanically induced amplification (MMIA) of
the output field, in our system, can be obtained in the blue detuned regime;
$\Delta_{m2}=-\omega_{b}$.

Figure 4: (Color online) Fano line shapes in the asymmetric absorption
Re$[\varepsilon_{out}]$ profiles are shown against the normalized probe
frequency $\delta/\omega_{b}$. (a) $\Delta_{m2}=0.7\omega_{b}$, $g_{2}=1.5$
MHz, $g_{1}=g_{mb}=0$, and (b) $\Delta_{m2}=0.7\omega_{b}$, $g_{1}=0$,
$g_{2}=1.5$ MHz, $G_{mb}=3.5$ MHz. (c) $\Delta_{m1,2}=0.7\omega_{b}$,
$g_{1}=g_{2}/2\pi=1.5$ MHz and $G_{mb}/2\pi=3.5$MHz, and (d)
$\Delta_{m1,2}=\omega_{b}$ , $g_{1}=g_{2}/2\pi=1.5$ MHz and $G_{mb}/2\pi=3.5$
MHz. In all panels, $g_{1}=g_{2}/2\pi=1.5$ MHz, $G_{mb}/2\pi=3.5$ MHz, and
rest of the parameters are give in Sec. III.

##  IV Fano resonances in the output field

In the following, we discuss the emergence and physical mechanism of the Fano
line shapes in the output spectrum. The shape of the Fano resonance is
distinctly different than the symmetric resonance curves in the EIT, MIT,
optomechanically induced transparency (OMIT) and MMIT windows Qu and Agarwal
(2013); Zhang  _et al._ (2016). Fano resonance has observed in the systems in
which EIT has reported by a suitable selection of the system parameters Zhang
_et al._ (2016); Qu and Agarwal (2013); Ullah  _et al._ (2018); Ullah (2019a);
Yasir and Liu (2016); Akram  _et al._ (2015); Zhang  _et al._ (2017). The
physical origin of Fano resonance in the systems having optomechanical-like
interactions has explained due to the presence of non-resonant interactions.
For example, in a standard optomechanical system, if the anti-Stokes process
is not resonant with the cavity frequency, asymmetric Fano shapes appear in
the spectrum Qu and Agarwal (2013); Ullah  _et al._ (2018); Ullah (2019a). In
our system, this corresponds to $\Delta_{m1}\neq\omega_{b}$, because instead
of a cavity mode, magnon mode $m_{1}$ is coupled with phonon mode via
optomechanical-like interaction.

Figure 5: (Color online). The transmission $|t_{p}|^{2}$ spectrum as a
function of normalized probe field frequency $\delta/\omega_{b}$ is shown for
different values of $g_{1}$. (a) $g_{1}/2\pi=0.5$ MHz (b) $g_{1}/2\pi=0.8$ MHz
(c) $g_{1}/2\pi=1.2$ MHz (d) $g_{1}/2\pi=1.5$ MHz. In all panels,
$g_{2}/2\pi=1.5$ MHz, $G_{mb}/2\pi=3.5$ MHz and the other parameters are given
in Sec. III.

The asymmetric Fano shapes can be seen in Figs. 4(a-c) for different non-
resonant cases, where the absorption spectrum of the output field as a
function of normalized detuning $\delta/\omega_{b}$ is shown. In Fig. 4(a), we
consider $g_{1}=g_{mb}=0$, and coupling of the magnon mode $m_{2}$ with the
cavity is non-zero. Due to the presence of non-resonant process
($\Delta_{m2}=0.7\omega_{m}$), the absorption spectrum of the symmetric MIT
(Fig. 2(a)) profile changes into asymmetric window profile, as shown in Fig.
4(a). Such asymmetric MIT band can be related to Fano-like resonance, emerging
frequently in optomechanical systems Qu and Agarwal (2013); Ullah  _et al._
(2018); Ullah (2019a); Yasir and Liu (2016); Akram  _et al._ (2015). If we
remove YIG1 and consider only YIG2 is coupled with the cavity mode, and
$\Delta_{m2}=0.7\omega_{m}$. We observe double Fano resonance in the output
spectrum, which is shown in Fig. 4(b). Similarly, in the presence of all three
couplings and $\Delta_{m1,m2}=0.7\omega_{m}$, the double Fano resonance goes
over to a triple Fano profile, as shown in Fig. 4(c). This is because the
cavity field can be build up by three coherent routes provided by the three
coupled systems (the magnons, cavity, and phonon modes), and that can
interfere with each other. The Fano resonances disappear when we consider a
resonant case $\Delta_{m1}=\Delta_{m2}=\omega_{b}$, as shown in Fig. 4(d).

##  V Numerical results for slow and fast lights

Here we investigate the transmission and group delay of the output signal, and
show the effect of the magnon-photon and magnon-phonon couplings on the
transmission spectrum. From Eq. (6), the rescaled transmission field
corresponding to the probe field can be expressed as

$t_{p}=\frac{\varepsilon_{p}-2\kappa_{a}a_{-}}{\varepsilon_{p}}.$ (7)

In Figs. 5(a-d), we plot the transmission spectrum of the probe field against
the scaled detuning $\delta/\omega_{b}$, for different values of $g_{1}$. It
is clear from Fig. 5(a), the transmission peak associated with the magnon-
photon coupling of YIG1 is smaller than the other two peaks. This is because
in Fig. 5(a) $g_{1}$ coupling is weaker than the other two interactions
$g_{2}$ and $G_{mb}$ present in the system. By increasing the coupling
strength $g_{1}$, the peak of the middle transparency profile grows up in
height and reaches close to unity, as shown in Figs. 5(b-c). In addition, Fig.
5(d) shows that the width of the transparency window can be increased at
higher higher values of the magnon-photon coupling $g_{1}$.

Figure 6: (Color online). The transmission $|t_{p}|^{2}$ spectrum as a
function of normalized probe field frequency $\delta/\omega_{b}$ is shown. (a)
$G_{mb}/2\pi=0.5$ MHz (b) $G_{mb}/2\pi=1.0$ MHz. In (c) $g_{2}/2\pi=0.4$ MHz,
and (d) $g_{2}/2\pi=0.8$ MHz. The other parameters are same as in Fig. 5.

In Figs. 6(a-b), the transmission spectrum of the probe field as a function of
dimensionless detuning is shown for different values of $G_{mb}$. In Figs.
6(a-b), we consider both $g_{1}$ and $g_{2}$ to be the same in the strong
coupling regime. However, the effective coupling
$\tilde{g}_{2}=g_{2}\alpha_{s}$ depends on the steady-state amplitude of the
cavity field $\alpha_{s}$ which depends on the $m_{2s}$. Consequently,
$\tilde{g}_{2}$ and $G_{mb}$ are related and it can be seen from Eq. (4). For
a smaller value of $G_{mb}$ in Fig. 6(a), we have two small peaks associated
with $g_{2}$ and $G_{mb}$, in addition, the third-highest peak is associated
with $g_{1}$. For a fixed value of $g_{mb}$, if we increase $G_{mb}$, it
increases $\tilde{g}_{1}$, and the peaks associated with these two couplings
become more visible, as shown in Fig. 6(b). Similarly, in Fig. 6(c-d), we
observe a similar increase in the height of two peaks associated with $g_{2}$
and $G_{mb}$, for the variation in $g_{2}$.

The phase $\phi_{t}$ of the transmitted probe field $t_{p}$ is given by the
relation $\phi_{t}=\text{Arg}[t_{p}]$. The plot of $\phi_{t}$ as a function of
normalized detuning $\delta/\omega_{b}$ is shown in Fig. 7. In the inset of
Fig. 7(a), we consider both $g_{1}$ and $g_{mb}$ are switched off, and only
$g_{2}$ is non-zero. This gives a conventional phase of the transmitted field
with a single MIT curve, which appears similar to the standard single OMIT
curve Qu and Agarwal (2013). In Fig. 7(b), we switch-off the YIG1 coupling
with the field ($g_{1}=0$), and the other two couplings are present
($g_{mb}\neq 0,g_{2}\neq 0$), due to which the single transparency window
splits into a double window. If we keep all three couplings non-zero, we get
triple transparency window which is shown in Fig. 7(c).

Figure 7: (Color online) The phase $\phi_{t}$ of the transmitted probe field
versus normalized detuning $\delta/\omega_{b}$ for different coupling
strengths. (a) $g_{1}=g_{mb}=0$, (b) $g_{1}=0$, $g_{2}/2\pi=1.5$ MHz,
$G_{mb}/2\pi=4$ MHz (c) $g_{1}/2\pi=g_{2}/2\pi=1.5$ MHz, and $G_{mb}/2\pi=4$
MHz. Rest of parameters are given in Sec. III.

Figure 8: (Color online) Group delay $\tau_{g}$ of the output probe field
against the amplitude of the magnetic field $B_{0}$ for (a) $g_{1}=0$, and (b)
$g_{1}/2\pi=1.5$ MHz. The other parameter are $g_{2}/2\pi=1.5$,
$G_{mb}/2\pi=3.5$ MHz, $\kappa_{b}/2\pi=100$ Hz,
$\kappa_{m1}/2\pi=\kappa_{m2}/2\pi=0.1$ MHz, $\kappa_{a}/2\pi=2.1$ MHz and
$\Omega_{d}=1.2$ THz.

The transmitted probe field phase is associated with the group delay
$\tau_{g}$ of the output field and it is defined as

$\tau_{g}=\frac{\partial\phi(\omega_{p})}{\partial\omega_{p}},$ (8)

which means a more rapid phase dispersion leads to a larger group delays and
vice versa. In addition, a negative slope of the phase represents a negative
group delay or fast light ($\tau_{g}<0$) whereas, a positive slope of the
transmitted field indicates positive group delay or slow light ($\tau_{g}>0$).
From Fig. 7, we observe that in the regime of the narrow transparency window,
there is a rapid variation in the probe phase, and this rapid phase dispersion
can lead to a significant group delay.

Fig. 8 shows the group delay $\tau_{g}$ can be tuned by the variation of the
bias magnetic field $B_{0}$ applied on YIG2. In the absence of YIG1 (Fig.
8(a)), we have a lower slope of Eq. (8), as a result, a maximum group delay of
$\tau_{g}=1$ ms is achieved. This group delay can be enhanced by one order of
magnitude once second YIG is introduced see Fig. 8(b). The slope of Fig. 8(b)
become steeper and the time delay for slow light is increased up to 13.8 ms.
This shows the two YIGs system is a good choice to observe a longer group
delay in a magnomechanical system while a single YIG system cannot do so.
Moreover, the numerical value of the group delay $\tau_{g}$ can be tuned from
positive (slow light) to negative (fast light) by tuning the magnon field
detuning $\Delta_{m1}=\omega_{b}$ to $\Delta_{m1}=-\omega_{b}$. Here, it is
worth mentioning that Fig. 8(b) can be switched into fast light with a maximum
group delay in the order of $\tau_{g}\approx-1.4$ ms in the presence of both
YIGs and we not show in the figure. This negative group delay for the fast
light propagation is one order of magnitude greater than a single YIG
magnomechanical system Kong  _et al._ (2019).

Figure 9: (Color online) The group delay $\tau_{g}$ of the transmitted probe
field as a function of the driving power $P_{d}$ for several values of the
magnon-photon couplings. (a) $\Delta_{m1}=\omega_{b}$, (b)
$\Delta_{m1}=-\omega_{b}$, and the other parameter are same as given in Fig.
8.

Finally, we investigate the control of group delay with the external microwave
driving power and magnon-photon couplings. For this purpose, in Figs. 9(a-b),
we plot $\tau_{g}$ against the driving power for different strengths of the
magnon-photon coupling of YIG1 with respect to the coupling frequency of YIG2.
Fig. 9(a) shows that the magnitude of the group delay increases with the
increase of $g_{1}$ corresponding to $g_{2}$, which indicates an enhanced
group delay of the transmitted probe field in a two-YIG system. We tuned the
coupling strength of YIG1 ($g_{1}$) for different values via keeping the
coupling strength of YIG2 ($g_{2}$) constant. This shows increasing the
magnon-photon coupling strength increases the group delay of the transmitted
probe field and vice versa. This helps us to obtain larger group delays at
relatively weak magnon-photon coupling strengths which is not otherwise
possible with a single YIG magnomechanical system Kong  _et al._ (2019).
Similar results can also be obtained by increasing the magnon-photon coupling
$g_{2}$ and fixing $g_{1}$. For the blue detuned regime
$\Delta_{m1}=-\omega_{b}$, group delay becomes negative. However, the effect
of magnon-photon couplings remains the same, as shown in Fig. 9(b). From Fig.
8 and Fig. 9, we see that two YIGs magnomechanical system provides not only
extra tunability, but also drastically enhances the group delays compared to
single YIG system studied in Ref. Kong  _et al._ (2019). Our system can be
used as a tunable switch, which can be controlled via different system
parameters, and our results are comparable with the existing proposals based
on various hybrid quantum systems Genes  _et al._ (2008); Ullah (2019b);
Tarhan  _et al._ (2013); Saif  _et al._ (2019).

##  VI Conclusion

We have investigated the transmission and absorption spectrum of a weak probe
field under a strong control field in a hybrid magnomechanical system in the
microwave regime. Due to the presence of a nonlinear phonon-magnon
interaction, we observed magnomechanically induced transparency (MMIT), and
the photon-magnon interactions lead to magnon induced transparency (MIT). We
found single MMIT, a result of the single-phonon process, and found two MIT
windows in the output probe spectra due to the presence of two magnon-photon
interactions. This is demonstrated by plotting the absorption, dispersion, and
transmission of the output field. We discussed the emergence of Fano
resonances in the output field spectrum of the probe field. These asymmetric
line shapes appeared due to the presence of anti-Stokes processes in the
system. We examined conditions of slow and fast light propagation in our
system, which can be controlled by different system parameters. It has shown
that in a two YIGs magnomechanical system, the tunability of the first
coupling strength (YIG1) corresponding to the coupling of the second YIG
(YIG2) has an immense effect on the slow and fast light and vice versa. This
not only helped to investigate larger group delays at a weak magnon-photon
coupling but also enhanced the group delay of the transmitted probe field,
which is not possible in a single YIG system. Our results suggest that this
system may find its applications to implement multi-band quantum memories Liu
_et al._ (2019).

## Acknowledgment

We acknowledge Prof. M. Cengiz Onbasli for fruitful discussions. We also thank
the anonymous referees for their valuable comments which greatly improved our
manuscript.

## References

  * Afzelius  _et al._ (2009) M. Afzelius, C. Simon, H. de Riedmatten, and N. Gisin, Phys. Rev. A 79, 052329 (2009).
  * Sinclair  _et al._ (2014) N. Sinclair, E. Saglamyurek, H. Mallahzadeh, J. A. Slater, M. George, R. Ricken, M. P. Hedges, D. Oblak, C. Simon, W. Sohler, and W. Tittel, Phys. Rev. Lett. 113, 053603 (2014).
  * Jobez  _et al._ (2016) P. Jobez, N. Timoney, C. Laplane, J. Etesse, A. Ferrier, P. Goldner, N. Gisin, and M. Afzelius, Phys. Rev. A 93, 032327 (2016).
  * McGloin  _et al._ (2000) D. McGloin, M. H. Dunn, and D. J. Fulton, Phys. Rev. A 62, 053802 (2000).
  * Ying  _et al._ (2014) K. Ying, Y. Niu, D. Chen, H. Cai, R. Qu, and S. Gong, Journal of Modern Optics 61, 631 (2014).
  * Bhushan  _et al._ (2019) S. Bhushan, V. S. Chauhan, D. M, and R. K. Easwaran, Physics Letters A 383, 125885 (2019).
  * Zhang  _et al._ (2015) D. Zhang, X.-M. Wang, T.-F. Li, X.-Q. Luo, W. Wu, F. Nori, and J. You, npj Quantum Information 1, 1 (2015).
  * Tabuchi  _et al._ (2014) Y. Tabuchi, S. Ishino, T. Ishikawa, R. Yamazaki, K. Usami, and Y. Nakamura, Phys. Rev. Lett. 113, 083603 (2014).
  * Goryachev  _et al._ (2014) M. Goryachev, W. G. Farr, D. L. Creedon, Y. Fan, M. Kostylev, and M. E. Tobar, Phys. Rev. Applied 2, 054002 (2014).
  * Huebl  _et al._ (2013) H. Huebl, C. W. Zollitsch, J. Lotze, F. Hocke, M. Greifenstein, A. Marx, R. Gross, and S. T. B. Goennenwein, Phys. Rev. Lett. 111, 127003 (2013).
  * Zhang  _et al._ (2014) X. Zhang, C.-L. Zou, L. Jiang, and H. X. Tang, Phys. Rev. Lett. 113, 156401 (2014).
  * Liu  _et al._ (2019) Z. Liu, H. Xiong, and Y. Wu, IEEE Access 7, 57047 (2019).
  * Kong  _et al._ (2019) C. Kong, B. Wang, Z.-X. Liu, H. Xiong, and Y. Wu, Opt. Express 27, 5544 (2019).
  * Zhang  _et al._ (2016) X. Zhang, C.-L. Zou, L. Jiang, and H. X. Tang, Science Advances 2 (2016), 10.1126/sciadv.1501286.
  * Fano (1961) U. Fano, Phys. Rev. 124, 1866 (1961).
  * Rybin  _et al._ (2009) M. V. Rybin, A. B. Khanikaev, M. Inoue, K. B. Samusev, M. J. Steel, G. Yushin, and M. F. Limonov, Phys. Rev. Lett. 103, 023901 (2009).
  * Xiao  _et al._ (2010) Y.-F. Xiao, M. Li, Y.-C. Liu, Y. Li, X. Sun, and Q. Gong, Phys. Rev. A 82, 065804 (2010).
  * Kaur  _et al._ (2016) S. Kaur, B. Yao, Y.-S. Gui, and C.-M. Hu, Journal of Physics D: Applied Physics 49, 475103 (2016).
  * Hisatomi  _et al._ (2016) R. Hisatomi, A. Osada, Y. Tabuchi, T. Ishikawa, A. Noguchi, R. Yamazaki, K. Usami, and Y. Nakamura, Phys. Rev. B 93, 174427 (2016).
  * Kim  _et al._ (2016) N.-H. Kim, J. Jung, J. Cho, D.-S. Han, Y. Yin, J.-S. Kim, H. J. M. Swagten, and C.-Y. You, Applied Physics Letters 108, 142406 (2016), https://doi.org/10.1063/1.4945685 .
  * Huang  _et al._ (2014) M. Huang, Y. Zhang, F. Li, L. Zhang, R. S. Ruoff, Z. Wen, and Q. Liu, Scientific reports 4, 3878 (2014).
  * Kittel (1948) C. Kittel, Phys. Rev. 73, 155 (1948).
  * Li  _et al._ (2018) J. Li, S.-Y. Zhu, and G. S. Agarwal, Phys. Rev. Lett. 121, 203601 (2018).
  * Li and Zhu (2019) J. Li and S.-Y. Zhu, New Journal of Physics 21, 085001 (2019).
  * Wang  _et al._ (2018) Y.-P. Wang, G.-Q. Zhang, D. Zhang, T.-F. Li, C.-M. Hu, and J. Q. You, Phys. Rev. Lett. 120, 057202 (2018).
  * Wang  _et al._ (2016) Y.-P. Wang, G.-Q. Zhang, D. Zhang, X.-Q. Luo, W. Xiong, S.-P. Wang, T.-F. Li, C.-M. Hu, and J. Q. You, Phys. Rev. B 94, 224410 (2016).
  * Xiong and Wu (2018) H. Xiong and Y. Wu, Applied Physics Reviews 5, 031305 (2018), https://doi.org/10.1063/1.5027122 .
  * Aspelmeyer  _et al._ (2014) M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, Reviews of Modern Physics 86, 1391 (2014).
  * Huang and Agarwal (2011) S. Huang and G. S. Agarwal, Phys. Rev. A 83, 043826 (2011).
  * Gardiner  _et al._ (2004) C. Gardiner, P. Zoller, and P. Zoller,  _Quantum noise: a handbook of Markovian and non-Markovian quantum stochastic methods with applications to quantum optics_ (Springer Science & Business Media, 2004).
  * Qu and Agarwal (2013) K. Qu and G. S. Agarwal, Phys. Rev. A 87, 063813 (2013).
  * Ullah  _et al._ (2018) K. Ullah, H. Jing, and F. Saif, Phys. Rev. A 97, 033812 (2018).
  * Ullah (2019a) K. Ullah, Eur. Phys. J. D 73, 267 (2019a).
  * Yasir and Liu (2016) K. A. Yasir and W.-M. Liu, Scientific reports 6, 22651 (2016).
  * Akram  _et al._ (2015) M. J. Akram, F. Ghafoor, and F. Saif, Journal of Physics B: Atomic, Molecular and Optical Physics 48, 065502 (2015).
  * Zhang  _et al._ (2017) S. Zhang, J. Li, R. Yu, W. Wang, and Y. Wu, Scientific reports 7, 1 (2017).
  * Genes  _et al._ (2008) C. Genes, D. Vitali, and P. Tombesi, Phys. Rev. A 77, 050307 (2008).
  * Ullah (2019b) K. Ullah, Chinese Physics B 28, 114209 (2019b).
  * Tarhan  _et al._ (2013) D. Tarhan, S. Huang, and O. E. Müstecaplıoğlu, Phys. Rev. A 87, 013824 (2013).
  * Saif  _et al._ (2019) F. Saif, K. Ullah, and S. Watanabe, AIP Conference Proceedings 2067, 020002 (2019).


# On the Normalizing Constant of
the Continuous Categorical Distribution

Elliott Gordon-Rodriguez
Columbia University
<EMAIL_ADDRESS>
&Gabriel Loaiza-Ganem11footnotemark: 1
Layer6 AI
<EMAIL_ADDRESS>
Andres Potapczynski
New York University
<EMAIL_ADDRESS>
&John P. Cunningham
Columbia University
<EMAIL_ADDRESS>
Equal contribution.

###### Abstract

Probability distributions supported on the simplex enjoy a wide range of
applications across statistics and machine learning. Recently, a novel family
of such distributions has been discovered: the _continuous categorical_. This
family enjoys remarkable mathematical simplicity; its density function
resembles that of the Dirichlet distribution, but with a normalizing constant
that can be written in closed form using elementary functions only. In spite
of this mathematical simplicity, our understanding of the normalizing constant
remains far from complete. In this work, we characterize the numerical
behavior of the normalizing constant and we present theoretical and
methodological advances that can, in turn, help to enable broader applications
of the continuous categorical distribution. Our code is available at
https://github.com/cunningham-lab/cb_and_cc/.

##  1 Introduction

The _continuous categorical_ (CC) distribution is defined by the following
density function (Gordon-Rodriguez et al., 2020a):

$\displaystyle\mathbb{x}\sim\mathcal{CC}(\boldsymbol{\lambda})\iff
p(\mathbb{x};\boldsymbol{\lambda})\propto\prod_{i=1}^{K}\lambda_{i}^{x_{i}}.$
(1)

Here, $\mathbb{x}$ denotes a simplex-valued random variable, and
$\boldsymbol{\lambda}$ denotes a simplex-valued parameter, in other
words:111Note that the $K$-simplex is also commonly defined as
$\Delta^{K}=\\{\mathbb{x}\in\mathbb{R}^{K}_{+}:\sum_{i=1}^{K}=1\\}$. The two
definitions are equivalent, however $\mathbb{S}^{K}$ is a subset of
$\mathbb{R}^{K-1}$ with positive Lebesgue measure, whereas $\Delta^{K}$ is a
subset of $\mathbb{R}^{K}$ with zero Lebesgue measure. For this reason, using
$\mathbb{S}^{K}$ will facilitate our later arguments involving integrals on
the simplex.

$\displaystyle\mathbb{x},\boldsymbol{\lambda}\in\mathbb{S}^{K}:=\left\\{\boldsymbol{x}\in\mathbb{R}_{+}^{K-1}:\sum_{i=1}^{K-1}x_{i}\leq
1\right\\},$ (2)

where we additionally define the $K$th coordinates as the remainder:

$\displaystyle x_{K}$ $\displaystyle=1-\sum_{i=1}^{K-1}x_{i}$ (3)
$\displaystyle\lambda_{K}$ $\displaystyle=1-\sum_{i=1}^{K-1}\lambda_{i}.$ (4)

It is natural to contrast the CC with the similar-looking Dirichlet
distribution:

$\displaystyle\mathbb{x}\sim\mathrm{Dirichlet}(\boldsymbol{\alpha})\iff
p(\mathbb{x};\boldsymbol{\alpha})\propto\prod_{i=1}^{K}x_{i}^{\alpha_{i}-1},$
(5)

where again $\mathbb{x}\in\mathbb{S}^{K}$, but
$\boldsymbol{\alpha}\in\mathbb{R}_{+}^{K}$ now denotes a positive
unconstrained parameter vector with one more dimension than $\mathbb{x}$.

While the densities in Eq. 1 and Eq. 5 look similar, they hide very different
normalizing constants. In the Dirichlet case, it is well known that
(Dirichlet, 1839):

$\displaystyle\int_{\mathbb{S}^{K}}\prod_{i=1}^{K}x_{i}^{\alpha_{i}-1}d\mu(\mathbb{x})=\frac{\prod_{i=1}^{K}\Gamma(\alpha_{i})}{\Gamma(\sum_{i=1}^{K}\alpha_{i})},$
(6)

where $\Gamma(\alpha)=\int_{0}^{\infty}t^{\alpha-1}e^{-t}dt$ denotes the gamma
function and $\mu(\cdot)$ is the Lebesgue measure. On the other hand, the
normalizing constant of the CC admits the following closed form (Gordon-
Rodriguez et al., 2020a):

$\displaystyle\int_{\mathbb{S}^{K}}\prod_{i=1}^{K}\lambda_{i}^{x_{i}}d\mu(\mathbb{x})=\sum_{k=1}^{K}\frac{\lambda_{k}}{\prod_{i\neq
k}\log{\frac{\lambda_{k}}{\lambda_{i}}}},$ (7)

which contains elementary operations only. In spite of its mathematical
simplicity, this normalizing constant can be numerically hard to compute,
particularly in high dimensions. Moreover, Eq. 7 breaks down under equality of
parameters, i.e., whenever $\lambda_{i}=\lambda_{k}$ for some $i\neq k$,
because the denominator evaluates to zero. These issues will be the primary
focus of our exposition, in particular:

  * •

In Section 2.1, we characterize the numerical behavior of our normalizing
constant. We demonstrate that vectorized computation can suffer from
catastrophic cancellation, the severity of which depends on the proximity
between parameter values.

  * •

In Section 2.2, we rederive the normalizing constant as an inverse Laplace
transform, which in turn can be be evaluated using numerical inversion
algorithms. We show that this alternative computation strategy exhibits good
numerical behavior in the regime where catastrophic cancellation is most
severe.

  * •

In Section 2.3, we propose an orthogonal computational approach based on a
recursive property of the normalizing constant.

  * •

In Section 3.2, we generalize Eq. 7 to arbitrary parameter values, i.e.,
including equality of parameters $\lambda_{i}=\lambda_{k}$ for any $i\neq k$.
The resulting formula depends on an expectation that can be computed using
automatic differentiation.

We conclude this section with some remarks. First, note that in the
1-dimensional case, the CC distribution reduces to the _continuous Bernoulli_
distribution (Loaiza-Ganem and Cunningham, 2019), which arose in the context
of generative models of images (Kingma and Welling, 2014; Bond-Taylor et al.,
2021) and provided the original inspiration for the CC family. More generally,
the CC is closely related to the categorical cross-entropy loss commonly used
in machine learning (Gordon-Rodriguez et al., 2020b).

We also note that the CC can be rewritten using the exponential family
canonical form:

$\displaystyle\mathbb{x}\sim\mathcal{CC}(\boldsymbol{\eta})\iff
p(\mathbb{x};\boldsymbol{\eta})=\frac{1}{C(\boldsymbol{\eta})}e^{\boldsymbol{\eta}^{\top}\mathbb{x}},$
(8)

where $\eta_{i}=\log(\lambda_{i}/\lambda_{K})$ is the natural parameter, which
conveniently becomes unconstrained real-valued. Note that, like with
$\mathbb{x}$ and $\boldsymbol{\lambda}$, we will drop the $K$th coordinate to
denote $\boldsymbol{\eta}=(\eta_{1},\dots,\eta_{K-1})$, since
$\eta_{K}=\log(1)=0$ is fixed.222In principle, we could let $\eta_{K}$ vary
together with $\eta_{1},\dots,\eta_{K-1}$; Eqs. 9 and 8 would still hold,
since the additional term $e^{\eta_{K}x_{K}}$ in the density would compensate
the change in $C(\boldsymbol{\eta})$. However, such a model would be
overparameterized as it would be invariant to a parallel shift across all the
$\eta_{i}$. For mathematical conciseness, we keep $\eta_{K}$ fixed at 0 and
work with $\boldsymbol{\eta}\in\mathbb{R}^{K-1}$. In this notation, the
normalizing constant becomes:

$\displaystyle
C(\boldsymbol{\eta})=\sum_{j=1}^{K}\frac{e^{\eta_{k}}}{\prod_{i\neq
k}\left(\eta_{k}-\eta_{i}\right)},$ (9)

which, again, is undefined whenever $\eta_{i}=\eta_{k}$ for some $i\neq k$.

##  2 Numerical computation of the normalizing constant

Eq. 9 can be vectorized efficiently as follows:

Algorithm 1 Vectorized computation of the normalizing constant

Input: A parameter vector $\boldsymbol{\eta}$
Output: The normalizing constant $C(\boldsymbol{\eta})$

1:Compute a $K\times K$ matrix of differences
$M=[\eta_{k}-\eta_{i}]_{i,k=1}^{K}$ and add to this the identity matrix $I$. A
tensor with an additional batch dimension can be used if necessary.

2:Take the product of the rows of $M+I$, using log space as necessary.

3:Multiply the resulting vector componentwise with the vector
$[e^{\eta_{k}}]_{k=1}^{K}$, and sum up the terms.

###  2.1 Catastrophic cancellation

Algorithm 1 is easy to code up and adds little computational overhead to most
models. However, the summation in Step 3 involves positive and negative
numbers, which can result in catastrophic cancellation (Goldberg, 1991). _We
stress that the log-sum-exp trick, while useful for preventing overflow,
cannot address catastrophic cancellation (see Section 2.4)._ For example,
consider the case $K=5$ with $\boldsymbol{\eta}=(1,2,3,4)$. In single-
precision floating-point, the summation in Eq. 9 evaluates to:

$\displaystyle\sum_{k=1}^{K}\frac{e^{\eta_{k}}}{\prod_{i\neq
k}\left(\eta_{k}-\eta_{i}\right)}$
$\displaystyle=-0.45304695+1.847264-3.3475895+2.2749228+0.04166667$
$\displaystyle=0.363217.$ (10)

Note that the output is an order of magnitude smaller than (at least one of)
the summands, and as a result we have lost one digit in precision. As we
increase the dimension of the CC, the cancellation becomes more severe. For
example, in the case $K=10$ with $\boldsymbol{\eta}=(1,2,3,4,5,6,7,8,9)$, the
same summation becomes:

$\displaystyle\sum_{k=1}^{K}\frac{e^{\eta_{k}}}{\prod_{i\neq
k}\left(\eta_{k}-\eta_{i}\right)}$ $\displaystyle=-6.7417699\times
10^{-5}-7.3304126\times 10^{-4}+4.6494300\times 10^{-3}$ $\displaystyle\ \ \ \
\ \ -1.8957691\times 10^{-2}+5.1532346\times 10^{-2}+9.3386299\times 10^{-2}$
$\displaystyle\ \ \ \ \ \ +1.0879297\times 10^{-1}-7.3932491\times
10^{-2}+2.2329926\times 10^{-2}$ $\displaystyle\ \ \ \ \ \ -2.7557318\times
10^{-6}$ $\displaystyle=3.5982\times 10^{-4}.$ (11)

We have now lost 3 digits, since the leading summand is 3 orders of magnitude
greater than the output. If we continue increasing $K$ in the same way, by
$K=20$ we will have lost 6 digits, and by $K=25$ we are no longer able to
compute $C(\boldsymbol{\eta})$ to even a single significant digit. If we were
to use double-precision floating-point instead, by $K=40$ we will have lost 13
digits, and by $K=50$ we can no longer compute $C(\boldsymbol{\eta})$ to a
single significant digit.

We can summarize the problem as follows: _when $C(\boldsymbol{\eta})$ is of a
much lower order of magnitude than the leading term in the summation of Eq. 9,
numerical computation fails due to catastrophic cancellation_. To complicate
things further, the relationship between the two orders of magnitude (that of
$C(\boldsymbol{\eta})$ and that of the leading summand) is nontrivial. This
relationship depends on the relative size of the exponential terms
$e^{\eta_{k}}$ and the products of the differences $\eta_{k}-\eta_{i}$, and is
not straightforward to analyze. However, if two elements of
$\boldsymbol{\eta}$ are particularly close to one another, meaning that
$|\eta_{j_{1}}-\eta_{j_{2}}|$ is close to 0 for some $j_{1}$ and $j_{2}$, it
follows that the corresponding summands:

$\displaystyle\frac{e^{\eta_{j_{1}}}}{\prod_{i\neq{j_{1}}}\left(\eta_{j_{1}}-\eta_{i}\right)},\
\ \mathrm{and}\ \
\frac{e^{\eta_{j_{2}}}}{\prod_{i\neq{j_{2}}}\left(\eta_{j_{2}}-\eta_{i}\right)},$
(12)

are (a) large in magnitude (due to the term $\eta_{j_{1}}-\eta_{j_{2}}$ in the
denominator), (b) of opposing sign (due to the difference in
$\eta_{j_{1}}-\eta_{j_{2}}$ versus $\eta_{j_{2}}-\eta_{j_{1}}$), and (c)
similar in absolute value (since the terms in the product are approximately
equal). Thus, it is likely that the terms in Eq. 12 are of leading order and
catastrophically cancel each other out. Note that, as the dimensionality $K$
increases, it becomes more likely that _some_ pair of the components are close
to one another, and therefore computation becomes harder.

Figure 1: Scaling behavior of $C(\boldsymbol{\eta})$ relative to its summands
(from Eq. 9), as the dimension $K$ varies. Each point represents a random draw
of $\eta_{i}\overset{iid}{\sim}N(0,1)$, for which we compute
$C(\boldsymbol{\eta})$. Note that catastrophic cancellation depends on the
difference in order of $C(\boldsymbol{\eta})$ and its summands; in the green
region this difference is at most 8 orders of magnitude, so that single-
precision floating-point is sufficient. In the yellow region, it is between 8
and 16 orders of magnitude, so that single-precision fails due to catastrophic
cancellation, but double-precision succeeds. In the red region, both fail.

Another helpful intuition can be obtained by reasoning from the integral:

$\displaystyle
C(\boldsymbol{\eta})=\int_{\mathbb{S}^{K}}e^{\boldsymbol{\eta}^{\top}\mathbb{x}}d\mu(\mathbb{x}).$
(13)

As $K$ increases, the Lebesgue measure of the simplex decays like $1/K!$.333As
can be seen, for example, by taking Eq. 6 with $\alpha_{i}=0$ for all $i$.
Therefore, assuming the components of $\boldsymbol{\eta}$ are $O(1)$, we have
that $e^{\boldsymbol{\eta}^{\top}\mathbb{x}}=O(1)$ also, and therefore
$C(\boldsymbol{\eta})=O(1/K!)$. Under this assumption, we also have that
$\eta_{k}-\eta_{i}=O(1)$, and therefore the summands in Eq. 9 cannot decay
factorially fast (they may, but need not, decay at most exponentially due to
the product of $K-1$ terms of constant order in the denominator). Thus, such a
regime implies catastrophic cancellation is inevitable for large enough $K$.
On the other hand, when all the components of $\boldsymbol{\eta}$ are far from
one another, we are spared of such cancellation and Algorithm 1 succeeds,
including in high dimensions. We demonstrate these behaviors empirically in
the following experiments (see Figures 1 and 2).

####  2.1.1 Experiments

To evaluate the effectiveness of Algorithm 1 for computing
$C(\boldsymbol{\eta})$, we first implemented an oracle capable of correctly
computing $C(\boldsymbol{\eta})$ to within 4 significant figures (at a
potentially large computational cost). Our oracle is an implementation of Eq.
9 with arbitrary-precision floating-point, using the mpmath library (Johansson
et al., 2013). In particular, for a given $\boldsymbol{\eta}$, we ensure the
level of precision is set appropriately by computing Eq. 9 repeatedly at
increasingly high precision, until the outputs converge (to 4 significant
figures). Equipped with this oracle, we then drew $\boldsymbol{\eta}$ vectors
from a normal prior for a variety of dimensions $K$, to analyze the behavior
of $C(\boldsymbol{\eta})$.

First, we took $\eta_{i}\overset{iid}{\sim}N(0,1)$ and compared the magnitude
of $C(\boldsymbol{\eta})$ to the magnitude of the largest summand in Eq. 9.
The results are plotted in Figure 1, where we observe that
$C(\boldsymbol{\eta})$ decays rapidly in $K$, whereas the same is not true of
its (largest) summands. Thus, under this prior, Algorithm 1 is unsuccessful
except in low-dimensional settings.

Next, we let the spread of $\boldsymbol{\eta}$ vary by drawing
$\eta_{i}\overset{iid}{\sim}N(0,\sigma^{2})$, where $\sigma$ ranges between
$10^{-2}$ and $10^{2}$. We then plotted, for each $\sigma$, the highest value
of $K$ such that the output of Algorithm 1 agreed with the Oracle to 3
significant figures. We repeated the procedure using single- and double-
precision floating-point (orange and pink lines in Figure 2, respectively), as
well as two Laplace inversion methods that will be discussed in Section 2.2.
As $\sigma$ increases, the parameter values tend to move away from one
another, making computation easier and allowing for much higher dimensions. On
the other hand, when $\sigma$ decreases, the parameter values come closer,
bringing us into the unstable regions of Eq. 9, and Algorithm 1 fails for all
but just a few dimensions.

Figure 2: Scaling behavior of the numerical accuracy for computing
$C(\boldsymbol{\eta})$, as a function of the spread of the parameter values.
For each $\sigma$ (x-axis), we draw
$\eta_{i}\overset{iid}{\sim}N(0,\sigma^{2})$ for $i=1,\dots,40$, and we
compute $C(\eta_{1},\dots,\eta_{K-1})$, for each $K=3,\dots,40$. We then show
on the y-axis the highest $K$ such that the numerical value of
$C(\eta_{1},\dots,\eta_{K-1})$ was equal to that obtained by the oracle (to 3
significant figures). The orange and pink lines show the result for Algorithm
1 using single- and double-precision, respectively. The purple and brown lines
show the result for computing $C(\boldsymbol{\eta})$ using the inverse Laplace
transform, i.e., Eq. 21, as discussed in Section 2.2. At each level of
$\sigma$, we show error bars over 10 random draws of $\boldsymbol{\eta}$.

###  2.2 Inverse laplace transform

In this Section, we show that $C(\boldsymbol{\eta})$ can be written as the
inverse Laplace transform of a function that does not suffer from catastrophic
cancellation. In particular, said function can be passed to Laplace inversion
methods (Davies and Martin, 1979) in order to evaluate $C(\boldsymbol{\eta})$
in the regime where Algorithm 1 fails.

Proposition: For a function $f:\mathbb{R}^{+}\to\mathbb{R}^{+}$, let
$\mathcal{L}[f](s)=\int_{0}^{\infty}f(t)e^{-st}dt$ denote the Laplace
transform. Define the function $c:\mathbb{R}^{+}\to\mathbb{R}^{+}$ by:

$\displaystyle c(t)$
$\displaystyle=\int_{\\{\mathbb{x}:\sum_{i=1}^{K-1}x_{i}\leq
t\\}}\prod_{i=1}^{K}e^{\eta_{i}x_{i}}d\mu(\mathbb{x}).$ (14)

Then the Laplace transform of $c$ is equal to:

$\displaystyle\mathcal{L}[c](s)=\prod_{i=1}^{K}\frac{1}{s-\eta_{i}}.$ (15)

Remark: The function $c$ includes the normalizing constant of the CC as a
special case $C(\boldsymbol{\eta})=c(1)$. More generally, we have that:

$\displaystyle c(t)=t^{K-1}C(\boldsymbol{\eta}/t)^{-1},$ (16)

as can be seen from letting $\tilde{x}_{i}=x_{i}/t$ in the integral 14.

Proof: Let $f_{\eta}(x)=e^{\eta x}$, so that:

$\displaystyle c(t)$
$\displaystyle=\int_{\\{\mathbb{x}:\sum_{i=1}^{K-1}x_{i}\leq
t\\}}\prod_{i=1}^{K}f_{\eta_{i}}(x_{i})d\mu(\mathbb{x})$
$\displaystyle=\int_{0}^{t}\int_{0}^{t-x_{1}}\cdots\int_{0}^{t-x_{1}-\cdots-
x_{K-2}}\prod_{i=1}^{K}f_{\eta_{i}}(x_{i})dx_{K-1}\cdots dx_{2}dx_{1}.$ (17)

Next we apply the transformation from Wolpert and Wolf (1995), defined as
$w_{1}=t$, $w_{k}=w_{k-1}-x_{k-1}$, or equivalently,
$w_{k}=t-\sum_{i=1}^{k-1}x_{i}$. We can then write our integral as a
convolution:

$\displaystyle
c(t)=\int_{0}^{w_{1}}\int_{0}^{w_{2}}\cdots\int_{0}^{w_{K-1}}\prod_{i=1}^{K}f_{\eta_{i}}(x_{i})dx_{K-1}\cdots
dx_{2}dx_{1}$
$\displaystyle=\int_{0}^{w_{1}}\int_{0}^{w_{2}}\cdots\int_{0}^{w_{K-2}}\prod_{i=1}^{K-2}f_{\eta_{i}}(x_{i})\int_{0}^{w_{K-1}}f_{\eta_{K-1}}(x_{K-1})f_{\eta_{K}}(w_{K-1}-x_{K-1})dx_{K-1}\cdots
dx_{2}dx_{1}$
$\displaystyle=\int_{0}^{w_{1}}\int_{0}^{w_{2}}\cdots\int_{0}^{w_{K-2}}\prod_{i=1}^{K-2}f_{\eta_{i}}(x_{i})(f_{\eta_{K-1}}\ast
f_{\eta_{K}})(w_{K-1})dx_{K-2}\cdots dx_{2}dx_{1}$
$\displaystyle=\int_{0}^{w_{1}}\cdots\int_{0}^{w_{K-3}}\prod_{i=1}^{K-3}f_{\eta_{i}}(x_{i})\int_{0}^{w_{K-2}}f_{\eta_{K-2}}(x_{K-2})(f_{\eta_{K-1}}\ast
f_{\eta_{K}})(w_{K-2}-x_{K-2})dx_{K-2}\cdots dx_{1}$
$\displaystyle=\int_{0}^{w_{1}}\int_{0}^{w_{2}}\cdots\int_{0}^{w_{K-3}}\prod_{i=1}^{K-3}f_{\eta_{i}}(x_{i})(f_{\eta_{K-2}}\ast
f_{\eta_{K-1}}\ast f_{\eta_{K}})(w_{K-2})dx_{K-3}\cdots dx_{2}dx_{1}$
$\displaystyle=\cdots$ $\displaystyle=(\circledast_{i=1}^{K}f_{\eta_{i}})(t).$
(18)

Next, since the Laplace transform of a convolution equals the product of the
Laplace transforms, we have that:

$\displaystyle\mathcal{L}[c](s)=\prod_{i=1}^{K}\mathcal{L}[f_{i}](s),$ (19)

but the univariate case is simply:

$\displaystyle\mathcal{L}[f_{i}](s)=\int_{0}^{\infty}e^{(\eta_{i}-s)t}dt=\frac{1}{s-\eta_{i}},$
(20)

and the result follows. ∎

Corollary: The normalizing constant of the continuous categorical distribution
can be written as the following inverse Laplace transform:

$\displaystyle
C(\boldsymbol{\eta})=\mathcal{L}^{-1}\left[\prod_{i=1}^{K}\frac{1}{s-\eta_{i}}\right](1).$
(21)

Proof: Take the inverse Laplace transform in Eq. 15 to find:

$\displaystyle
c(t)=\mathcal{L}^{-1}\left[\prod_{i=1}^{K}\frac{1}{s-\eta_{i}}\right](t).$
(22)

Taking $t=1$ gives the desired result. ∎

Unlike Eq. 9, the product in Eq. 15 does not suffer from catastrophic
cancellation, nor does it diverge whenever $\eta_{j_{1}}\approx\eta_{j_{2}}$
for some $j_{1}\neq j_{2}$. The corresponding Laplace inversion, i.e., Eq. 21,
provides an alternative method to compute our normalizing constant
$C(\boldsymbol{\eta})$. Numerous numerical algorithms exist for inverting the
Laplace transform; see (Cohen, 2007) for a survey. We note, however, that
inverting the Laplace transform is generally a hard problem (Epstein and
Schotland, 2008).

Empirically, we found some modest empiricasuccess in computing Eq. 21
numerically. We tested three inversion algorithms, due to Talbot (Talbot,
1979), Stehfest (Stehfest, 1970), and De Hoog (De Hoog et al., 1982). The
experimental setup was identical to that of Section 2.1.1, and the results are
incorporated into Figure 2. We found De Hoog’s method to be the most effective
on our problem, whereas Talbot’s always failed and is omitted from the Figure.
Importantly, De Hoog’s method showed some success in the regime where
Algorithm 1 failed, meaning it could be used as a complementary technique.
However, no inversion method was able to scale beyond $K=30$.

###  2.3 Inductive approach

In this section, we provide an alternative algorithm to compute our
normalizing constant. This algorithm will be based on the following recursive
property of $C(\boldsymbol{\eta})$, which was implicitly used as part of the
proof of Eq. 9 (Gordon-Rodriguez et al., 2020a).

Proposition: Define the subvector notation
$\boldsymbol{\eta}_{:k}=(\eta_{1},\dots,\eta_{k-1})$, and make the dependence
on $K$ explicit by writing
$C_{K}(\boldsymbol{\eta})=C(\eta_{1},\dots,\eta_{K-1})$. We also use the
notation
$\boldsymbol{\eta}-\eta_{k}=(\eta_{1}-\eta_{k},\dots,\eta_{K-1}-\eta_{k})$.
Then:

$\displaystyle
C_{K}(\boldsymbol{\eta})=\frac{e^{\eta_{K-1}}C_{K-1}\left(\boldsymbol{\eta}_{:(K-1)}-\eta_{K-1}\right)-C_{K-1}\left(\boldsymbol{\eta}_{:(K-1)}\right)}{\eta_{K-1}}.$
(23)

Proof: We start from the integral definition of the left hand side:

$\displaystyle C_{K}(\text{\boldmath$\eta$})$
$\displaystyle=\int_{{\mathbb{S}^{K-1}}}e^{\boldsymbol{\eta}^{\top}\mathbb{x}}d\mu$
$\displaystyle=\int_{0}^{1}\int_{0}^{1-x_{1}}\cdots\int_{0}^{1-x_{1}-\cdots-
x_{K-2}}e^{\sum_{i=1}^{K-1}\eta_{i}x_{i}}dx_{K-1}\cdots dx_{2}dx_{1}.$ (24)

For the innermost integral, we have:

$\displaystyle\int_{0}^{1-x_{1}-\cdots-
x_{K-2}}e^{\sum_{i=1}^{K-1}\eta_{i}x_{i}}dx_{K-1}$
$\displaystyle=e^{\sum_{i=1}^{K-2}\eta_{i}x_{i}}\int_{0}^{1-x_{1}-\cdots-
x_{K-2}}e^{\eta_{K-1}x_{K-1}}dx_{K-1}$
$\displaystyle=e^{\sum_{i=1}^{K-2}\eta_{i}x_{i}}\left(\frac{e^{\eta_{K-1}(1-x_{1}-\cdots-
x_{K-2})}-1}{\eta_{K-1}}\right)$
$\displaystyle=\frac{e^{\eta_{K-1}}e^{\sum_{i=1}^{K-2}(\eta_{i}-\eta_{K-1})x_{i}}-e^{\sum_{i=1}^{K-2}\eta_{i}x_{i}}}{\eta_{K-1}}$
$\displaystyle=\frac{e^{\eta_{K-1}}e^{(\boldsymbol{\eta}_{:(K-1)}-\eta_{K-1})^{\top}\mathbb{x}_{:(K-1)}}-e^{\boldsymbol{\eta}_{:(K-1)}^{\top}\mathbb{x}_{:(K-1)}}}{\eta_{K-1}},$
(25)

and the result follows by linearity of the integral. ∎

Remark: The base case $K=2$ corresponds to the univariate continuous Bernoulli
distribution, which admits a straightforward Taylor expansion that is useful
around the unstable region $\eta\approx 0$ (Loaiza-Ganem and Cunningham,
2019):

$\displaystyle
C_{2}(\boldsymbol{\eta}_{:2})=\frac{e^{\eta_{1}}-1}{\eta_{1}}=\frac{(1+\eta_{1}+\frac{1}{2!}\eta_{1}^{2}+\cdots)-1}{\eta_{1}}=1+\frac{1}{2!}\eta_{1}+\cdots.$
(26)

In words, Eq. 23 is stating that we can compute $C(\boldsymbol{\eta})$ for the
full parameter vector $\boldsymbol{\eta}=(\eta_{1},\dots,\eta_{K-1})$ by first
computing it for the lower-dimensional parameter vectors:

$\displaystyle\boldsymbol{\eta}_{:K-1}$
$\displaystyle=(\eta_{1},\dots,\eta_{K-2}),$
$\displaystyle\boldsymbol{\eta}_{:K-1}-\eta_{K-1}$
$\displaystyle=(\eta_{1}-\eta_{K-1},\dots,\eta_{K-2}-\eta_{K-1}).$

Substituting these back into Eq. 23, we have that:

$\displaystyle C_{K-1}(\boldsymbol{\eta}_{:K-1})$
$\displaystyle=\frac{e^{\eta_{K-2}}C_{K-2}(\boldsymbol{\eta}_{1:K-2}-\eta_{K-2})-C_{K-2}(\boldsymbol{\eta}_{1:K-2})}{\eta_{K-2}},$
$\displaystyle C_{K-1}(\boldsymbol{\eta}_{:K-1}-\eta_{K-1})$
$\displaystyle=\frac{e^{\eta_{K-2}-\eta_{K-1}}C_{K-2}(\boldsymbol{\eta}_{1:K-2}-\eta_{K-2})-C_{K-2}(\boldsymbol{\eta}_{1:K-2}-\eta_{K-1})}{\eta_{K-2}-\eta_{K-1}}.$

Note that we are now left with not 4, but 3 new parameter vectors to recurse
on: $\boldsymbol{\eta}_{1:K-3}$, $\boldsymbol{\eta}_{1:K-3}-\eta_{K-1}$, and
$\boldsymbol{\eta}_{1:K-3}-\eta_{K-2}$. Repeating the argument $K-2$ times and
working backwards we obtain Algorithm 2 for computing the normalizing
constant.

Algorithm 2 Inductive computation of the normalizing constant

Input: A parameter vector $\boldsymbol{\eta}$
Output: The normalizing constant $C(\boldsymbol{\eta})$

1:Initialize $\mathbb{c}=(1,\dots,1)\in\mathbb{R}^{K-1}$ and set
$\mathbb{\tilde{c}}=\mathbb{c}$.

2:for $k=1,2,\dots,K-1$ do

3: for $i=1,\dots,K-k$ do

4: Set $\xi_{i}=\eta_{k}-\eta_{k+i}$

5: Set $\tilde{c}_{i}=\frac{e^{\xi_{i}}c_{1}-c_{i+1}}{\xi_{i}}$

6: end for

7: Set $\mathbb{c}=\mathbb{\tilde{c}}$

8:end for

9:return $c_{1}$

We find the numerical properties of Algorithm 2 to perform similarly to
Algorithm 1, suffering from the same cancellation issues in high dimensions.
Nevertheless, we hope this alternative scheme may help to inspire further
numerical improvements. For example, since the floating-point behavior of
Algorithm 2 depends on the ordering of the elements of $\boldsymbol{\eta}$, it
may be possible do design a reordering scheme that improves the overall
precision of the algorithm; such reorderings have been explored in the context
of sampling algorithms for the CC (Gordon-Rodriguez et al., 2020a). Other
possibilities include Kahan summation (Kahan, 1965) or the compensated Horner
algorithm (Langlois and Louvet, 2007); we leave their study to future work.

###  2.4 Overflow

We conclude this section with a small remark on numerical overflow in the
context of Algorithm 1. In high dimensions and with high $\eta$ values,
overflow can occur for the terms in the summand, due to a large $e^{\eta_{j}}$
term or a small denominator. This can be addressed by re-writing the
normalizing constant in a form that allows us to take advantage of the log-
sum-exp trick:

$\displaystyle\log C(\text{\boldmath$\eta$})$
$\displaystyle=\log\left(\sum_{k=1}^{K}\frac{{e^{\eta_{k}}}}{\prod_{i\neq
k}\left(\eta_{k}-\eta_{i}\right)}\right)$
$\displaystyle=\log\left(\sum_{k=1}^{K}\frac{{e^{\eta_{k}}}}{\mathrm{sign}\left(\prod_{i\neq
k}\left(\eta_{k}-\eta_{i}\right)\right)\prod_{i\neq
k}|\eta_{k}-\eta_{i}|}\right)$
$\displaystyle=\log\left(\sum_{k=1}^{K}\mathrm{sign}\left(\prod_{i\neq
k}\left(\eta_{k}-\eta_{i}\right)\right)\exp\left(\eta_{k}-\sum_{i\neq
k}\log\left|\eta_{k}-\eta_{i}\right|\right)\right).$

##  3 Normalizing constant with repeated parameters

Whenever we have an equality between any pair of parameters, Eq. 9 is
undefined and, indeed, its proof by Gordon-Rodriguez et al. (2020a) breaks
down. In this Section, we derive a counterpart to Eq. 9 for the case when 2 or
more elements of $\boldsymbol{\eta}$ are equal to one another. Note that, for
mathematical convenience, we shall now denote
$\boldsymbol{\eta}\in\mathbb{R}^{K}$, where the $K$th component $\eta_{K}$ is
now included in the vector $\boldsymbol{\eta}$. As discussed in Section 1,
this component can be taken as fixed at 0, or it can be treated as an
additional free parameter, resulting in an overparameterized model (our
results will remain correct nevertheless).

###  3.1 A simple example

First, we illustrate the main idea of our argument using an example with $K=3$
and $\eta_{1}=\eta_{2}\neq\eta_{3}=0$.444Note that the case $\eta_{3}\neq 0$
and the case $\eta_{1}\neq\eta_{2}=\eta_{3}$ are mathematically equivalent
(albeit more algebraically cumbersome), since we can permute the elements of
$\boldsymbol{\eta}$ and shift by a constant without loss of generality. By
definition, the normalizing constant is then:

$\displaystyle
C_{3}(\boldsymbol{\eta})=\int_{0}^{1}\int_{0}^{1-x_{1}}e^{\eta_{1}x_{1}+\eta_{2}x_{2}}dx_{2}dx_{1}=\int_{0}^{1}\int_{0}^{1-x_{1}}e^{\eta_{1}(x_{1}+x_{2})}dx_{2}dx_{1}.$
(27)

We apply the change of variables $u=x_{1}+x_{2}$, $v=x_{1}$ to obtain:

$\displaystyle
C_{3}(\boldsymbol{\eta})=\int_{0}^{1}\int_{0}^{u}e^{\eta_{1}u}dvdu=\int_{0}^{1}ue^{\eta_{1}u}du.$
(28)

Note that the last integral corresponds to the expectation of a univariate CC
random variable, an idea we now generalize.

###  3.2 General formula

We start by proving the following lemma, which relates $C(\boldsymbol{\eta})$,
where some elements of $\boldsymbol{\eta}$ are repeated (potentially many
times), to $C(\boldsymbol{\eta}^{\prime})$, where $\boldsymbol{\eta}^{\prime}$
collapses the repeated elements of $\boldsymbol{\eta}$ onto a single
coordinate. We again use subvector notation
$\mathbb{x}_{k:}=(x_{k},x_{k+1},\dots,x_{K-1})$.

Lemma: Let $\mathbb{x}\sim\mathcal{CC}_{K}(\boldsymbol{\eta})$ with
$\eta_{1}=\eta_{2}=\dots=\eta_{k}$, for $1\leq k\leq K-1$, and let
$f:\mathbb{R}^{K-k-1}\to\mathbb{R}$. Then:

${C_{K}(\boldsymbol{\eta})}\mathbb{E}_{\mathbb{x}\sim\mathcal{CC}_{K}(\boldsymbol{\eta})}[f(\mathbb{x}_{(k+1):})]={C_{K-k+1}(\boldsymbol{\eta}_{k:})}\mathbb{E}_{\mathbb{u}\sim\mathcal{CC}_{K-k+1}(\boldsymbol{\eta}_{k:})}\left[\dfrac{u_{1}^{k-1}}{(k-1)!}f(\mathbb{u}_{2:})\right].$
(29)

Remark: We are not assuming that the last $K-k$ coordinates of
$\boldsymbol{\eta}$ are all different, we are simply assuming that the first
$k$ are identical. Note also that the positions of the coordinates could be
arbitrary and need not be the first $k$ ones; we can always use this lemma to
collapse repeated parameter values provided that the function $f$ does not
depend on the corresponding coordinates (we can simply relabel the coordinates
by a suitable permutation without loss of generality).

Proof: By definition of the left hand side:

$\displaystyle C_{K}$
$\displaystyle(\boldsymbol{\eta})\mathbb{E}_{\mathbb{x}\sim\mathcal{CC}_{K}(\boldsymbol{\eta})}[f(\mathbb{x}_{(k+1):})]=\int_{0}^{1}\int_{0}^{1-x_{1}}\cdots\int_{0}^{1-x_{1}-\cdots-
x_{K-2}}f(\mathbb{x}_{(k+1):})e^{\boldsymbol{\eta}^{\top}\mathbb{x}}dx_{K-1}\cdots
dx_{2}dx_{1}$
$\displaystyle=\int_{0}^{1}\int_{0}^{1-x_{1}}\cdots\int_{0}^{1-x_{1}-\cdots-
x_{K-2}}f(\mathbb{x}_{(k+1):})e^{\eta_{k}(x_{1}+\dots+x_{k})+\boldsymbol{\eta}_{(k+1):}^{\top}\mathbb{x}_{(k+1):}}dx_{K-1}\cdots
dx_{2}dx_{1}.$ (30)

Consider the following change of variable (note this is just like Section 3.1,
but with more bookkeeping):

$\begin{cases}u_{1}=x_{1}+x_{2}+\cdots+x_{k}\\\ u_{2}=x_{k+1}\\\
u_{3}=x_{k+2}\\\ \hskip 16.0pt\vdots\\\ u_{K-k}=x_{K-1}\\\ v_{1}=x_{1}\\\
v_{2}=x_{2}\\\ \hskip 16.0pt\vdots\\\ v_{k-1}=x_{k-1}\end{cases}$ (31)

This change of variable amounts to an invertible linear transformation with
the property that the absolute value of the determinant of its Jacobian is
$1$, so that we have:

$\displaystyle C_{K}(\boldsymbol{\eta})$
$\displaystyle\mathbb{E}_{\mathbb{x}\sim\mathcal{CC}_{K}(\boldsymbol{\eta})}[f(\mathbb{x}_{(k+1):})]=\int_{0}^{1}\int_{0}^{1-u_{1}}\cdots\int_{0}^{1-u_{1}-u_{2}-\cdots-
u_{K-k}}\int_{0}^{u_{1}}\int_{0}^{u_{1}-v_{1}}\cdots$ (32)
$\displaystyle\dots\int_{0}^{u_{1}-v_{1}-v_{2}-\dots-
v_{k-2}}f(\mathbb{u}_{2:})e^{\eta_{k}u_{1}+\boldsymbol{\eta}_{(k+1):}^{\top}\mathbb{u}_{2:}}dv_{k-1}\cdots
dv_{2}dv_{1}du_{K-k}\cdots du_{2}du_{1}.$

Note that the change of variables is such that the integrand does not depend
on $v_{1},v_{2},\dots,v_{k-1}$. Therefore:

$\displaystyle C_{K}(\boldsymbol{\eta})$
$\displaystyle\mathbb{E}_{\mathbb{x}\sim\mathcal{CC}_{K}(\boldsymbol{\eta})}[f(\mathbb{x}_{(k+1):})]$
(33)
$\displaystyle=\int_{0}^{1}\int_{0}^{1-u_{1}}\cdots\int_{0}^{1-u_{1}-u_{2}-\cdots-
u_{K-k}}g(u_{1})f(\mathbb{u}_{2:})e^{\boldsymbol{\eta}_{k:}^{\top}\mathbb{u}}du_{K-k}\cdots
du_{2}du_{1},$

where:

$\displaystyle
g(u_{1})=\int_{0}^{u_{1}}\int_{0}^{u_{1}-v_{1}}\cdots\int_{0}^{u_{1}-v_{1}-v_{2}-\cdots-
v_{k-2}}dv_{k-1}\cdots dv_{2}dv_{1}.$ (34)

But this is simply the Lebesgue measure of a simplex inscribed in the
hypercube $[0,u_{1}]^{K-1}$, so that
$g(u_{1})=u_{1}^{k-1}\mu(\mathbb{S}^{K})=u_{1}^{k-1}/(k-1)!$ (this can also be
seen by changing variables to $\tilde{v}_{i}=v_{i}/u_{i}$, or by applying Eq.
16). Multiplying and dividing by $C_{K-k+1}(\boldsymbol{\eta}_{k:})$ gives the
desired result. ∎

We can now derive a formula for the normalizing constant for an arbitrary
parameter vector $\boldsymbol{\eta}$.

Corollary: Let $\boldsymbol{\eta}\in\mathbb{R}^{K}$ contain $D\leq K$ unique
elements. Assume, without loss of generality, that
$\boldsymbol{\eta}=(\eta_{1},\dots,\eta_{1},\eta_{2},\dots,\eta_{2},\dots,\eta_{D},\dots,\eta_{D})$,
where each coordinate is repeated $1\leq r_{i}\leq K$ times, with
$\sum_{i=1}^{D}r_{i}=K$. Then:

$C_{K}(\boldsymbol{\eta})={C_{D}(\eta_{1},\eta_{2},\dots,\eta_{D})}{\mathbb{E}_{\mathbb{u}\sim\mathcal{CC}_{D}(\eta_{1},\eta_{2},\dots,\eta_{D})}\left[\displaystyle\prod_{i=1}^{D}\dfrac{u_{i}^{r_{i}-1}}{(r_{i}-1)!}\right]}$
(35)

Proof: The result follows by applying the lemma $D$ times. First, we apply the
lemma with $f(\cdot)=1$:

$\displaystyle
C_{K}(\boldsymbol{\eta})={C_{K-r_{1}+1}(\boldsymbol{\eta}^{\prime})}{\mathbb{E}_{\mathbb{u}\sim\mathcal{CC}_{K-r_{1}+1}(\boldsymbol{\eta}^{\prime})}\left[\dfrac{u_{1}^{r_{1}-1}}{(r_{1}-1)!}\right]},$
(36)

where
$\boldsymbol{\eta}^{\prime}=\boldsymbol{\eta}_{r_{1}:}=(\eta_{1},\eta_{2},\dots,\eta_{2},\dots,\eta_{D},\dots,\eta_{D})$,
i.e., we have collapsed the first parameter value onto a single coordinate.
Next, we apply the lemma a second time on the new expectation term to collapse
the $\eta_{2}$ values, this time using $f(u_{1})=u_{1}^{r_{1}-1}/(r_{1}-1)!$,
which does not depend on the $\eta_{2}$ coordinates:

$\displaystyle{C_{K-r_{1}+1}(\boldsymbol{\eta}^{\prime})}\mathbb{E}$
${}_{\mathbb{u}\sim\mathcal{CC}_{K-r_{1}+1}(\boldsymbol{\eta}^{\prime})}\left[\dfrac{u_{1}^{r_{1}-1}}{(r_{1}-1)!}\right]$
(37)
$\displaystyle={C_{K-r_{1}-r_{2}+2}(\boldsymbol{\eta}^{\prime\prime})}\mathbb{E}_{\mathbb{u}\sim\mathcal{CC}_{K-r_{1}-r_{2}+2}(\boldsymbol{\eta}^{\prime\prime})}\left[\dfrac{u_{1}^{r_{1}-1}}{(r_{1}-1)!}\dfrac{u_{2}^{r_{2}-1}}{(r_{2}-1)!}\right],$

where
$\boldsymbol{\eta}^{\prime\prime}=(\eta_{1},\eta_{2},\eta_{3},\dots,\eta_{3},\dots,\eta_{D},\dots,\eta_{D})$.
Repeating $D$ times yields the desired result. ∎

Note that Eq. 35 can be computed using known results. The term
$C_{D}(\eta_{1},\dots,\eta_{D})$ can be evaluated with Eq. 9, since all the
parameter values are distinct. The expectation term can be computed by
differentiating the moment generating function of
$\mathbb{u}\sim\mathcal{CC}_{D}(\eta_{1},\eta_{2},\dots,\eta_{D})$, as
discussed by Gordon-Rodriguez et al. (2020a). Note that in real data
applications exact equality between parameter values may never occur, and it
is unclear how close the elements of $\boldsymbol{\eta}$ should be in order to
warrant applying Eq. 35. Nevertheless, Eq. 35 is important for theoretical
completeness.

##  4 Discussion

The normalizing constant of the CC distribution is essential to applications,
being a necessary prerequisite for evaluating likelihoods, optimizing models,
and simulating samples alike. Computing this normalizing constant is
nontrivial, and doing so in high dimensions remains an open problem. Our work
represents a significant step toward this goal, improving our understanding of
the numerical properties of different computation techniques, as well as
advancing the underlying theory and algorithms. In addition, we hope our
results will help to inspire further advances and to develop increasingly
robust numerical techniques that will ultimately enable the use of the CC
distribution with arbitrary parameter values on high-dimensional problems.

## References

  * Bond-Taylor et al. (2021) Sam Bond-Taylor, Adam Leach, Yang Long, and Chris G Willcocks.  Deep generative modelling: A comparative review of vaes, gans, normalizing flows, energy-based and autoregressive models.  _arXiv preprint arXiv:2103.04922_ , 2021.
  * Cohen (2007) Alan M Cohen.  _Numerical methods for Laplace transform inversion_ , volume 5.  Springer Science & Business Media, 2007.
  * Davies and Martin (1979) Brian Davies and Brian Martin.  Numerical inversion of the laplace transform: a survey and comparison of methods.  _Journal of computational physics_ , 33(1):1–32, 1979.
  * De Hoog et al. (1982) Frank R De Hoog, JH Knight, and AN Stokes.  An improved method for numerical inversion of laplace transforms.  _SIAM Journal on Scientific and Statistical Computing_ , 3(3):357–366, 1982.
  * Dirichlet (1839) Peter Gustav Lejeune Dirichlet.  Sur une nouvelle méthode pour la détermination des intégrales multiples.  _Journal de Mathématiques, Ser I, 4_ , pages 164–168, 1839.
  * Epstein and Schotland (2008) Charles L Epstein and John Schotland.  The bad truth about laplace’s transform.  _SIAM review_ , 50(3):504–520, 2008.
  * Goldberg (1991) David Goldberg.  What every computer scientist should know about floating-point arithmetic.  _ACM computing surveys (CSUR)_ , 23(1):5–48, 1991\.
  * Gordon-Rodriguez et al. (2020a) Elliott Gordon-Rodriguez, Gabriel Loaiza-Ganem, and John Cunningham.  The continuous categorical: a novel simplex-valued exponential family.  In _International Conference on Machine Learning_ , pages 3637–3647. PMLR, 2020a.
  * Gordon-Rodriguez et al. (2020b) Elliott Gordon-Rodriguez, Gabriel Loaiza-Ganem, Geoff Pleiss, and John Patrick Cunningham.  Uses and abuses of the cross-entropy loss: Case studies in modern deep learning.  In _Proceedings on "I Can’t Believe It’s Not Better!" at NeurIPS Workshops_ , volume 137 of _Proceedings of Machine Learning Research_ , pages 1–10. PMLR, 2020b.
  * Johansson et al. (2013) Fredrik Johansson et al.  _mpmath: a Python library for arbitrary-precision floating-point arithmetic (version 0.18)_ , December 2013.  http://mpmath.org/.
  * Kahan (1965) William Kahan.  Pracniques: further remarks on reducing truncation errors.  _Communications of the ACM_ , 8(1):40, 1965.
  * Kingma and Welling (2014) Diederik P. Kingma and Max Welling.  Auto-encoding variational bayes.  In _2nd International Conference on Learning Representations, ICLR 2014, Banff, AB, Canada, April 14-16, 2014, Conference Track Proceedings_ , 2014.
  * Langlois and Louvet (2007) Philippe Langlois and Nicolas Louvet.  How to ensure a faithful polynomial evaluation with the compensated horner algorithm.  In _18th IEEE Symposium on Computer Arithmetic (ARITH’07)_ , pages 141–149. IEEE, 2007.
  * Loaiza-Ganem and Cunningham (2019) Gabriel Loaiza-Ganem and John P Cunningham.  The continuous bernoulli: fixing a pervasive error in variational autoencoders.  In _Advances in Neural Information Processing Systems_ , pages 13266–13276, 2019.
  * Stehfest (1970) Harald Stehfest.  Algorithm 368: Numerical inversion of laplace transforms [d5].  _Communications of the ACM_ , 13(1):47–49, 1970\.
  * Talbot (1979) Alan Talbot.  The accurate numerical inversion of laplace transforms.  _IMA Journal of Applied Mathematics_ , 23(1):97–120, 1979.
  * Wolpert and Wolf (1995) David H Wolpert and David R Wolf.  Estimating functions of probability distributions from a finite set of samples.  _Physical Review E_ , 52(6):6841, 1995.