# Observation of magnon cross-Kerr effect in cavity magnonics

Wei-Jiang Wu  Interdisciplinary Center of Quantum Information and Zhejiang
Province Key Laboratory of Quantum Technology and Device, Department of
Physics and State Key Laboratory of Modern Optical Instrumentation, Zhejiang
University, Hangzhou 310027, China  Da Xu<EMAIL_ADDRESS>Interdisciplinary
Center of Quantum Information and Zhejiang Province Key Laboratory of Quantum
Technology and Device, Department of Physics and State Key Laboratory of
Modern Optical Instrumentation, Zhejiang University, Hangzhou 310027, China
<EMAIL_ADDRESS>Jie Qian  Interdisciplinary Center of Quantum Information and
Zhejiang Province Key Laboratory of Quantum Technology and Device, Department
of Physics and State Key Laboratory of Modern Optical Instrumentation,
Zhejiang University, Hangzhou 310027, China  Jie Li  Interdisciplinary Center
of Quantum Information and Zhejiang Province Key Laboratory of Quantum
Technology and Device, Department of Physics and State Key Laboratory of
Modern Optical Instrumentation, Zhejiang University, Hangzhou 310027, China
Yi-Pu Wang<EMAIL_ADDRESS>Interdisciplinary Center of Quantum
Information and Zhejiang Province Key Laboratory of Quantum Technology and
Device, Department of Physics and State Key Laboratory of Modern Optical
Instrumentation, Zhejiang University, Hangzhou 310027, China  J. Q. You
<EMAIL_ADDRESS>Interdisciplinary Center of Quantum Information and Zhejiang
Province Key Laboratory of Quantum Technology and Device, Department of
Physics and State Key Laboratory of Modern Optical Instrumentation, Zhejiang
University, Hangzhou 310027, China

###### Abstract

When there is a certain amount of field inhomogeneity, the biased
ferrimagnetic crystal will exhibit the higher-order magnetostatic (HMS) mode
in addition to the uniform-precession Kittel mode. In cavity magnonics, we
show both experimentally and theoretically the cross-Kerr-type interaction
between the Kittel mode and HMS mode. When the Kittel mode is driven to
generate a certain number of excitations, the HMS mode displays a
corresponding frequency shift, and vice versa. The cross-Kerr effect is caused
by an exchange interaction between these two spin-wave modes. Utilizing the
cross-Kerr effect, we realize and integrate a multi-mode cavity magnonic
system with only one yttrium iron garnet (YIG) sphere. Our results will bring
new methods to magnetization dynamics studies and pave a way for novel cavity
magnonic devices by including the magnetostatic mode-mode interaction as an
operational degree of freedom.

##  I INTRODUCTION

Nonlinear effects are a large class of physical phenomena that have been
continuously and actively explored in optics Boyd ; Autere2018 ; Zhang2019 ,
plasmonics Kauranen2012 ; Mesch2016 ; Zhong2020 , and acoustics Hamilton1998 ;
Fang2017 ; Darabi2019 , as well as in other oscillations and wave systems
Nayfeh2008 ; Kartashov2019 ; Zuo2020 . Nonlinearity and nonlinear interactions
are weak and difficult to be detected in some systems, but they can also
become sufficiently strong to be the dominant factor in some other systems.
The cross-Kerr effect is one of the subtle nonlinear interactions between
fields and waves that can exist in natural ions Ding2017 , atoms He2014 ;
Xia2018 ; Sinclair2019 , and artificial mesoscopic architectures, such as
superconducting circuits Hoi2013 ; Kounalakis2018 ; Vrajitoarea2020 . The
cross-Kerr effect is a nonlinear shift in the frequency of a resonator as a
function of the number of excitations in another mode that interacts with the
resonator. Understanding these nonlinear interactions is not only of
fundamental importance, but also useful in various applications. In practice,
the cross-Kerr effect can be utilized to implement quantum logic gates
Turchette1995 ; Brod2016 , prepare entangled photons Lukin2000 ; Sheng2008 ,
and carry out quantum non-demolition measurements Wang2019 ; Dassonneville2020
.

Cavity magnonics, on the other hand, has gradually demonstrated its unique
advantages in fundamental and applied research over the last few years
Huebl2013 ; Tabuchi2014 ; Zhang2014 ; Goryachev2014 ; Bai2015 ; Cao2015 ;
Zhang2015 ; Tabuchi2015 ; Haigh2016 ; Maier-Flaig2016 ; A.Osada2016 ;
X.Zhang2018 ; Zhang2017 ; PuWang2018 ; A.Osada2018 ; JieLi2018 ; Harder2018 ;
Y.P.Wang2019 ; Lachance-Quirion2020 ; Y.Yang2020 ; J.Xu2020 ; S.P.Wolski2020 ;
Jinwei2021 ; Xufeng2021 ; Potts2021 ; D.Zhang2015 . It is expected to be a
critical component of hybrid quantum systems LachanceQuirion2019  and quantum
network nodes Jie2021 . The commonly used physical implementation of a cavity
magnonic system consists of a microwave cavity and ferro-(ferri-)magnets, such
as YIG. The spin collective excitation modes are spin waves in the
ferrimagnetic spin ensemble. The most widely studied spin wave mode is the
uniform precession mode of spins (i.e., the Kittel mode), which has the
largest magnetic dipole moment and can strongly interact with the microwave in
both the quantum limit Huebl2013 ; Tabuchi2014  and room-temperature classical
regime Zhang2014 . Under the bias of an external magnetic field with certain
inhomogeneity, besides the Kittel mode, there are also higher-order
magnetostatic (HMS) modes in the spin ensemble J.F.Dillon1958 ;
P.C.Fletcher1959 . Previous studies often ignored the influence of these
magnetostatic modes on the observation and application of Kittel mode in
experiments. However, the HMS modes actually have nontrivial spin textures.
Their non-zero orbital angular momentum feature can be used to demonstrate
nonreciprocal Brillouin light scattering in optomagnonics and other chiral
optics A.Osada2018 ; J.Graf2018 ; J.A.Haigh2018 ; Osada_2018 . In the case of
YIG sphere supporting optical whispering-gallery modes, the optomagnonic
coupling strength can be enhanced for HMS modes due to their reduced mode
volume and the localized distribution on the boundary J.Graf2018 ;
A.Gloppe2019 , which increases the microwave-to-optical photon transduction
efficiency in the cavity optomagnonic system.

In this work, we focus on the cavity magnonic system with the cross-Kerr
interaction between the Kittel mode and HMS mode, used as a new control degree
of freedom. This cross-Kerr interaction is investigated both theoretically and
experimentally in a coupled 3D cavity and YIG sphere system, where the YIG
sphere supports both the Kittel mode and HMS mode. A large frequency detuning
between these two modes is found. When the Kittel mode or HMS mode is pumped
to generate a certain number of magnon excitations, the pumped mode gains a
frequency shift due to the self-Kerr effect PuWang2018 ; Y.P.Wang2016 ; GQ2019
; Kong2019 ; Xiao2021 ; Ruichang2021 ; Gurevich2020 . Meanwhile, the frequency
of the undriven mode also shifts due to the cross-Kerr effect between these
two modes. Then, the ratio of the self-Kerr and cross-Kerr coefficients can be
obtained. Furthermore, using the cross-Kerr coefficient as a ruler, we find
that the magnetostatic mode’s self-Kerr coefficient is greater than that of
the Kittel mode. This offers additional controllable degrees of freedom in
cavity magnonics without increasing system’s complexity. Our study will also
bring new ideas to the fundamental magnetostatics studies and hybrid magnonic
operations Y.Li2020 ; B.Z.R ; Yuan .

Figure 1: (a) Schematic diagram of the experimental setup. The YIG sphere is
coupled to the cavity mode $TE_{102}$. The cavity is placed in the bias
magnetic field $B_{0}$ (blue arrow). The microwave source connected to the
loop antenna provides a drive field perpendicular to the bias magnetic field.
The vector network analyzer (VNA) measures the cavity transmission spectrum
$|S_{21}|^{2}$. Inset: in the YIG sphere, the Kittel mode and HMS mode are
denoted by two sub-magnetizations. (b) Schematic diagram of the self-Kerr and
cross-Kerr effect. In the upper panel, the drive field applied on the Kittel
mode will excite a certain number of magnons and enlarge the precession angle
of the Kittel-mode magnetization $M_{K}$. The self-Kerr effect will cause a
frequency shift of the Kittel mode. The cross-Kerr effect will result in a
frequency shift of the HMS mode without changing the precession angle of the
HMS-mode magnetization $M_{H}$. A similar effect occurs by exchanging the
roles of the two modes, as shown in the lower panel.

##  II SYSTEM AND THEORETICAL MODEL

The experimental setup of the cavity magnonic system is illustrated in Fig.
1(a). It consists of the microwave field in a 3D copper cavity (internal
dimension is 44.0$\times$20.0$\times$6.0 mm3 ), the magnon modes in a single
crystal YIG sphere (1 mm in diameter), and the drive field provided by a
microwave (MW) source. The YIG sphere is glued on the top of a beryllium
ceramic rod and mounted at the antinode of the magnetic field of the cavity
mode $TE_{102}$ (${\omega_{c}}/2\pi=10.07$ GHz). The whole sample is placed in
a bias magnetic field $B_{0}$ applied along the $\left[{110}\right]$ crystal
axis, which results in a negative self-Kerr coefficient PuWang2018 ;
Gurevich2020 . The bias magnetic field, the magnetic field of the cavity mode
$TE_{102}$, and the drive field magnetic component are perpendicular to each
other at the position of the YIG sphere. This configuration maximizes the
cavity-YIG coupling strength and drive efficiency. The drive field is loaded
to port 3, which is terminated with a loop antenna. The vector network
analyzer (VNA) is used for measuring the microwave transmission spectrum
$|S_{21}|^{2}$ via port 1 and port 2. The probe power is much weaker than the
drive power. It should be noted that the loop antenna plays multiple roles:
(i) radiates the drive field; (ii) works as a dissipative channel for the
magnon modes; (iii) creates an inhomogeneous magnetic field, which gives rise
to the emergence of HMS modes.

Figure 2: Drive on the Kittel mode. (a) The transmission spectra are measured
versus magnet coil current. The Kittel mode (red dashed line) and HMS mode
(green dashed line) are coupled with the cavity mode. (b) The transmission
spectra are measured versus magnet coil current, while a drive field of 9.8
GHz and 25 dBm is applied on the Kittel mode. The Kittel mode has a frequency
shift $\Delta_{k,s}$ due to the self-Kerr effect, and the HMS mode has a even
larger frequency shift $\Delta_{h,c}$ caused by the cross-Kerr effect. (c) The
frequency shifts of the Kittel mode and HMS mode are extracted from (b), and
plotted versus $\delta_{k}=\omega_{k}-\omega_{d}$. The experimental results
are fitted by Eq. (2) and
${\Delta_{h,c}}=\frac{K_{cross}}{K_{k,s}}\Delta_{k,s}$, which also gives the
ratio $K_{cross}/K_{k,s}=2.5$.

Owing to the inhomogeneity of the magnetic field, apart from the Kittel mode,
the HMS mode is also excited. The two modes are both supported by the
collective spin motion. More spin moments contribute to the spin uniform
precession mode (Kittel mode), yielding a large dipole moment for the Kittel
mode. Therefore, the couplings between the Kittel mode and other
electromagnetic fields are usually strong. However, fewer spin moments
contribute to the dipole of the HMS mode, so the coupling between the HMS mode
and cavity mode is weak in cavity magnonics LachanceQuirion2019 ; C.Zhang2021
. Both the Kittel mode and HMS mode are nonlinear modes on account of the
magnetocrystalline anisotropy. The magnetocrystalline anisotropy energy
measures the energy cost when tuning the magnetization orientation from the
easy axis to the hard axis. Mathematically, the energy term is proportional to
$S_{z}^{2}$ ($S_{z}$ is the macrospin operator) and a self-Kerr term
$K_{s}m^{\dagger}mm^{\dagger}m$ can be obtained, where $K_{s}$ is the Kerr
coefficient and $m^{\dagger}$ ($m$) is the creation (annihilation) operator of
the magnon mode Y.P.Wang2016 ; GQ2019 . When the magnon mode is pumped to
generate a certain number of excitations, the precession angle of the
magnetization increases. This process leads to the change of the
magnetocrystalline anisotropy energy, which results in the frequency shift of
the nonlinear mode. In previous studies, the interaction between the Kittel
mode and HMS mode has rarely been investigated. Here, we model the two magnon
modes with two sub-magnetizations ($M_{K}$ and $M_{H}$), which nonlinearly
interact with each other. The interaction term is proportional to $S_{K}\cdot
S_{H}$ (subscript $K$ and $H$ represent the Kittel mode and HMS mode), which
gives a two-mode cross-Kerr term. The total Hamiltonian of the coupled system
is (see Appendix A)

$\begin{split}H/\hbar&=\omega_{c}a^{\dagger}a+\omega_{k}b^{\dagger}b+\omega_{h}c^{\dagger}c\\\
&+K_{k,s}b^{\dagger}bb^{\dagger}b+K_{h,s}c^{\dagger}cc^{\dagger}c+K_{cross}b^{\dagger}bc^{\dagger}c\\\
&+g_{k}\left(a^{\dagger}b+b^{\dagger}a\right)+g_{h}\left(a^{\dagger}c+c^{\dagger}a\right)+g_{k,h}\left(b^{\dagger}c+c^{\dagger}b\right)\\\
&+\Omega_{k}(b^{\dagger}e^{-i{\omega_{d}}t}+be^{i{\omega_{d}}t})+\Omega_{h}(c^{\dagger}e^{-i{\omega_{d}}t}+ce^{i{\omega_{d}}t})\end{split}$
(1)

where $a^{\dagger}$ ($a$), $b^{\dagger}$ ($b$) and $c^{\dagger}$ ($c$) are the
creation (annihilation) operators of the cavity mode at frequency
$\omega_{c}$, the Kittel mode at frequency $\omega_{k}$, and the HMS mode at
frequency $\omega_{h}$, respectively; $K_{k,s}$ and $K_{h,s}$ are the self-
Kerr coefficients of the Kittel mode and HMS mode; $K_{cross}$ is the cross-
Kerr coefficient between the two magnon modes; $g_{k}$ ($g_{h}$) is the
coupling strength between the Kittel (HMS) mode and cavity mode; $g_{k,h}$ is
the coherent coupling strength between the Kittel mode and HMS mode (see
Appendix A); ${\Omega_{k}}$ (${\Omega_{h}}$) is the drive-field strength on
the Kittel (HMS) mode, and $\omega_{d}$ is the drive-field frequency.

The observable effects of the self-Kerr and cross-Kerr effects are
schematically illustrated in Fig.1(b). The drive field applied on the Kittel
mode enlarges the precession angle of the magnetization $M_{K}$, which also
corresponds to the generation of a certain number of magnon excitations. In
the case of Kittel mode with a negative self-Kerr coefficient and under the
mean-field approximation, the frequency of the Kittel mode will have a red-
shift $\Delta_{k,s}$ PuWang2018 ; Ruichang2021 , equals to $2K_{k,s}\langle
b^{\dagger}b\rangle$. Here the extra subscript $s$ represent the self-Kerr.
Simultaneously, the cross-Kerr effect will also induce a frequency shift of
the undriven HMS mode, and the precession angle of the magnetization $M_{H}$
is, however, unchanged. The frequency shift caused by the cross-Kerr effect is
also proportional to the Kittel-mode magnon excitation number, which is
$\Delta_{h,c}=2K_{cross}\langle b^{\dagger}b\rangle$. Here the latter
subscript $c$ represents the cross-Kerr. A similar mechanism occurs when the
drive field is applied solely on the HMS mode, as shown in the bottom panel of
Fig. 1(b). It is worth mentioning that the scale ratios of the frequency
shifts drawn in the figure are not set arbitrarily but roughly correspond to
the experimental results discussed in the following section.

##  III EXPERIMENTAL RESULTS

The transmission spectra of the coupled system are measured versus the magnet
coil current as shown in Fig. 2(a). The large anti-crossing indicates the
strong coupling between the Kittel mode and cavity mode, and the coherent
coupling strength is 40.5 MHz. At the left part of the transmission mapping,
when the magnet coil current is about 4.75 A, we find a small coupling
feature, which is attributed to the HMS mode. The coupling strength is 2 MHz.
The detuning between the Kittel mode and HMS mode is 301 MHz.

###  III.1 Drive the Kittel mode

In Fig. 2(b), we apply the drive field with a power of 25 dBm at the fixed
frequency of 9.8 GHz on the Kittel mode and measure the transmission spectra
versus the magnetic coil current. Typical bistable frequency jumps can be
observed in both the Kittel mode and HMS mode, induced by the self-Kerr effect
and cross-Kerr effect, respectively. This measurement is equivalent to
sweeping the frequency of the drive field while fixing the frequencies of the
magnon modes. We define the drive field detuning as
$\delta_{k}=\omega_{k}-\omega_{d}$. It should be noted that our system works
in the double-dispersive regime, where $\Lambda_{c,k}=\omega_{c}-\omega_{k}\gg
g_{k}$ and $\Lambda_{h,k}=\omega_{h}-\omega_{k}\gg g_{k,h}$. The coherent
couplings induced frequency shifts approximately equal
$g_{k}^{2}/\Lambda_{c,k}$ and $g_{k,h}^{2}/\Lambda_{h,k}$, respectively
Y.P.Wang2016 ; GQ2019 . In the dispersive regime, they are all small terms.
Also, the coupling strength between the HMS mode and cavity mode is relatively
small. In this case, we can intently pay attention to the Kerr nonlinearity
induced effects. Then we extract and plot the frequency shifts of the two
magnon modes versus $\delta_{k}$ in Fig. 2(c). The driven Kittel mode has a
maximum frequency shift of about -60 MHz, and the undriven HMS mode has a
synchronous frequency shift of about -150 MHz.

The frequency shift $\Delta_{k,s}$ of the driven Kittel mode can be fitted by
the following equation Y.P.Wang2016 ; GQ2019 :

$\left[{{{\left({{\Delta_{k,s}}+{\delta_{k}}}\right)}^{2}}+{(\gamma_{k}/2)}^{2}}\right]{\Delta_{k,s}}-{c_{k}}{P_{d}}=0,$
(2)

where ${\gamma_{k}/2\pi}=11.6$ MHz is the Kittel mode linewidth, $c_{k}$ is
the coefficient of drive efficiency on the Kittel mode, and $P_{d}$ is the
drive power. For the undriven HMS mode, the frequency shift can be fitted by
${\Delta_{h,c}}=\frac{K_{cross}}{K_{k,s}}\Delta_{k,s}$. The fitting curves are
shown as solid black lines in Fig. 2(c), which are in good agreement with the
experimental results. We get the ratio of $K_{cross}/K_{s,k}=$2.5, which
indicates that the cross-Kerr nonlinearity is 2.5 times larger than the self-
Kerr nonlinearity of the Kittel mode.

Figure 3: Drive on the HMS mode. (a) The transmission spectra are measured
versus magnet coil current, while a drive field of 10.1 GHz and 25 dBm is
applied on the HMS mode. The HMS mode has a frequency shift $\Delta_{h,s}$ due
to the self-Kerr effect, and the Kittel mode has a smaller frequency shift
$\Delta_{k,c}$ caused by the cross-Kerr effect. (b) The frequency shifts of
the Kittel mode and HMS mode are extracted from (a), and plotted versus
$\delta_{h}=\omega_{h}-\omega_{d}$. The experimental results are fitted by Eq.
(3) and $\Delta_{k,c}=\frac{K_{cross}}{K_{h,s}}\Delta_{h,s}$, which gives the
ratio $K_{cross}/K_{h,s}=0.49$.

###  III.2 Drive the HMS mode

The cross-Kerr effect implies that the HMS mode excitation will also induce a
frequency shift of the Kittel mode. Next, we apply the drive field (25 dBm) at
a fixed frequency of 10.1 GHz on the HMS mode and sweep the magnetic field.
The measured transmission mapping is shown in Fig. 3(a). As expected, the HMS
mode has a frequency shift due to the self-Kerr effect, and simultaneously the
cross-Kerr effect shifts the Kittel mode frequency. We extract the frequency
shifts of these two modes versus the drive field detuning
$\delta_{h}=\omega_{h}-\omega_{d}$ as shown in Fig.3(b). Similarly, the self-
Kerr effect induced frequency shift $\Delta_{h,s}$ of the HMS mode obeys the
following equation:

$\left[{{{\left({{\Delta_{h,s}}+{\delta_{h}}}\right)}^{2}}+{(\gamma_{h}/2)}^{2}}\right]{\Delta_{h,s}}-{c_{h}}{P_{d}}=0,$
(3)

where ${\gamma_{h}/2\pi}=5.0$ MHz is the HMS mode linewidth, $c_{h}$ is the
coefficient of drive efficiency on the HMS mode. The frequency shift
$\Delta_{k,c}$ of the Kittel mode induced by the cross-Kerr effect can be
further fitted by ${\Delta_{k,c}}=\frac{K_{cross}}{K_{h,s}}\Delta_{h,s}$. The
fitting curves are shown as black lines in Fig. 3(b). We get
$K_{cross}/K_{h,s}=0.50$.

From the above two sets of experiments, we get $K_{h,s}/K_{k,s}=5$, which
shows that the HMS mode’s self-Kerr nonlinearity is five times larger than
that of the Kittle mode. In any case, the cross-Kerr effect plays a vital role
in the nonlinear process, and it acts as a ruler telling us the ratio between
the nonlinear coefficients of different spin wave modes. More specifically, by
comparing the amounts of the frequency shifts observed in the two sets of
experiments where the drive powers are both 25 dBm, we can conclude that
pumping the HMS mode is more challenging than pumping the Kittel mode with the
same loop antenna, which may be caused by the small dipole moment of the HMS
mode.

Figure 4: The ratios between the cross-Kerr and the two self-Kerr coefficients
at different drive field frequencies. The ratios are found to be nearly
constant. The mean value of $K_{cross}/K_{k,s}$ is 2.49 and the mean value of
$K_{cross}/K_{h,s}$ is 0.53.

###  III.3 Stable feature under various drive frequencies

Finally, we carry out a series of experiments where the drive field frequency
is changed. We measure the ratios of the cross-Kerr to the two self-Kerr
coefficients. As shown in Fig.4, the ratios keep consistent at different drive
field frequencies, which imply that the observed nonlinearity is the intrinsic
property of this ferrimagnetic crystal. This stable feature provides more
guarantees for the information conversion and storage among different modes in
a single YIG sphere.

##  IV CONCLUSION

The cross-Kerr effect is an important nonlinearity, which describes the effect
of excitations in one mode on another mode’s resonant frequency. In this work,
we investigate the cross-Kerr effect in a ferrimagnetic spin ensemble, in
which one Kittel mode and one higher-order magnetostatic mode co-exist. The
cross-Kerr interaction between them plus their individual self-Kerr effects
give rise to fascinating phenomena. In the experiment, we drive the Kittel
model and the higher-order magnetostatic mode, respectively, and obtain the
ratios between the cross-Kerr coefficient and their self-Kerr coefficients.
Our experiment opens up a new path to study nonlinear effects in magnetic
materials and spin ensembles. The mutual interaction between different
magnetostatic modes in a single YIG sample can also provide new degrees of
freedom for cavity spintronics and cavity magnonics Roadmap . We anticipate
that this finding may stimulate more designs and applications of cavity
magnonics.

###### Acknowledgements.

The authors thank Y. H. Li for helpful discussions. This work is supported by
the National Natural Science Foundation of China (Grants No. 11934010, No.
U1801661, and No. 12174329), Zhejiang Province Program for Science and
Technology (Grant No. 2020C01019), the Fundamental Research Funds for the
Central Universities (No. 2021FZZX001-02), and the China Postdoctoral Science
Foundation (No. 2019M660137).

Figure 5: (a) Original transmission spectrum mapping. (b) The transmission
spectrum $|S_{21}|^{2}$ correspond to the black dashed line labelled in (a) at
a certain bias magnetic field. (c) The Kerr effects induced frequency shifts
are extracted. We subtract the linear frequency shift caused by the external
bias field from all data points. The horizontal axis is changed to
${\delta_{k}}={\omega_{k}}-{\omega_{d}}$ for fitting with Eq. (2).

##  APPENDIX A INTERACTION BETWEEN THE KITTEL MODE AND HMS MODE

First, we derive the Hamiltonian of the magnon modes without including the
cross-Kerr effect. Under the bias magnetic field $B_{0}$, the Hamiltonian of
the Kittel mode is written as (HMS mode has a similar form) Y.P.Wang2016

${H_{k}}=-\int\limits_{{V}}{{\bf{M}}_{K}}\cdot{{\bf{B}}_{0}}d\tau-\frac{{{\mu_{0}}}}{2}\int\limits_{{V}}{{\bf{M}}_{K}\cdot{{\bf{H}}_{\rm{an}}}d\tau},$
(4)

where the first term represents the Zeeman energy and the second term is the
magnetocrystalline anisotropy energy.
${\bf{M}}_{K}=\left({{M_{K,x}},{M_{K,y}},{M_{K,z}}}\right)$ is the sub-
magnetization corresponding to the Kittel mode. $V$ is the volume of the YIG
sphere, $\mu_{0}$ is the vacuum permeability, and $\mathbf{H}_{\rm{an}}$ is
the anisotropic field due to the magnetocrystalline anisotropy in the YIG
crystal.

We adopt the direction of the bias magnetic $\mathbf{B}_{0}$ as the $z$
direction ($\mathbf{B}_{0}=B_{0}\mathbf{e}_{z}$). When the [110] crystal axis
of the YIG sphere is aligned along the bias magnetic field, the anisotropic
field is given by stancil2009

$\mathbf{H}_{\rm{an}}=\frac{3K_{\rm{an}}M_{x}}{\mu_{0}M^{2}}\mathbf{e}_{x}+\frac{9K_{\rm{an}}M_{y}}{4\mu_{0}M^{2}}\mathbf{e}_{y}+\frac{K_{\rm{an}}M_{z}}{\mu_{0}M^{2}}\mathbf{e}_{z},$
(5)

where $K_{\rm{an}}$ is the first-order magnetocrystalline anisotropy constant.
Then, the Hamiltonian of the Kittel mode has the form

$\displaystyle H_{k}$ $\displaystyle=$ $\displaystyle-
B_{0}M_{K,z}V-\frac{3K_{\rm{an}}M_{K,x}}{2M_{K}^{2}}\mathbf{e}_{x}$ (6)
$\displaystyle-$
$\displaystyle\frac{9K_{\rm{an}}M_{K,y}}{8M_{K}^{2}}\mathbf{e}_{y}-\frac{K_{\rm{an}}M_{K,z}}{2M_{K}^{2}}\mathbf{e}_{z}.$

By using the macrospin operator
$\mathbf{S}=\mathbf{M}V/\hbar\gamma\equiv(S_{x},S_{y},S_{z})$, where $\gamma$
is the electronic gyromagnetic ratio Soykal2010 , the Hamiltonian $H_{k}$ can
be written as

$\displaystyle H_{k}/\hbar$ $\displaystyle=$ $\displaystyle-\gamma
B_{0}S_{K,z}-\frac{3\hbar K_{\rm{an}}\gamma^{2}}{2M_{K}^{2}V}S_{K,x}^{2}$ (7)
$\displaystyle-$ $\displaystyle\frac{9\hbar
K_{\rm{an}}\gamma^{2}}{8M_{K}^{2}V}S_{K,y}^{2}-\frac{\hbar
K_{\rm{an}}\gamma^{2}}{2M_{K}^{2}V}S_{K,z}.$

We harness the Holstein-Primakoff transformation Holstein1940 :
$S_{K}^{+}=\sqrt{2S_{K}-b^{\dagger}b}b$,
$S_{K}^{-}=b^{\dagger}\sqrt{2S_{K}-b^{\dagger}b}$, and
$S_{K,z}=S_{K}-b^{{\dagger}}b$, where $S_{K}$ is the total spin number, and
$S_{K}^{\pm}=S_{K,x}\pm iS_{K,y}$. The excited magnon number is small compared
with the total spin number $S_{K}$, and the higher-order terms can be
neglected in the Holstein-Primakoff transformation. In this case, we have
$S^{+}\approx\sqrt{2S}b$ and $S^{-}\approx b^{{\dagger}}\sqrt{2S}$. Under the
rotating-wave approximation, we finally obtain the Hamiltonian

${H_{k}}/\hbar={\omega_{k^{\prime}}}{b^{\dagger}}b+{K_{k,s}}{b^{\dagger}}b{b^{\dagger}}b,$
(8)

where

${K_{k,s}}=13\hbar K_{\rm{an}}\gamma^{2}/(16M_{K}^{2}V),$

and the value of $K_{\rm{an}}$ can be found in Ref. stancil2009 .

For the interaction Hamiltonian between the Kittel mode and HMS mode, we
consider these two modes as two sets of submagnetization, and the interaction
between them can be written as

$\begin{split}H_{k,h}&=\beta\int\limits_{V}{\textbf{M}_{K}\cdot{\textbf{M}}_{H}}d\tau\\\
&=\frac{\beta\hbar^{2}\gamma^{2}}{V}\textbf{S}_{K}\cdot\textbf{S}_{H}\\\
&=\frac{\beta\hbar^{2}\gamma^{2}}{V}(S_{K,x}S_{H,x}+S_{K,y}S_{H,y}+S_{K,z}S_{H,z}),\end{split}$
(9)

where $\beta$ is a coefficient that measures the mode overlapping between the
Kittel mode and HMS mode. This mode overlap originates from the partial local
spins shared by the Kittel mode and other spin wave modes A.Gloppe2019 ;
stancil2009 . Under the Holstein-Primakoff transformation and rotating-wave
approximation, we can get

$\begin{split}{H_{k,h}}&=\frac{\beta\hbar^{2}{\gamma^{2}}}{V}[-S_{H}{b^{\dagger}}b-S_{K}{c^{\dagger}}c\\\
&+{b^{\dagger}}b{c^{\dagger}}c+\sqrt{{S_{K}}{S_{H}}}({b^{\dagger}}c+b{c^{\dagger}})].\\\
\end{split}$ (10)

The first and second terms produce constant frequency shifts in the Kittel
mode and HMS mode, respectively. The third term refers to the coherent
coupling between the Kittel mode and HMS mode. In our experiment, the detuning
between the Kittel mode and HMS mode is sufficiently large, and the system is
in the dispersive regime. If we tune the Kittel mode close to the HMS mode, we
can predict level repulsion. A related effect has been observed in a previous
study Paul .

The total Hamiltonian of the cavity magnonic system becomes

$\begin{split}H/\hbar&=\omega_{c}a^{\dagger}a+\omega_{k}b^{\dagger}b+\omega_{h}c^{\dagger}c\\\
&+K_{k,s}b^{\dagger}bb^{\dagger}b+K_{h,s}c^{\dagger}cc^{\dagger}c+K_{cross}b^{\dagger}bc^{\dagger}c\\\
&+g_{k}\left(a^{\dagger}b+b^{\dagger}a\right)+g_{h}\left(a^{\dagger}c+c^{\dagger}a\right)+g_{k,h}\left(b^{\dagger}c+c^{\dagger}b\right)\\\
&+\Omega_{k}(b^{\dagger}e^{-i{\omega_{d}}t}+be^{i{\omega_{d}}t})\\\
&+\Omega_{h}(c^{\dagger}e^{-i{\omega_{d}}t}+ce^{i{\omega_{d}}t}),\end{split}$
(11)

where
${\omega_{k}}={\omega_{k^{\prime}}}-{\beta\hbar{\gamma^{2}}{S_{H}}}/{V},~{}{\omega_{h}}={\omega_{h^{\prime}}}-{\beta\hbar{\gamma^{2}}{S_{K}}}/{V},~{}{K_{cross}}={\beta\hbar{\gamma^{2}}}/{V},$
and $g_{k,h}=\beta\hbar{\gamma^{2}}\sqrt{{S_{K}}{S_{H}}}/{V}$.

##  APPENDIX B DATA PROCESSING METHOD

The transmission spectrum mapping consists of a sequence of transmission
spectra measured at different bias magnetic fields, as shown in Fig. 5(a). The
individual transmission spectrum is depicted in Fig.5(b). The anti-resonances
correspond to various modes; for example, the green triangle dot represents
the HMS mode, and the circle red dot denotes the Kittel mode. We extract the
data points from the sequence of transmission spectra and plot them in the
left panel of Fig.5(c). To obtain a plot of frequency shift versus driving
detuning that is easy to fit with our theoretical model, we subtract the
linear component that increases with the external bias field from all data
points. Then, the remaining component is the contribution of frequency shift
caused by the driving field and Kerr effects. The middle panel of Fig. 5(c)
depicts the linear frequency shift caused by the increased bias magnetic
field. In order to fit the experimental data with Eq. (2), we replace the
horizontal axis with drive detuning ${\delta_{k}}={\omega_{k}}-{\omega_{d}}$.

## References

  * (1) R. W. Boyd, Nonlinear optics (Academic press, London, 2020).
  * (2) A. Auere, H. Jussila, Y. Dai, Y. Wang, H. Lipsanen, and Z. Sun, Nonlinear optics with 2d layered materials, Advanced Materials 30, 1705963(2018).
  * (3) X. Zhang, Q.-T. Cao, Z. Wang, Y.-x. Liu, C.-W. Qiu, L. Yang, Q. Gong, and Y.-F. Xiao, Symmetry-breaking-induced nonlinear optics at a microcavity surface, Nature Photonics 13, 21 (2019).
  * (4) M. Kauranen and A. V. Zayats, Nonlinear plasmonics, Nature Photonics 6, 737 (2012)
  * (5) M. Mesch, B. Metzger, M. Hentschel, and H.Giessen, Nonlinear Plasmonic Sensing, Nano Letters 16, 3155 (2016)
  * (6) J.-H. Zhong, J. Vogelsang, J.-M. Yi, D. Wang, L. Wittenbecher, S. Mikaelsson, A. Korte, A. Chimeh, C. L. Arnold, P. Schaaf, E. Runge, Anne L’ Huillier, A. Mikkelsen, and C. Lienau, Nonlinear plasmon-exciton coupling enhances sum-frequency generation from a hybrid metal/semiconductor nanostructure, Nature Communications 11, 1464 (2020).
  * (7) M. F. Hamilton, D. T. Blackstock, Nonlinear acoustics, (Academic press, San Diego, 1998).
  * (8) X. Fang, J. Wen, B. Bonello, J. Yin, and D. Yu, Ultra-low and ultra-broad-band nonlinear acoustic metamaterials, Nature Communications 8, 1288 (2017).
  * (9) A. Darabi, L. Fang, A. Mojahed, M. D. Fronk, A. F. Vakakis, and M. J. Leamu, Broadband passive nonlinear acoustic diode, Phys. Rev. B 99, 214305 (2019).
  * (10) A. H. Nayfeh and D. T. Mook, Nonlinear oscillations (John Wiley & Sons, New York, 2008).
  * (11) Y. V. Kartashov, G. E. Astrakharchik, B. A. Malomed, and L. Torner, Frontiers in multidimensional self-trapping of nonlinear fields and matter, Nature Reviews Physics 1, 185 (2019).
  * (12) Y. Zuo, W. Yu, C. Liu, X. Cheng, R. Qiao, J.liang, X. Zhou, J. Wang, M. Wu, Y. Zhao, P. Gao, S. Wu, Z. Sun, K. Liu, X. Bai, and Z. Liu, Optical fibres with embedded two-dimensional materials for ultrahigh nonlinearity, Nature Nanotechnology 15, 987 (2020).
  * (13) S. Ding, G. Maslennikov, R. Habl’́utzel, and D. Matsukevich, Cross-Kerr Nonlinearity for Phonon Counting, Phys. Rev. Lett. 119,193602 (2017).
  * (14) B. He, A. V. Sharypov, J. Sheng, C. Simon, and M. Xiao, Two-photon dynamics in coherent Rydberg atomic ensemble, Phys. Rev. Lett. 112, 133606 (2014).
  * (15) K. Xia, F. Nori, and M. Xiao, Cavity-free optical isolators and circulators using a chiral cross-Kerr nonlinearity, Phys. Rev. Lett. 121, 203602 (2018).
  * (16) J. Sinclair, D. Angulo, N. Lupu-Gladstein, K. Bonsma-Fisher, and A. M. Steinberg, Observation of a large, resonant, cross-Kerr nonlinearity in a cold Rydberg gas, Phys. Rev. Research , 033193 (2019).
  * (17) I.-C. Hoi, A. F. Kockum, T. Palomaki, T. M. Stace, B. Fan, L. Tornberg, S. R. Sathyamoorthy, G. Johansson, P. Delsing, and C. M. Wilson, Giant Cross–Kerr Effect for Propagating Microwaves Induced by an Artificial Atom,Phys. Rev. Lett. 111, 053601 (2013).
  * (18) M. Kounalakis, C. Dickel, A. Bruno, N. K. Langford and G. A. Steele, Tuneable hopping and nonlinear cross-Kerr interactions in a high-coherence superconducting circuit, npj Quantum Information 4, 38 (2018).
  * (19) A. Vrajitoarea, Z. Huang, P. Groszkowski, J. Koch, and A. A. Houck, Quantum control of an oscillator using a stimulated Josephson nonlinearity, Nature Physics 16, 211 (2020)
  * (20) Q. A. Turchette, C. J. Hood, W. Lange, H. Mabuchi, and H. J. Kimble, Measurement of Conditional Phase Shifts for Quantum Logic, Phys. Rev. Lett. 75, 4710 (1995).
  * (21) D. J. Brod and J. Combes, Passive CPHASE Gate via Cross-Kerr Nonlinearities, Phys. Rev. Lett. 117, 080502 (2016).
  * (22) M. D. Lukin and A. Imamoǧlu, Nonlinear Optics and Quantum Entanglement of Ultraslow Single Photons, Phys. Rev. Lett. 84, 1419 (2000).
  * (23) Y. B. Sheng, F. G. Deng, and H. Y. Zhou, Efficient polarization-entanglement purification based on parametric down-conversion sources with cross-Kerr nonlinearity, Phys. Rev. A 77, 042308 (2008).
  * (24) X. Wang, A. Miranowicz, and F. Nori, Ideal Quantum Nondemolition Readout of a Flux Qubit without Purcell Limitations, Phys. Rev. Applied 12, 064037 (2019).
  * (25) R. Dassonneville, T. Ramos, V. Milchakov, L. Planat, É. Dumur, F. Foroughi, J. Puertas, S. Leger, K. Bharadwaj, J. Delaforce, C. Naud, W. Hasch-Guichard, J. J. García-Ripoll, N. Roch, and O. Buisson, Fast High-Fidelity Quantum Nondemolition Qubit Readout via a Nonperturbative Cross-Kerr Coupling, Phys. Rev. X 10, 011045 (2020).
  * (26) H. Huebl, C. W. Zollitsch, J. Lotze, F. Hocke, M. Greifenstein, A. Marx, R. Gross, and S. T. B. Goennenwein, High Cooperativity in Coupled Microwave Resonator Ferrimagnetic Insulator Hybrids, Phys. Rev. Lett. 111, 127003 (2013).
  * (27) Y. Tabuchi, S. Ishino, T. Ishikawa, R. Yamazaki, K. Usami, and Y. Nakamura, Hybridizing Ferromagnetic Magnons and Microwave Photons in the Quantum Limit, Phys. Rev. Lett. 113, 083603 (2014).
  * (28) X. Zhang, C.-L. Zou, L. Jiang, and H. X. Tang, Strongly Coupled Magnons and Cavity Microwave Photons, Phys. Rev. Lett. 113, 156401 (2014).
  * (29) M. Goryachev, W. G. Farr, D. L. Creedon, Y. Fan, M. Kostylev, and M. E. Tobar, High-Cooperativity Cavity QED with Magnons at Microwave Frequencies, Phys. Rev. Applied 2, 054002 (2014).
  * (30) L. Bai, M. Harder, Y. P. Chen, X. Fan, J. Q. Xiao, and C.-M. Hu, Spin pumping in electrodynamically coupled magnon-photon systems, Phys. Rev. Lett. 114, 227201 (2015).
  * (31) Y. Cao, P. Yan, H. Huebl, S. T. B. Goennenwein, and G. E. W. Bauer, Exchange magnon-polaritons in microwave cavities, Phys. Rev. B 91, 094423 (2015).
  * (32) X. Zhang, C.-L. Zou, N. Zhu, F. Marquardt, L. Jiang, and H. X. Tang, Magnon dark modes and gradient memory, Nature Communications 6, 8914 (2015)
  * (33) Y. Tabuchi, S. Ishino, A. Noguchi, T. Ishikawa, R. Yamazaki, K. Usami, and Y. Nakamura, Coherent coupling between a ferromagnetic magnon and a superconducting qubit, Science 349, 405 (2015).
  * (34) D. Zhang, X. M. Wang, T. F. Li, X. Q. Luo, W. Wu, F. Nori, and J. Q. You, Cavity quantum electrodynamics with ferromagnetic magnons in a small yttrium-iron-garnet sphere, npj Quantum Information 1, 15014 (2015).
  * (35) J. A. Haigh, A. Nunnenkamp, A. J. Ramsay, and A. J. Ferguson, Triple-Resonant Brillouin Light Scattering in Magneto-Optical Cavities, Phys. Rev. Lett. 117, 133602 (2016).
  * (36) H. Maier-Flaig, M. Harder, R. Gross, H. Huebl, and S. T. B. Goennenwein, Spin pumping in strongly coupled magnon-photon systems, Phys. Rev. B 94, 054433 (2016).
  * (37) A. Osada, R. Hisatomi, A. Noguchi, Y. Tabuchi, R. Yamazaki, K. Usami, M. Sadgrove, R. Yalla, M. Nomura, and Y. Nakamura, Cavity Optomagnonics with Spin-Orbit Coupled Photons, Phys. Rev. Lett. 116, 223601 (2016).
  * (38) X. Zhang, N. Zhu, C.-L. Zou, and H. X. Tang, Optomagnonic Whispering Gallery Microresonators, Phys. Rev. Lett. 117, 123605 (2018).
  * (39) D. Zhang, X.-Q. Luo, Y.-P. Wang, T.-F. Li and J. Q. You, Observation of the exceptional point in cavity magnon-polaritons, Nature Communications 8, 1368 (2017).
  * (40) Y.-P. Wang, G.-Q. Zhang, D. Zhang, T.-F. Li, C.-M. Hu, and J. Q. You, Bistability of Cavity Magnon Polaritons, Phys. Rev. Lett. 120, 057202 (2018).
  * (41) J. Li, S.-Y. Zhu, and G. S. Agarwal, Magnon-Photon-Phonon Entanglement in Cavity Magnomechanics, Phys. Rev. Lett. 121, 203601 (2018).
  * (42) A. Osada, A. Gloppe, R. Hisatomi, A. Noguchi, R. Yamazaki, M. Nomura, Y. Nakamura, and K. Usami, Brillouin Light Scattering by Magnetic Quasivortices in Cavity Optomagnonics, Phys. Rev. Lett. 120, 133602 (2018).
  * (43) M. Harder, Y. Yang, B. M. Yao, C. H. Yu, J. W. Rao, Y. S. Gui, R. L. Stamps, and C.-M. Hu, Level Attraction Due to Dissipative Magnon-Photon Coupling, Phys. Rev. Lett. 121, 137203 (2018).
  * (44) Y.-P. Wang, J. W. Rao, Y. Yang, P.-C. Xu, Y. S. Gui, B. M. Yao, J. Q. You, and C.-M. Hu, Nonreciprocity and Unidirectional Invisibility in Cavity Magnonics, Phys. Rev. Lett. 123, 127202 (2019).
  * (45) D. Lachance-Quirion, S. P. Wolski, Y. Tabuchi, S. Kono, K. Usami, and Y. Nakamura, Entanglement-based single-shot detection of a single magnon with a superconducting qubit, Science 367, 425 (2020).
  * (46) Y. Yang, Y.-P. Wang, J. W. Rao, Y. S. Gui, B. M. Yao, W. Lu, and C.-M. Hu, Unconventional Singularity in Anti-Parity-Time Symmetric Cavity Magnonics, Phys. Rev. Lett.125, 147202 (2020).
  * (47) J. Xu, C. Zhong, X. Han, D. Jin, L. Jiang, and X. Zhang, Floquet Cavity Electromagnonics, Phys. Rev. Lett. 125, 237201 (2020).
  * (48) S. P. Wolski, D. Lachance-Quirion, Y. Tabuchi, S. Kono, A. Noguchi, K. Usami, and Y. Nakamura, Dissipation-Based Quantum Sensing of Magnons with a Superconducting Qubit, Phys. Rev. Lett. 125, 117701 (2020).
  * (49) J. W. Rao, P. C. Xu, Y. S. Gui, Y. P. Wang, Y. Yang, Bimu Yao, J. Dietrich, G. E. Bridges, X. L. Fan, D. S. Xue and C.-M. Hu, Interferometric control of magnon-induced nearly perfect absorption in cavity magnonics, Nature Communications 12, 1933 (2021).
  * (50) J. Xu, C. Zhong, X. Han, D. Jin, L. Jiang, and X. Zhang, Coherent Gate Operations in Hybrid Magnonics, Phys. Rev. Lett. 126, 207202 (2021).
  * (51) C. A. Potts, E. Varga, V. A. S. V. Bittencourt, S. V. Kusminskiy, and J. P. Davis, Dynamical Backaction Magnomechanics, Phys. Rev. X 11, 031053 (2021).
  * (52) D. Lachance-Quirion, Y. Tabuchi, A. Gloppe, K. Usami, and Y. Nakamura, Hybrid quantum systems based on magnonics, Applied Physics Express 12, 070101 (2019).
  * (53) J. Li, Y.-P. Wang, W.-J. Wu, S.-Y. Zhu, and J.Q. You, Quantum Network with Magnonic and Mechanical Nodes, PRX Quantum 2, 040344 (2021).
  * (54) J. F. Dillon, Jr. Magnetostatic Modes in Ferrimagnetic Spheres, Phys. Rev. 112, 59 (1958).
  * (55) P. C. Fletcher and R. O. Bell, Ferrimagnetic Resonance Modes in Spheres, Journal of Applied Physics 30, 687 (1959).
  * (56) J. Graf, H. Pfeifer, F. Marquardt, and S. V. Kusminskiy, Cavity optomagnonics with magnetic textures: Coupling a magnetic vortex to light, Phys. Rev. B 98, 241406 (2018).
  * (57) J. A. Haigh, N. J. Lambert, S. Sharma, Y. M. Blanter, G. E. W. Bauer, and A. J. Ramsay, Selection rules for cavity-enhanced Brillouin light scattering from magnetostatic modes, Phys. Rev. B 97, 214423 (2018)
  * (58) A. Osada, A. Gloppe, Y. Nakamura, and K. Usami, Orbital angular momentum conservation in Brillouin light scattering within a ferromagnetic sphere, New Journal of Physics, 20, 103018 (2018).
  * (59) A. Gloppe, R. Hisatomi, Y. Nakata, Y. Nakamura, and K. Usami, Resonant Magnetic Induction Tomography of a Magnetized Sphere, Phys. Rev. Applied 12, 014061 (2019).
  * (60) 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, Magnon Kerr effect in a strongly coupled cavity-magnon system, Phys. Rev. B 94, 224410 (2016).
  * (61) G.-Q. Zhang, Y.-P. Wang and J. Q. You, Theory of the magnon Kerr effect in cavity magnonics, Science China Physics, Mechanics & Astronomy 62, 987511 (2019).
  * (62) C. Kong, H. Xiong, and Y. Wu, Magnon-Induced Nonreciprocity Based on the Magnon Kerr Effect, Phys. Rev. Applied 12, 034001 (2019).
  * (63) M. X. Bi, X. H. Yan, Y. Zhang, and Y. Xiao, Tristability of cavity magnon polaritons, Phys. Rev. B 103, 104411 (2021).
  * (64) R.-C. Shen, Y.-P. Wang, J. Li, S.-Y. Zhu, G. S. Agarwal, and J. Q. You, Long-Time Memory and Ternary Logic Gate Using a Multistable Cavity Magnonic System, Phys. Rev. Lett. 127, 183202 (2021).
  * (65) A. G. Gurevich and G. A Melkov, Magnetization oscillations and waves (CRC press, Boca Raton, 2020).
  * (66) Y. Li, W. Zhang, V. Tyberkevych, W.-K Kwok, A Hoffmann, and V. Novosad, Hybrid magnonics: Physics, circuits, and applications for coherent information processing, Journal of Applied Physics 128, 130902 (2020).
  * (67) B. Z. Rameshti, S. V. Kusminskiy, J. A. Haigh, K. Usami, D. Lachance-Quirion, Y. Nakamura, C.-M. Hu, H. X. Tang, G. E. W. Bauer, Y. M. Blanter, Cavity Magnonics, arXiv:2106.09312
  * (68) H. Y. Yuan, Y. Cao, A. Kamra, R. A. Duine, and P. Yan, Quantum magnonics: when magnon spintronics meets quantum information science, arXiv:2111.14241
  * (69) C. Zhang, C. Jia, Y. Shi, C. Jiang, D. Xue, C. K. Ong, and G. Chai, Nonreciprocal multimode and indirect couplings in cavity magnonics, Phys. Rev. B 103, 184427 (2021).
  * (70) A. V. Chumak, et al. Roadmap on Spin-Wave Computing, arXiv:2111.00365
  * (71) D. D. Stancil and A. Prabhakar, _Spin Waves: Theory and Applications_ (Springer, New York, 2009).
  * (72) Ö. O. Soykal and M. E. Flatté, Strong Field Interactions between a Nanomagnet and a Photonic Cavity, Phys. Rev. Lett. 104, 077202 (2010).
  * (73) T. Holstein and H. Primakoff, Field Dependence of the Intrinsic Domain Magnetization of a Ferromagnet, Phys. Rev. 58, 1098 (1940).
  * (74) P. C. Fletcher and I. H. Solt Jr, Journal of Applied Physics 30, S181 (1959).


* Bicer  _et al._ (2014) M. Bicer _et al._ (TLEP Design Study Working Group),  _Proceedings, 2013 Community Summer Study on the Future of U.S. Particle Physics: Snowmass on the Mississippi (CSS2013): Minneapolis, MN, USA, July 29-August 6, 2013_ , JHEP 01, 164 (2014), arXiv:1308.6176 [hep-ex] .
  * Liu  _et al._ (2017) Z. Liu, L.-T. Wang, and H. Zhang, Chin. Phys. C41, 063102 (2017), arXiv:1612.09284 [hep-ph] .
  * Curtin  _et al._ (2018) D. Curtin _et al._ ,  _Long-Lived Particles at the Energy Frontier: The MATHUSLA Physics Case_ , Tech. Rep. (2018) arXiv:1806.07396 [hep-ph] .
  * Shtabovenko  _et al._ (2016) V. Shtabovenko, R. Mertig, and F. Orellana, Comput. Phys. Commun. 207, 432 (2016), arXiv:1601.01167 [hep-ph] .
  * Sirunyan  _et al._ (2018b) A. M. Sirunyan _et al._ (CMS), Phys. Rev. Lett. 120, 071802 (2018b), arXiv:1709.05543 [hep-ex] .
  * Vernieri (2017) C. Vernieri (CMS),  _Proceedings, 2017 European Physical Society Conference on High Energy Physics (EPS-HEP 2017): Venice, Italy, July 5-12, 2017_ , PoS EPS-HEP2017, 349 (2017), arXiv:1711.10508 [hep-ex] .
  * Collaboration (2018) C. Collaboration (CMS),  _Search for low mass resonances decaying into bottom quark-antiquark pairs in pp collisions at $\sqrt{s}=13~{}\mathrm{TeV}$_, Tech. Rep. (2018).


# How about that Bayes: Bayesian techniques and the simple pendulum

Matthew Heffernan<EMAIL_ADDRESS>Department of Physics, McGill
University, 3600 University Street, Montreal, QC, H3A 2T8, Canada

###### Abstract

Physics increasingly relies on Bayesian techniques for rigorous data analysis
and model-to-data comparison. This paper describes how these methods can be
implemented to answer questions of relevance to the teaching of introductory
physics concepts. It demonstrates the proper application of statistical
approaches to quantitatively assess when more complex calculations are
necessitated by data. The case study of a simple pendulum allows for the
incorporation of the precision commonly found in introductory physics
laboratories in order to rigorously assess model preference between the small
angle and exact expressions for the period. This extends the current simple
pendulum literature’s focus on individual assessments of different
approximations to the exact solution to place such discussions on a data-
driven foundation and provides guidance to instructors for teaching laboratory
design.

##  I Introduction

Introductory physics textbooks teach that the period of a simple pendulum
displaced by a small angle may be found by applying the approximation
$\sin{\theta}\approx\theta$, which holds when $\theta$ is small. However, the
period of a pendulum also depends on the initial angular displacement.
Attempts to increase accuracy in approximations of the pendulum period have so
far focused on more closely approximating the full non-analytic integral
expression for the period so students in undergraduate laboratories can
investigate the dependence of the period on initial angular displacement Lima
and Arun (2006). Only recently has the idea of model comparison and complexity
been applied to the pendulum problem Blais (2020).

Guidance for the use of the small angle approximation in simple pendula has
not previously been rigorously quantified using realistic uncertainties
measured in teaching laboratories. Often, the only suggestion is that after
$\sim 15^{\circ}$, the difference between $\sin{\theta}$ and $\theta$ in
radians exceeds 1% and the small angle approximation should no longer be used
Giancoli (2016). This does not take into account the variety of measurement
uncertainties detailed in introductory physics laboratory materials or if more
elaborate formulae are demanded by the data. This study demonstrates that
students with both an exact and approximate analytical formula for the period
of a pendulum could successfully use Bayesian tools to establish a preference
between the two at moderate displacements. This can also be extended to the
case of choosing between various more complex approximations. This study
remedies a gap in the literature by quantitatively motivating more complex
formula in teaching laboratories and will provide guidance to instructors and
their students in restricting the initial angular displacement or employing a
more complex formula (c.f. Hinrichsen (2020) for a recent review).

This paper will use the case of the period of a simple pendulum to demonstrate
the techniques of Bayesian data analysis to differentiate the small angle
approximation and a numerical calculation of the exact integral equation.
First, an introduction to Bayesian analysis is provided and the exact
expression for the period of a simple pendulum is derived as well as a common
simple approximation. It is then established that it is possible to detect the
impact of the approximation using realistic measurement uncertainties. Given
this, prior distributions for any relevant parameters are discussed. Prior
predictive checks are made, followed by inference to fit the models to data.
The posterior distribution is then analyzed and posterior predictive checks
are made and are found to be consistent with expectations. Finally, it is
determined at what initial angular displacement weak, moderate, and strong
preference for the exact calculation are found given realistic, but precise,
measurements. This preference can be used to inform recommendations of maximum
initial angular displacement in laboratory materials and textbooks.

##  II Bayesian Analysis

###  II.1 Introduction

Bayesian inference is an increasingly-commonly used tool in physics research,
e.g. gravitational-wave astronomy Smith  _et al._ (2020) and recently in
ultra-relativistic heavy ion collisions Everett  _et al._ (2020), among many
others von Toussaint (2011). However, the techniques and process can be
readily applied to answer pedagogical questions in a principled manner. In
Bayesian statistics, probability is interpreted as a degree of belief in an
event. Bayes’ theorem,

$P(H|d,I)=\frac{P(d|H,I)\hskip 7.22743ptP(H|I)}{P(d|I)},$ (1)

may be derived from the product rule of probability and may be used to solve
statistical inverse problems and infer the probability distribution of model
parameters. Here, $P(H|I)$ is called the “prior” and represents the prior
belief in the hypothesis $H$. $P(d|H,I)$ is the “likelihood” and is the
likelihood for the data $d$ to be true given the hypothesis $H$. $P(d|I)$ is
called the “Bayes evidence” and is typically treated as a normalization
constant such that $\int P(H|d,I)dH=1$. It quantifies how likely one believes
the data to be given other information $I$, which can include the model in
question in which a particular choice of parameters forms a hypothesis $H$.
$P(H|d,I)$ is called the “posterior” and quantifies the belief in any
hypothesis $H$ given data $d$. $I$ denotes any other additional information
such as the choice of physics model or any other reasonable expectations that
are not informed by $d$. These concepts will be demonstrated quantitatively
using the instructive case of the simple pendulum, outlined further in the
next subsection.

Bayes’ Theorem can be numerically evaluated using computational tools such as
Markov Chain Monte Carlo in order to compute the posterior for hypotheses such
as parameter values. This then yields an inferential estimate of the
parameters of the model as well as their structure, which in turn allows for
the precise quantification of the degree of belief. An accessible introduction
to Bayesian inference can be found in Trotta (2008). A guide to Bayesian
modeling in a rigorous workflow may be found in Gelman  _et al._ (2020);
Betancourt (2019).

Bayes’ Theorem and related techniques will be used to ask when, given
realistic data, it is possible to detect weak, moderate, and strong preference
for the exact expression for the period of a pendulum (Eq. 7) relative to the
expression for the period of a pendulum in the small angle approximation (Eq.
4). To do so, realistically-generated pseudodata and principled prior
expectations will be used in order to evaluate Bayes’ Theorem (Eq. 1). This
study employs the use of the Gaussian likelihood function, usually represented
in its logarithmic form, which assumes that uncertainties are normally
distributed around a central value.

$\ln
P(d|H,I)=-\frac{1}{2}\sum_{i=1}^{M}\left[\ln(2\pi\sigma^{2})+\frac{(d_{i}-Y_{\rm
Model})^{2}}{\sigma^{2}}\right],$ (2)

where $Y_{\rm Model}$ is the model calculation of the period and
$\sigma^{2}=\sigma_{\rm data}^{2}+\sigma_{\rm model}^{2}$. In this study, only
one parameter – the gravitational acceleration $g$ – will be varied and thus a
particular value of $g$ constitutes a hypothesis $H$ given the choice of model
or other assumptions $I$. Generalizing this inferential framework to higher
dimensions is trivial.

###  II.2 Models: The period of a simple pendulum

For completeness, the derivation of both the small angle approximation and the
exact expression for the formula of a simple pendulum are included. An ideal
simple pendulum has a bob of mass $m$ suspended from a frictionless attachment
point by a massless rod of length $L$ and moves in a single angular dimension
$\theta$. In the absence of nonconservative forces such as friction and drag,
the equation of motion is

$\ddot{\theta}=-\frac{g\sin(\theta)}{L}$ (3)

where $g$ is the gravitational acceleration. In this study, a fixed value
$L=1\text{m}$ is used as representative of home-made pendula. For small
displacements, $\sin{\theta}\approx\theta$ and Eq. 3 can be identified as a
simple harmonic oscillator. The period of a simple harmonic oscillator is
$T=2\pi\sqrt{\frac{m}{k}}$ and it is trivial to identify $k=\frac{mg}{L}$.
Thus, in this regime,

$T_{0}=2\pi\sqrt{\frac{L}{g}}.$ (4)

An exact expression for the period of a simple pendulum is readily derived
using the conservation of energy. The expression is simplified by assuming the
pendulum to be initially stationary and at its initial angular displacement
$\theta_{0}$. The zero of potential energy is chosen as when $\theta(t)=0$.
The energy conservation equation is thus

$\displaystyle mgL\left(1-\cos\theta_{0}\right)$ $\displaystyle=$
$\displaystyle\frac{1}{2}mL^{2}\dot{\theta}^{2}+mgL(1-\cos\theta).$ (5)

It is then straightforward to solve for $\dot{\theta}$ and integrate $\theta$
from 0 to $\theta_{0}$, corresponding to one-quarter of the period. This
yields

$T=2\sqrt{2}\sqrt{\frac{L}{g}}\int_{0}^{\theta_{0}}\frac{1}{\sqrt{\cos\theta-\cos\theta_{0}}}d\theta.$
(6)

This integral is improper when $\theta=\theta_{0}$, but a simple substitution
$\cos\theta=1-2\sin^{2}(\theta/2)$ may be combined with a change of variables
$\sin\phi=\frac{\sin(\theta/2)}{\sin(\theta_{0}/2)}$ to yield

$T=4\sqrt{\frac{L}{g}}\int_{0}^{\pi/2}\frac{1}{\sqrt{1-\sin^{2}(\theta_{0}/2)\sin^{2}(\phi/2)}}d\phi$
(7)

Lima and Arun (2006); Schery (1976) and may be evaluated numerically using
tools such as SciPy Virtanen  _et al._ (2020).

The small angle approximation (Eq. 4) and the exact expression (Eq. 7) are the
two models considered in this work.

###  II.3 Generating realistic measurements

Realistic measurements are generated by making a set of simple choices
informed by experience in undergraduate teaching laboratories. In these
laboratories, standard advice about the uncertainty is that the uncertainty of
a measurement is equal to half of the smallest increment of the measurement
device. A common measuring device in teaching laboratories includes meter
sticks with milimetre gradations, corresponding to an uncertainty of $\pm
0.0005$ m. Recently, the use of smartphone accelerometers in teaching
laboratories has become more common communication with instructor (2020). The
time step in these measurements is significantly more precise than the
precision of a stopwatch and is not impacted by human reaction times. Common
measurements with accelerometers in recent courses at McGill University had
uncertainty of approximately $\pm 0.02$ s communication with instructor
(2020). Typical advice given to students is to time up to 10 periods of a
pendulum in order to reduce the relative error on the measurement, a process
replicated in the model used in this work. This study assumes that
measurements are normally distributed with standard deviation corresponding to
the uncertainty in timing alone. It further assumes that there are no
systematic biases in the measurements. Generated data is used for this study
in order to calculate an idealized limit without systematic bias and without
needing to account for complexities such as different pendulum lengths and a
variety of other effects already discussed. Generated data used in model
validation is termed “pseudodata” while generated data used for the final
inference is termed “proxy data”.

In order to determine when the small angle approximation becomes insufficient
using realistic measurements, realistic choices must be made. Different values
of the “true” gravitational acceleration and timing precision may be specified
when generating data; it is then possible to specify the initial angular
displacement and produce pseudo- and proxy data that reflects the precision of
data in introductory teaching laboratories.

With this choice for pseudodata generation, it is possible to plot raw model
predictions from the two methods of calculating the period of a simple
pendulum. This is shown in Fig. 1, where it is clear that the models become
increasingly differentiable as $\theta_{0}$ increases. This lends credence to
the study and suggests that as $\theta_{0}$ increases, progressively stronger
preference for a more complex model will be found.

Figure 1: Model predictions at fixed length $L=1$ m and gravitational
acceleration $g=9.8065$ $m/s^{2}$. Uncertainties shown include realistic
student measurement uncertainty of the pendulum length.

###  II.4 Defining the prior state of knowledge

As the initial angular displacement and period can be easily measured and the
length of the pendulum has been fixed, the only free model parameter in both
the equations using the small angle approximation and the exact expression for
the period of a pendulum is the gravitational acceleration $g$. In order to
not bias the result with existing knowledge about the solution, it is
important to explicitly construct a prior with physically-motivated but
relatively weakly-informative constraints. This allows for the model to
provide constraint to parameters via the likelihood rather than providing
artificially strong constraint via the prior. It is reasonable to assert _a
priori_ that for the period to be rigorously defined, $g$ cannot be zero,
cannot be negative, and cannot be infinite. It is also possible to ascribe
some degree of physical intuition to the problem. Gravitational acceleration
on the surface of a body scales with the mass of the body and gravitational
acceleration on the Moon is approximately 1.625 $m/s^{2}$ Hirt and
Featherstone (2012) and the surface gravity of Jupiter is 24.79 $m/s^{2}$ NASA
(2021). Thus, proceeding with confidence that the mass of the Earth is likely
between that of the Moon and Jupiter, a weakly-informative prior can be
constructed that allows for some probability of surprise. A common choice for
a weakly-informative prior on positive definite quantities is the inverse
gamma distribution, which can be tuned so that only 1% of the probability
falls below 1.625 $m/s^{2}$ and above 24.79 $m/s^{2}$. Thus, no positive real
value has been ruled out, although some values have been determined to be
unlikely. This prior expectation for $g$ is shown explicitly in Fig. 2.

Figure 2: The proposed prior probability density function (PDF) of the
gravitational acceleration $g$ measured using a simple pendulum used in a
teaching physics laboratory. This PDF is $P(H|I)$ in Eq. 1 where each value of
$g$ is a hypothesis $H$.

###  II.5 Prior predictive checks

Prior predictive checking is an important step in gaining an understanding of
model behavior Box (1980). In prior predictive checks, the computational model
is evaluated with draws from the prior distribution. This is then used to
assess the underlying model and determine the presence of non-trivial
features. In the case of the simple models shown here, few of these features
exist, but the step is still important in order to demonstrate a rigorous
modeling procedure. Predictions for the period with $g$ sampled according to
the prior distribution are shown in Fig. 3 while the normalized distribution
of these prior predictions is shown in Fig. 4.

Figure 3: Predictions for the period with $g$ sampled from the prior and
initial angular displacement sampled from a uniform distribution $U(0,\pi/4)$.
Figure 4: Normalized distribution of prior predictions for the period with
varying $g$.

These are compatible with the prior knowledge defined in Sec. II.4. It is
expected that measured periods of the pendulum are between approximately 1 and
5 seconds, while longer and shorter periods and corresponding values of $g$
are still theoretically possible, but disfavored. Additionally, the absolute
difference between the two models is small relative to the impact of the
gravitational acceleration, meaning that the result is highly sensitive to the
parameter $g$.

###  II.6 Model validation

Model validation, also known as closure testing or empirical coverage, is a
means of assessing whether the model has sufficient discriminating power to
constrain parameters given data. This is established by generating pseudodata
using the models themselves and seeing if the underlying truth can be
recovered. Pseudodata is generated with a gravitational acceleration
$g=3.72076$ $m/s^{2}$ corresponding to that found on Mars Hirt  _et al._
(2012) and a timing uncertainty of 0.05 s over 10 periods. Initial
displacements are chosen with arbitrary spacing,
$\theta_{0}=\left[0.05,0.2,0.35,0.4,\pi/5\right]$ radians, and calculations of
the corresponding periods are shown in Fig. 5.

Figure 5: Pseudodata used for model validation.

The models are then both calibrated to this closure psuedodata and the
posterior distributions are shown in Fig. 6, where it can be clearly seen that
the exact formula for the period is superior to the small angle approximation.
The posterior is evaluated numerically using Markov Chain Monte Carlo (MCMC).
For each parameter value considered, the likelihood (Eq. 2) is calculated and
the value of the posterior can be computed. While this is simple in a small
number of dimensions, this is impractical for high-dimensional applications
and computing the posterior directly becomes computationally inefficient. MCMC
algorithms are efficient sampling techniques such that the equilibrium
distribution of samples corresponds to a specified target distribution – in
this case, the posterior specified in Eq. 1. A chain of these samples can then
be analyzed as draws from the posterior. This study uses a numerical
implementation of MCMC in the ptemcee
package111https://github.com/willvousden/ptemcee Vousden  _et al._ (2015);
Foreman-Mackey  _et al._ (2013). It is possible to differentiate between the
exact model and the small angle posterior distributions, with the latter
biased to smaller values of $g$.

Figure 6: Posterior distributions for $g$ in model validation using the exact
formula for the period and the small angle approximation. The true value is
shown in green and vertical dashed lines denote the 16th, 50th, and 84th
percentile of the samples in increasing $g$.

It is important to consider the Markov Chain Monte Carlo behavior to ensure
that the chains converged to their target distributions and that samples are
indeed drawn from the target distribution. This is most simply analyzed by
considering autocorrelation, which quantifies how correlated samples are
within a certain “lag”. True independent samples of the posterior should be
uncorrelated. The autocorrelation for the exact and small angle approximation
MCMC chains are shown in Fig. 7, which show that the autocorrelation rapidly
drops to 0, showing that the chains are uncorrelated and are sampling the
target distribution successfully and are not influenced by their previous
location. This degree of uncorrelated MCMC samples is a common feature of the
0-temperature chain produced by the parallel tempering method of ptemcee.

Figure 7: Autocorrelation of MCMC chains for model validation. The rapid drop
and subsequent small-scale fluctuations around 0 are indicative of a
decorrelated chain, meaning that the chain is sampling from the target
distribution (the posterior).

The final check that must be performed before confidently progressing to the
final inference is to consider the posterior predictive distribution, Fig. 8.
Predictions can be made from the posterior distribution to form the posterior
predictive distribution, which helps determine if the model truly behaves as
anticipated by considering the distribution of results. As expected from Fig.
6, the small angle approximation attempts to fit the data by overshooting at
low initial angular displacements and undershooting at large ones while the
exact formula for the pendulum accurately captures the angular dependence.
With this consistent picture that matches expectations, it is now possible to
proceed to the inferential stage.

Figure 8: Violin plot of the prior and posterior predictive distributions for
$g$ using the exact formula for the period (upper panel) and the small angle
approximation (lower panel). The posterior distributions are those shown in
Fig. 6.

###  II.7 Inference

To demonstrate and validate this technique, proxy data for the pendulum with
$g=9.8065$ $m/s^{2}$ is generated at arbitrary intervals:
$\theta_{0}=[0.01,0.2,0.4,0.5,\pi/4]$ radians (Fig. 9). Proxy data is used
rather than student measurements as this study considers an idealized case
free from potential systematic or statistical errors arising from real-world
or finite-size effects. From the model validation study, larger
differentiation between the two models should be expected as the relative
uncertainty of the data is smaller. While proxy data is used in this study for
the purpose of attempting to determine a maximally-restrictive result for when
the small angle approximation to the period of a pendulum becomes
insufficient, in general experimental data is used for this stage.

Figure 9: Proxy data at selected $\theta_{0}$.

The posterior is again evaluated numerically using ptemcee. The posterior
distribution for the gravitational acceleration $g$ is shown in Fig. 10; the
exact model is able to correctly infer the gravitational acceleration with
high precision ($\pm\sim$ 0.001 $m/s^{2}$ in the 68% credible region) while
the small angle approximation achieves comparable precision but with a
remarkable $\sim 0.1$ $m/s^{2}$ bias, causing the two posterior distributions
to be broadly distinct.

Figure 10: Posterior distributions for $g$ using the exact formula for the
period and the small angle approximation. The true value is shown in green and
vertical dashed lines denote the 16th, 50th, and 84th percentile of the
samples in increasing $g$.

###  II.8 Posterior predictive checks

The source of this bias can be clearly seen in the posterior predictive
distribution (Fig. 11). This distribution is simply the posterior distribution
evaluated in each model of the period (small angle approximation or exact).
The predictions from the exact model calculation of the posterior are tightly
constrained to the data while the small angle model compensates for
undershooting at large angular displacement by introducing a bias and
overshooting small displacements and undershooting large displacements. A
common introductory physics experiment is to calculate $g$ using a simple
pendulum given measured vales of $L$ and the period. By using the small angle
approximation expression of the period, the bias on $g$ seen in Fig. 10 is not
accounted for by standard guidance on uncertainties.

Figure 11: Violin plot of the prior and posterior predictive distributions for
$g$ using the exact formula for the period (upper panel) and the small angle
approximation (lower panel). The posterior distributions are those shown in
Fig. 10.

This clearly establishes that it is possible to bias results by using the
small angle approximation of the period of a pendulum and that if one includes
large initial angular displacements, the impact of the approximation is
clearly detectable while accounting for realistic uncertainties.

###  II.9 Model selection

Now that it has been established that a difference is detectable with a finite
number of measurements in $\theta_{0}$, it is time to turn to the tools of
Bayesian model selection to quantify when a preference for the exact
expression for the period is weak, moderate, and strong.

Bayesian model selection is evaluated via the Bayes Factor. The Bayes factor
is the ratio of Bayes evidences (see Eq. 1) between two models.

$B_{01}=\frac{P(d|M_{0})}{P(d|M_{1})}$ (8)

where models $M_{i}$ are subscripted 0 and 1. The Bayes factor $B_{01}$ is
used to quantify model preference. A Bayes factor greater than one represents
an increase of support for $M_{0}$ relative to $M_{1}$ and directly
corresponds to the odds of $M_{0}$:$M_{1}$. A more comprehensive introduction
to Bayesian model selection and its advantage over frequentist methods may be
found in Trotta (2008). From here, the exact formula for the period is
subscripted 0 and the small angle approximation is subscripted 1.

Empirical scales are used to determine when there is weak, moderate, and
strong evidence for $M_{0}$ vs. $M_{1}$. In this work, the Jeffreys’ Scale is
used (Table. 1).

$\left|\ln B_{01}\right|$ | Odds | Probability | Strength of evidence
---|---|---|---
$<1.0$ | $\leq 3:1$ | $<0.750$ | Inconclusive
1.0 | $\sim 3:1$ | 0.750 | Weak evidence
2.5 | $\sim 12:1$ | 0.923 | Moderate evidence
5.0 | $\sim 150:1$ | 0.993 | Strong evidence
Table 1: The Jeffreys’ Scale, reproduced from Trotta (2008).

Data is generated for the model selection study as follows. It is unreasonable
to assume that realistic introductory physics students would be able to make
measurements at increments finer than $1^{\circ}$ with a one-meter-long
pendulum or at more than 12 increments in initial angular displacement
$\theta_{0}$. Typically, students make fewer measurements in teaching
laboratories and thus this study represents an attempt to find an idealized,
maximally-restrictive constraint on initial angular displacement $\theta_{0}$
for the use of the applicability of the small angle approximation expression
of the period.

12 proxy data points are generated if the maximum $\theta_{0}$ is larger than
12∘. If the maximum $\theta_{0}$ is less than 12∘, data is generated in 1∘
increments. The Bayes factor is then calculated and it is possible to
interpolate between calculations with their uncertainty to determine at what
$\rm{Max}[\theta_{0}]$ the evidence for the exact model becomes weak,
moderate, or strong.

Figure 12: $\ln$ Bayes Factor demonstrating preference for the exact
expression for the period of the pendulum over the small angle approximation.
Horizontal lines denote thresholds for weak, moderate, and strong preference.

The $\ln$ Bayes factor is shown in Fig. 12. As anticipated, as the maximum
initial angular displacement $\theta_{0}$ increases, there is progressively
stronger preference for the exact expression. Conversely, at small
$\theta_{0}$, the data does not exhibit a preference between the two models
and it may be stated confidently that the small angle approximation is valid
and the models are indistinguishable at the level of the experimental
precision this work considers. Preference for the exact calculation at maximum
initial angular displacement $\theta_{0}$ is given in Table. 2

Strength of evidence | $\left|\ln B_{01}\right|$ |  $\theta_{0}$ (rad) |  $\theta_{0}$ (deg)
---|---|---|---
Weak | 1.0 | $0.192^{+0.008}_{-0.009}$ | $11.0\pm 0.5$
Moderate | 2.5 | $0.238\pm 0.005$ | $13.6\pm 0.3$
Strong | 5.0 | $0.282\pm 0.003$ | $16.2\pm 0.2$
Table 2: Model preference thresholds where subscript 0 denotes the exact model
and subscript 1 denotes the small angle approximation. Conversion to degrees
is provided for convenience.

##  III Conclusions

Using the techniques of Bayesian inference, it is possible to set rigorous
data-driven constraints on the use of the small angle approximation in
introductory physics experiments. To my knowledge, this is the first instance
in the literature of modern computational and statistical techniques being
applied to this problem and exemplifies a consistent way to apply the
techniques of research physics to common pedagogical approximations. The
heuristic guidance of constraining simple pendula to $\sim 15^{\circ}$ is
between the inferential results for moderate and strong preference. The data
does not exhibit a strong preference for a more complex calculation until
$16.2^{\circ}\pm 0.2^{\circ}$ under highly idealized conditions and for a
pendulum with length 1m. This suggests that the current guidance is consistent
with rigorous statistical expectations, but may be overly-restrictive for
shorter pendula in which the same absolute timing uncertainty contributes to
larger relative error in the period. Future works should thus continue with
the recommendation of restricting initial angular displacement to $\sim
15^{\circ}$ when employing the small angle approximation formula for the
period.

As this analysis was performed with a specific set of model assumptions
regarding the length of the pendulum, the timing accuracy of the pseudo and
proxy data, and no nonconservative forces, it represents general guidance for
laboratory manuals and textbooks.

##  IV Acknowledgements

This work was supported in part by the Natural Sciences and Engineering
Research Council of Canada. Thanks are due to C. Gale, N. Provatas, M. Singh,
N. Fortier, and R. Yazdi for critical reading of the manuscript and to D.
Everett for both a critical reading of the manuscript and detailed discussions
on the implementation of Bayesian techniques and their interpretation. R.
Furnstahl’s resources on Bayesian inference were instrumental in the
development of this work. C. Gale, S. Jeon, and K. Ragan provided guidance,
useful discussions, and support to this study.

## References

  * Lima and Arun (2006) F. M. S. Lima and P. Arun, “An accurate formula for the period of a simple pendulum oscillating beyond the small angle regime,” American Journal of Physics 74, 892–895 (2006), https://doi.org/10.1119/1.2215616 .
  * Blais (2020) Brian S. Blais, “Model comparison in the introductory physics laboratory,” The Physics Teacher 58, 209–213 (2020), https://doi.org/10.1119/1.5145420 .
  * Giancoli (2016) Douglas C Giancoli,  _Physics: principles with applications_ (Boston: Pearson, 2016).
  * Hinrichsen (2020) Peter F Hinrichsen, “Review of approximate equations for the pendulum period,” European Journal of Physics 42, 015005 (2020).
  * Smith  _et al._ (2020) Rory J E Smith, Gregory Ashton, Avi Vajpeyi, and Colm Talbot, “Massively parallel Bayesian inference for transient gravitational-wave astronomy,” Monthly Notices of the Royal Astronomical Society 498, 4492–4502 (2020), https://academic.oup.com/mnras/article-pdf/498/3/4492/33798799/staa2483.pdf .
  * Everett  _et al._ (2020) D. Everett _et al._ (JETSCAPE), “Multi-system Bayesian constraints on the transport coefficients of QCD matter,” (2020), arXiv:2011.01430 [hep-ph] .
  * von Toussaint (2011) Udo von Toussaint, “Bayesian inference in physics,” Rev. Mod. Phys. 83, 943–999 (2011).
  * Trotta (2008) Roberto Trotta, “Bayes in the sky: Bayesian inference and model selection in cosmology,” Contemporary Physics 49, 71–104 (2008), https://doi.org/10.1080/00107510802066753 .
  * Gelman  _et al._ (2020) Andrew Gelman, Aki Vehtari, Daniel Simpson, Charles C. Margossian, Bob Carpenter, Yuling Yao, Lauren Kennedy, Jonah Gabry, Paul-Christian Bürkner, and Martin Modrák, “Bayesian workflow,” (2020), arXiv:2011.01808 [stat.ME] .
  * Betancourt (2019) M. Betancourt, “Falling (in love with principled modeling),”  (2019).
  * Schery (1976) Stephen D. Schery, “Design of an inexpensive pendulum for study of large‐angle motion,” American Journal of Physics 44, 666–670 (1976), https://doi.org/10.1119/1.10352 .
  * Virtanen  _et al._ (2020) Pauli Virtanen _et al._ , “Scipy 1.0: Fundamental algorithms for scientific computing in python,” Nature Methods (2020).
  * communication with instructor (2020) Private communication with instructor, “McGill University Physics 101 lab manual,” (2020).
  * Hirt and Featherstone (2012) C. Hirt and W.E. Featherstone, “A 1.5km-resolution gravity field model of the Moon,” Earth and Planetary Science Letters 329-330, 22–30 (2012).
  * NASA (2021) NASA, “By the numbers: Jupiter,”  (2021).
  * Box (1980) G. E. Box, “Sampling and Bayes’ inference in scientific modelling and robustness,” (1980).
  * Hirt  _et al._ (2012) C. Hirt, S.J. Claessens, M. Kuhn, and W.E. Featherstone, “Kilometer-resolution gravity field of Mars: MGM2011,” Planetary and Space Science 67, 147–154 (2012).
  * Vousden  _et al._ (2015) W. D. Vousden, W. M. Farr, and I. Mandel, “Dynamic temperature selection for parallel tempering in Markov chain Monte Carlo simulations,” Monthly Notices of the Royal Astronomical Society 455, 1919–1937 (2015), https://academic.oup.com/mnras/article-pdf/455/2/1919/18514064/stv2422.pdf .
  * Foreman-Mackey  _et al._ (2013) Daniel Foreman-Mackey, David W. Hogg, Dustin Lang, and Jonathan Goodman, “emcee: The MCMC hammer,” Publications of the Astronomical Society of the Pacific 125, 306–312 (2013).