††thanks<EMAIL_ADDRESS>
# Optical sensing of magnons via the magnetoelastic displacement

Zhi-Yuan Fan  Interdisciplinary Center of Quantum Information, State Key
Laboratory of Modern Optical Instrumentation, and Zhejiang Province Key
Laboratory of Quantum Technology and Device, Department of Physics, Zhejiang
University, Hangzhou 310027, China  Rui-Chang Shen  Interdisciplinary Center
of Quantum Information, State Key Laboratory of Modern Optical
Instrumentation, and Zhejiang Province Key Laboratory of Quantum Technology
and Device, Department of Physics, Zhejiang University, Hangzhou 310027, China
Yi-Pu Wang  Interdisciplinary Center of Quantum Information, State Key
Laboratory of Modern Optical Instrumentation, and Zhejiang Province Key
Laboratory of Quantum Technology and Device, Department of Physics, Zhejiang
University, Hangzhou 310027, China  Jie Li  Interdisciplinary Center of
Quantum Information, State Key Laboratory of Modern Optical Instrumentation,
and Zhejiang Province Key Laboratory of Quantum Technology and Device,
Department of Physics, Zhejiang University, Hangzhou 310027, China  J. Q. You
Interdisciplinary Center of Quantum Information, State Key Laboratory of
Modern Optical Instrumentation, and Zhejiang Province Key Laboratory of
Quantum Technology and Device, Department of Physics, Zhejiang University,
Hangzhou 310027, China

###### Abstract

We show how to measure a steady-state magnon population in a magnetostatic
mode of a ferromagnet or ferrimagnet, such as yttrium iron garnet. We adopt an
optomechanical approach and utilize the magnetoelasticity of the ferromagnet.
The magnetostrictive force dispersively couples magnons to the deformation
displacement of the ferromagnet, which is proportional to the magnon
population. By further coupling the mechanical displacement to an optical
cavity that is resonantly driven by a weak laser, the magnetostrictively
induced displacement can be sensed by measuring the phase quadrature of the
optical field. The phase shows an excellent linear dependence on the magnon
population for a not very large population, and can thus be used as a
‘magnometer’ to measure the magnon population. We further study the effect of
thermal noises, and find a high signal-to-noise ratio even at room
temperature. At cryogenic temperatures, the resolution of magnon excitation
numbers is essentially limited by the vacuum fluctuations of the phase, which
can be significantly improved by using a squeezed light.

##  I Introduction

Ever since the successful demonstration of the strong coupling between
microwave cavity photons and ferrimagnetic magnons in yttrium iron garnet
(YIG) S1 ; S2 ; S3 , the field of cavity magnonics has attracted considerable
attention and achieved significant progress NakaRev ; NatRev ; RMP . Such a
strongly coupled system has become a new platform for studying rich and
stimulating phenomena belonging to cavity quantum electrodynamics. This
benefits largely from the distinct advantages of YIG, such as the high spin
density, low damping rate, and rich nonlinear interactions between excitations
in the crystal.

Many studies in cavity magnonics require the estimation of the magnon
population in a magnetostatic mode of a ferromagnet or ferrimagnet. For
example, it can be used to determine the effective dispersive coupling
strength in either the magnon-qubit interaction Nak17 ; Nak20S  or the magnon-
phonon interaction JL18 ; JL19 . The method of measuring the transmission (or
reflection) spectrum of the microwave cavity in the conventional cavity-magnon
system S1 ; S2 ; S3  is not able to infer the magnon population because their
linear beamsplitter coupling is independent of the magnon population. One can,
however, implement state tomography of microwave photons that have a
beamsplitter interaction with magnons. In this way, the magnon population can
be determined, since the population is just the square of the amplitude in the
phase space. The tomography of magnonic states has been realized in Ref. prb21
. There are approaches for measuring a large magnon population based on, e.g.,
electromagnetic induction Melkov , the inverse spin-Hall effect invSH1 ;
invSH2 ; invSH3 ; prb21 , Brillouin light scattering BLS1 ; BLS2 , etc.
Alternatively, one may consider utilizing the magnon Kerr nonlinearity, where
the Kerr effect leads to a frequency shift of the magnon mode dependent on the
magnon population YP16 ; YP18 ; YP21 . Thus the magnon population can be
inferred by measuring the frequency shift. This approach works also for a
large magnon population (under a strong pump). For a small magnon population,
it can be realized by dispersively coupling magnons to a superconducting qubit
Nak17 ; Nak20S ; Nak20L , where either the qubit frequency Nak17  or dephasing
Nak20L  depends on the magnon population. Therefore, by measuring the
frequency shift or coherence time of the qubit, the magnon excitations can be
resolved. These approaches can achieve very high resolutions of magnons, but
typically work at very low temperatures.

Here, we introduce a new approach by adopting the magnetostrictive
(magnetoelastic) coupling between magnons and the vibrational motion of a
ferrimagnet, e.g. YIG. The magnetostrictive force leads to the geometric
deformation of the ferrimagnet, forming vibrational modes (phonons) Kittle .
We consider the situation where the phonon frequency is much smaller than the
magnon frequency JL18 ; Tang ; Davis , such that their dispersive interaction
becomes dominant, leading to the deformation displacement proportional to the
magnon population. This is analogous to the displacement of a mechanical
oscillator induced by the radiation pressure in optomechanics OMrmp .
Therefore, by measuring the displacement one can infer the magnon population.
To do it, we adopt an optomechanical approach by coupling the displacement to
an optical cavity. When the cavity is resonantly driven, a high-sensitivity
measurement of the mechanical displacement can be realized by measuring the
phase of the cavity field. To minimize the disturbance to the displacement due
to radiation pressure, we use a weak probe light. This guarantees that the
frequency shift of the optical cavity is mainly caused by the
magnomechanically induced displacement. For a not very large magnon
population, we find an excellent linear dependence of the optical phase on the
magnon population. Therefore, in this linear regime the optical phase can act
as a meter for the magnon population. We further study the effect of thermal
noises on our protocol, and find that the signal-to-noise ratio (SNR) can be
very high even at room temperature. The resolution of magnons is essentially
limited by the vacuum fluctuations of the phase when working at cryogenic
temperatures, e.g., tens of millikelvin. The resolution can be significantly
improved by feeding a squeezed light into the cavity, which can suppress the
phase noise to be well below the standard quantum limit.

This article is organized as follows. In Sec. II, we introduce our model and
provide its Hamiltonian and quantum Langevin equations. We then solve the
equations and obtain steady-state solutions for the classical averages, which
lead to the main result of the work: the magnon population can be measured by
the optical phase based on their linear dependence. In Sec. III, we study the
effect of thermal noises on our protocol, which lead to phase fluctuations and
a finite SNR. We derive the analytical expression of the phase noise and show
that a high SNR can still be achieved in the presence of thermal noises, even
at room temperature. We discuss how to improve the SNR and the resolution of
magnons. Finally, we summarize in Sec. IV.

Figure 1: Sketch of the protocol for optically measuring magnon population in
a ferromagnet. By using a homodyne setup, the magnetostrictively induced
displacement, proportional to the magnon population, is read out in the phase
quadrature of the optical field. BS: beamsplitter. See the text for a detailed
description.

##  II The protocol

We consider a magnomechanical system, e.g., in a YIG sphere note , which
consists of a magnon mode (e.g., the uniform-precession Kittel mode Kittel ),
and a phonon mode of deformation vibration, as sketched in Fig. 1. The magnon
mode represents the collective motion (spin wave) of a large number of spins
in the ferrimagnet, and is activated by placing the ferrimagnet in a uniform
bias magnetic field and applying, e.g., via a loop coil Naka16 , a microwave
drive field with its magnetic component perpendicular to the bias field. The
magnon mode couples to a deformation phonon mode of the ferrimagnet via the
magnetostrictive interaction. We consider the situation where the frequency of
the phonon mode is much lower than that of the magnon mode JL18 ; Tang ; Davis
, such that they couple via a radiation pressure-like dispersive interaction.
Such a nonlinear coupling has been recognized as a cornerstone of many quantum
protocols JieNJP ; Tan ; Kong ; Ding ; JL20 ; Davis20 ; Yang ; QST ; NSR ; HFW
; TB ; Jing ; JHW ; HJ ; CLi . Relevant experiments have demonstrated
magnomechanically induced transparency/absorption Tang  and mechanical
cooling/lasing Davis .

We aim to probe the magnons in the ferrimagnet, and particularly measure its
steady-state population. To this end, we adopt an optomechanical approach by
coupling the magnetostrictively induced displacement to an optical cavity.
This can be realized, e.g., by attaching a high-reflectivity mirror pad SG  to
the surface of the ferrimagnet. The mirror is so small that it does not
appreciably affect the mechanical properties of the ferrimagnet. Therefore,
through this way the magnomechanical displacement can be probed by light
exploiting the optomechanical interaction, which is then used to determine the
magnon population. The corresponding Hamiltonian of the magnomechanical system
combined with the optical cavity reads

$H=H_{0}+H_{\rm int}+H_{d},$ (1)

where

$\begin{split}H_{0}/\hbar=&\,\omega_{m}m^{\dagger}m+\frac{\omega_{b}}{2}\left(q^{2}+p^{2}\right)+\omega_{c}c^{\dagger}c,\\\
H_{\rm int}/\hbar=&\,g_{mb}m^{\dagger}mq-g_{cb}c^{\dagger}cq,\\\
H_{d}/\hbar=&\,i\Omega\left(m^{\dagger}e^{-i\omega_{0}t}-me^{i\omega_{0}t}\right)+iE\left(c^{\dagger}e^{-i\omega_{L}t}-ce^{i\omega_{L}t}\right),\end{split}$
(2)

are the free Hamiltonian, the interaction Hamiltonian, and the Hamiltonian of
the drives of the magnon and cavity modes, respectively. $m$ and $c$
($m^{\dagger}$ and $c^{\dagger}$) are the annihilation (creation) operators of
the magnon and cavity modes, respectively, satisfying the commutation relation
$[j,j^{\dagger}]=1$ ($j=m,c$). $q$ and $p$ ($[q,p]=i$) denote the
dimensionless position and momentum of the deformation vibrational mode,
modeled as a mechanical oscillator, which simultaneously couples to the magnon
mode via the magnetostrictive interaction and to the cavity field via the
radiation pressure. Both interactions are of the nonlinear dispersive type,
and the corresponding bare coupling rates are $g_{mb}$ and $g_{cb}$.
$\omega_{m}$, $\omega_{b}$ and $\omega_{c}$ are the resonant frequencies of
the magnon, mechanical and cavity modes, respectively. The frequency of the
magnon mode can be adjusted by varying the bias magnetic field $H_{0}$ via
$\omega_{m}=\gamma H_{0}$, with the gyromagnetic ratio $\gamma/2\pi=28$ GHz/T
for the YIG. In typical cavity magnonic experiments S1 ; S2 ; S3 ,
$\omega_{m}/2\pi$ is about 10 GHz. The magnon mode is driven by a microwave
magnetic field with amplitude $B_{0}$ and frequency $\omega_{0}$. The
corresponding Rabi frequency $\Omega=\frac{\sqrt{5}}{4}\gamma\sqrt{N}B_{0}$
for a YIG sphere JL18 , with $N$ the total number of spins. The cavity is
weakly driven by a laser with frequency $\omega_{L}$ and power $P_{L}$, and
the corresponding coupling strength
$E=\sqrt{\kappa_{c}P_{L}/\hbar\omega_{L}}$, with $\kappa_{c}$ the cavity decay
rate.

By including the dissipation and input noise of each mode, and working in the
interaction picture with respect to
$\hbar\omega_{0}m^{\dagger}m+\hbar\omega_{L}c^{\dagger}c$, the quantum
Langevin equations (QLEs) of the whole system are given by

$\begin{split}\dot{m}=&-\left(i\Delta_{m}+\frac{\kappa_{m}}{2}\right)m-ig_{mb}mq+\Omega+\sqrt{\kappa_{m}}m^{\rm
in},\\\
\dot{c}=&-\left(i\Delta_{c}+\frac{\kappa_{c}}{2}\right)c+ig_{cb}cq+E+\sqrt{\kappa_{c}}c^{\rm
in},\\\
\dot{q}=&\omega_{b}p,\,\,\,\dot{p}=-\omega_{b}q-\gamma_{b}p-g_{mb}m^{\dagger}m+g_{cb}c^{\dagger}c+\xi,\end{split}$
(3)

where $\Delta_{m}=\omega_{m}-\omega_{0}$ ($\Delta_{c}=\omega_{c}-\omega_{L}$)
is the detuning of magnon (cavity) mode with respect to its drive field, and
$\kappa_{m}$, $\kappa_{c}$ and $\gamma_{b}$ ($m^{\rm in}$, $c^{\rm in}$ and
$\xi$) are the dissipation rates (input noises) of the magnon, cavity, and
mechanical modes, respectively. The input noise operators have a zero
expectation value and the following nonzero correlation functions GC :
$\left\langle m^{\rm in}(t)m^{\rm
in\dagger}(t^{\prime})\right\rangle=\left[n_{m}(\omega_{m})+1\right]\delta(t-t^{\prime})$,
$\left\langle m^{\rm in\dagger}(t)m^{\rm
in}(t^{\prime})\right\rangle=n_{m}(\omega_{m})\delta(t-t^{\prime})$, and
$\left\langle c^{\rm in}(t)c^{\rm
in\dagger}(t^{\prime})\right\rangle=\delta(t-t^{\prime})$, and a Markovian
approximated $\delta$-correlated mechanical noise:
$\left\langle\xi(t)\xi(t^{\prime})+\xi(t^{\prime})\xi(t)\right\rangle/2\simeq\gamma_{b}\left[2n_{b}(\omega_{b})+1\right]\delta(t-t^{\prime})$,
which is the case for a large mechanical quality factor
$Q=\omega_{b}/\gamma_{b}\gg 1$ Kac . The mean thermal excitation number
$n_{j}(\omega_{j})=1/[\mathrm{exp}(\hbar\omega_{j}/k_{B}T)-1]$ ($j=m,b$) at an
environmental temperature $T$.

Under continuous drives, the system will evolve to a steady state under the
parameters constrained by the stability condition (which will be discussed in
Sec. III). When the amplitudes of the magnon and cavity modes are sufficiently
large, $\left|\left\langle m\right\rangle\right|,\,\left|\left\langle
c\right\rangle\right|\gg 1$, one can linearize the nonlinear dynamics around
the steady-state averages by writing each mode operator as a classical average
plus a quantum fluctuation operator, $k=\left\langle k\right\rangle+\delta k$
($k=m,c,q,p$), and neglecting the small second-order fluctuation terms. As a
result, the QLEs (3) are separated into two sets of linear equations for
classical averages and quantum fluctuations, respectively. By solving the set
of equations for classical averages, we obtain

$\begin{split}\left\langle
m\right\rangle&=\frac{\Omega}{i\tilde{\Delta}_{m}+\frac{\kappa_{m}}{2}},\,\,\,\,\,\,\,\,\,\,\,\,\,\left\langle
c\right\rangle=\ \frac{E}{i\tilde{\Delta}_{c}+\frac{\kappa_{c}}{2}},\\\ \
\left\langle q\right\rangle=\,&\Big{(}g_{cb}\left|\left\langle
c\right\rangle\right|^{2}-g_{mb}\left|\left\langle
m\right\rangle\right|^{2}\Big{)}/\omega_{b},\,\,\,\,\,\left\langle
p\right\rangle=0,\end{split}$ (4)

where $\tilde{\Delta}_{m}=\Delta_{m}+g_{mb}\left\langle q\right\rangle$ and
$\tilde{\Delta}_{c}=\Delta_{c}-g_{cb}\left\langle q\right\rangle$ are the
effective detunings. We consider a resonant laser drive,
$\omega_{L}\simeq\omega_{c}$ ($\Delta_{c}\simeq 0$), which corresponds to an
optimal situation for realizing high-sensitivity detection of the mechanical
position by measuring the phase of the cavity field JOSAB . Under the resonant
drive, the Stokes and anti-Stokes scattering probabilities are equal and the
optomechanical interaction resembles a quantum nondemolition interaction
Grangier . Furthermore, the cavity is weakly driven, leading to
$g_{cb}\left|\left\langle c\right\rangle\right|^{2}\ll
g_{mb}\left|\left\langle m\right\rangle\right|^{2}$, such that the
displacement is dominantly induced by magnetostriction, i.e., $\left\langle
q\right\rangle\simeq-\frac{g_{mb}}{\omega_{b}}\left|\left\langle
m\right\rangle\right|^{2}$, which shows a linear dependance of $\left\langle
q\right\rangle$ on the magnon population $N_{m}=|\left\langle
m\right\rangle|^{2}$. This also ensures that the radiation pressure will yield
a negligible backaction on the magnon population via the mediation of the
mechanical oscillator.

Figure 2: Steady-state average of (a) the phase quadrature $\left|\left\langle
Y_{c}\right\rangle\right|$ and (b) the amplitude quadrature $\left\langle
X_{c}\right\rangle$ versus magnon population $\left|\left\langle
m\right\rangle\right|^{2}$. See the text for other parameters.

The tight connection between the mechanical displacement and the optical phase
(under a resonant drive) offers the possibility to observe this linear
dependance in the phase quadrature by homodyning the cavity output field (see
Fig. 1). This can be seen from the expressions of the phase and amplitude
quadratures of the cavity field, given by

$\left\langle Y_{c}\right\rangle={-}\\!\sqrt{2}\
\frac{E\tilde{\Delta}_{c}}{\tilde{\Delta}_{c}^{2}+\frac{\kappa_{c}^{2}}{4}},\,\,\,\,\,\,\left\langle
X_{c}\right\rangle=\sqrt{2}\
\frac{\frac{1}{2}E\kappa_{c}}{{\tilde{\Delta}_{c}}^{2}+\frac{\kappa_{c}^{2}}{4}},$
(5)

where $\tilde{\Delta}_{c}=-g_{cb}\left\langle
q\right\rangle\simeq\frac{g_{cb}\,g_{mb}}{\omega_{b}}\left|\left\langle
m\right\rangle\right|^{2}$, implying that the cavity frequency shift is mainly
caused by the magnetostrictively induced displacement. Under the condition
that this frequency shift is much smaller than the cavity linewidth
$\tilde{\Delta}_{c}\ll\kappa_{c}$, we can approximately obtain

$\left\langle Y_{c}\right\rangle\simeq{-}4\\!\sqrt{2}\
\frac{E\tilde{\Delta}_{c}}{\kappa_{c}^{2}}\simeq{-}4\\!\sqrt{2}\,\frac{Eg_{cb}g_{mb}}{\kappa_{c}^{2}\omega_{b}}\left|\left\langle
m\right\rangle\right|^{2},\\\ $ (6)

and a constant amplitude average $\left\langle
X_{c}\right\rangle=\frac{2\sqrt{2}E}{\kappa_{c}}$. Clearly, Eq. (6) shows a
linear dependence of the phase $\left\langle Y_{c}\right\rangle$ on the magnon
population for a fixed laser power (more precisely, $E$). Though the condition
$\tilde{\Delta}_{c}\ll\kappa_{c}$ leads to an excellent linear dependence, it
places a limit on the maximum number of the magnon population that can be
probed linearly. For a sufficiently large population (causing
$\tilde{\Delta}_{c}$ to be comparable with $\kappa_{c}$), a nonlinear
dependence starts to emerge in both $\left\langle Y_{c}\right\rangle$ and
$\left\langle X_{c}\right\rangle$. This is seen in Figs. 2(a) and 2(b).
Although in this nonlinear regime, one can measure both the phase and
amplitude quadratures to determine the magnon population, we focus on the
linear regime where one only needs to measure the phase quadrature, and the
dependance is rather straightforward. We have used the following parameters
for Fig. 2: $\omega_{b}/2\pi=10$ MHz, $g_{mb}/2\pi=1$ Hz,
$\kappa_{c}/2\pi=100$ kHz, $g_{cb}/2\pi=10$ Hz, and a laser with power
$P_{L}=1$ $\mu$W and wavelength $\lambda_{L}=1064$ nm.

Figure 3: Steady-state $\left|\left\langle Y_{c}\right\rangle\right|$ versus
magnon population $\left|\left\langle m\right\rangle\right|^{2}$ at different
laser powers: (a) and (d) $P_{L}=10^{-7}$ W, (b) $P_{L}=10^{-8}$ W, (c)
$P_{L}=10^{-9}$ W. We take $\kappa_{c}/2\pi=100$ kHz in (a)-(c) and
$\kappa_{c}/2\pi=1$ MHz in (d). The rest of the parameters are the same as in
Fig.2.

Note that there is also a lower limit of measuring range of the magnon
population determined by the perturbation condition $g_{cb}\left|\left\langle
c\right\rangle\right|^{2}\,{\ll}\,g_{mb}\left|\left\langle
m\right\rangle\right|^{2}$ for a fixed laser power. In Figs. 3(a)-3(c), we
show the phase average $\left|\left\langle Y_{c}\right\rangle\right|$ versus
the magnon population for various laser powers. For each power, the measuring
range is chosen to meet both the conditions for perturbation and the linear
dependence. It shows that a larger power yields a larger $\left|\left\langle
Y_{c}\right\rangle\right|$ (which increases the SNR, to be discussed in Sec.
III), but the measuring range is narrowed. In contrast, increasing
$\kappa_{c}$ broadens the measuring range, but meantime reduces
$\left|\left\langle Y_{c}\right\rangle\right|$ (thus the SNR), as revealed by
comparing Figs. 3(a) and 3(d). So, there is a tradeoff between the SNR and the
measuring range by changing $P_{L}$ or $\kappa_{c}$. The two conditions can be
written in the following compatible form:

$\frac{g_{cb}^{2}}{\omega_{b}}\left|\left\langle
c\right\rangle\right|^{2}\ll\frac{g_{cb}\,g_{mb}}{\omega_{b}}\left|\left\langle
m\right\rangle\right|^{2}\ll\kappa_{c},$ (7)

which sets an upper limit for the coupling strength
$E\lll\kappa_{c}\sqrt{\kappa_{c}\omega_{b}}/\big{(}2g_{cb}\big{)}$, and thus
for the laser power. This leads to an upper bound for $\left|\left\langle
Y_{c}\right\rangle\right|/N_{m}$ from the result of Eq. (6), i.e.,

$\frac{\left|\left\langle Y_{c}\right\rangle\right|}{N_{m}}\lll
2\sqrt{2}\frac{g_{mb}}{\sqrt{\kappa_{c}\omega_{b}}}.$ (8)

It determines the highest possible SNR (for a given phase noise) and thus the
resolution of magnons in our protocol. The upper bound suggests adopting a
smaller cavity linewidth and a vibrational mode with a lower frequency but a
larger magnon-phonon (dispersive) coupling.

##  III Effect of thermal noises and signal-to-noise ratio

In the preceding section, we have established a linear dependence of the
steady-state phase average on the magnon population. Although the phase
average can be considered as the signal, there is also phase noise due to the
presence of thermal noises from the environment. In particular, the phonon
mode has a much lower frequency (typically in MHz for hundreds of microns
sized YIG spheres Tang ; Davis ), and thus possesses the dominate thermal
noise of the system. The corresponding SNR thus becomes a key indicator, which
eventually determines the resolution of magnons of our approach. In what
follows, we investigate the phase noise and the SNR in the steady state, and
demonstrate that our protocol is still effective under substantial thermal
noises.

The set of equations for quantum fluctuations can be written, in the
quadrature form, as

$\begin{split}\delta\dot{X}_{c}&=-\frac{\kappa_{c}}{2}\delta
X_{c}+\tilde{\Delta}_{c}\delta Y_{c}+\sqrt{\kappa_{c}}X_{c}^{\rm in},\\\
\delta\dot{Y}_{c}&=-\tilde{\Delta}_{c}\delta X_{c}-\frac{\kappa_{c}}{2}\delta
Y_{c}+G_{c}\delta q+\sqrt{\kappa_{c}}Y_{c}^{\rm in},\\\
\delta\dot{X}_{m}&=-\frac{\kappa_{m}}{2}\delta X_{m}+\tilde{\Delta}_{m}\delta
Y_{m}+{\rm Im}G_{m}\delta q+\sqrt{\kappa_{m}}X_{m}^{\rm in},\\\
\delta\dot{Y}_{m}&=-\tilde{\Delta}_{m}\delta X_{m}-\frac{\kappa_{m}}{2}\delta
Y_{m}-{\rm Re}G_{m}\delta q+\sqrt{\kappa_{m}}Y_{m}^{\rm in},\\\
\delta\dot{q}&=\omega_{b}\delta p,\\\ \delta\dot{p}&=-\omega_{b}\delta
q\,{-}\gamma_{b}\delta p\,{-}{\rm Re}G_{m}\delta X_{m}\,{-}{\rm Im}G_{m}\delta
Y_{m}\,{+}G_{c}\delta X_{c}\,{+}\xi,\\\ \end{split}$ (9)

where the quantum fluctuations of the quadratures are defined as $\delta
X_{k}\,\,{=}\,\,\frac{1}{\sqrt{2}}\left(\delta k+\delta k^{\dagger}\right)$,
$\delta Y_{k}\,\,{=}\,\,\frac{i}{\sqrt{2}}\left(\delta k^{\dagger}-\delta
k\right)$, and the corresponding quadratures of input noises $X_{k}^{\rm
in}=\frac{1}{\sqrt{2}}\left(k^{\rm in}+k^{\rm in\dagger}\right)$, $Y_{k}^{\rm
in}\,\,{=}\,\,\frac{i}{\sqrt{2}}\left(k^{\rm in\dagger}-k^{\rm in}\right)$
($k\,\,{=}\,\,c,m$). We consider a real optomechanical coupling
$G_{c}\,{=}\,\sqrt{2}g_{cb}\left\langle c\right\rangle$, since $\left\langle
c\right\rangle\,{\simeq}\,\left\langle c\right\rangle^{*}$ when
$\tilde{\Delta}_{c}\,\,{\ll}\,\,\kappa_{c}$, and a complex magnomechanical
coupling $G_{m}\,{=}\,\sqrt{2}g_{mb}\left\langle m\right\rangle$, as in
general $\left\langle m\right\rangle\,{\neq}\,\left\langle
m\right\rangle^{*}$. The equations (9) can be rewritten in a compact matrix
form as

$\dot{v}(t)={\cal A}v(t)+n(t),$ (10)

where $v(t)=\left[\delta X_{c}(t),\,\delta Y_{c}(t),\,\delta X_{m}(t),\,\delta
Y_{m}(t),\,\delta q(t),\,\delta p(t)\right]^{T}$,
$n(t)=\left[\\!\sqrt{\kappa_{c}}X^{\rm in}_{c}(t),\sqrt{\kappa_{c}}Y^{\rm
in}_{c}(t),\sqrt{\kappa_{m}}X^{\rm in}_{m}(t),\sqrt{\kappa_{m}}Y^{\rm
in}_{m}(t),0,\xi(t)\right]^{T}$, and the drift matrix ${\cal A}$ is given by

${\cal A}=\begin{pmatrix}-\frac{\kappa_{c}}{2}&\tilde{\Delta}_{c}&0&0&0&0\\\
-\tilde{\Delta}_{c}&-\frac{\kappa_{c}}{2}&0&0&G_{c}&0\\\
0&0&-\frac{\kappa_{m}}{2}&\tilde{\Delta}_{m}&\textup{Im}G_{m}&0\\\
0&0&-\tilde{\Delta}_{m}&-\frac{\kappa_{m}}{2}&-\textup{Re}G_{m}&0\\\
0&0&0&0&0&\omega_{b}\\\
G_{c}&0&-\textup{Re}G_{m}&-\textup{Im}G_{m}&-\omega_{b}&-\gamma_{b}\end{pmatrix}.$
(11)

The system becomes stable when $t\to\infty$ if all the eigenvalues of the
drift matrix ${\cal A}$ have negative real parts. This is equivalent to the
stability condition obtained from the Routh-Hurwitz criterion RH , but the
inequalities become quite involved for the present tripartite system. All of
the results presented in this work satisfy this condition and are thus in the
steady state.

The equations (9) can be solved conveniently in the frequency domain by taking
the Fourier transform of each equation. After some algebra, we obtain the
following solution for the phase fluctuation $\delta Y_{c}(\omega)$:

$\begin{split}&\delta
Y_{c}(\omega)=\chi_{c}(\omega)\,\chi_{c}^{*}(-\omega)\,\times\\\
&\left[-\tilde{\Delta}_{c}\\!\sqrt{\kappa_{c}}X_{c}^{\rm
in}(\omega)+\left(\frac{\kappa_{c}}{2}{-}i\omega\right)\left(\\!\\!\sqrt{\kappa_{c}}Y_{c}^{\rm
in}(\omega)+G_{c}\chi_{b}(\omega)\Xi_{\rm
back}(\omega)\right)\right],\end{split}$ (12)

where $\Xi_{\rm back}(\omega)$ denotes the optomechanical backaction noise
from the mechanical oscillator, given by

$\Xi_{\rm
back}(\omega)=\frac{\xi(\omega)+\\!\sqrt{\kappa_{c}}\,\chi_{c}(\omega)\chi_{c}^{*}(-\omega)\,{\cal
C}(\omega)-\\!\sqrt{\kappa_{m}}\,\chi_{m}(\omega)\chi_{m}^{*}(-\omega){\cal
M}(\omega)}{1-\chi_{b}(\omega)\left(G_{c}^{2}\,\tilde{\Delta}_{c}\,\chi_{c}(\omega)\chi_{c}^{*}(-\omega)+\left|G_{m}\right|^{2}\tilde{\Delta}_{m}\,\chi_{m}(\omega)\chi_{m}^{*}(-\omega)\right)},$
(13)

which further contains three noise sources: $i$) the noise directly from the
mechanical thermal bath $\xi(\omega)$, $ii$) the optomechanical backaction
noise from the cavity field ${\cal C}(\omega)$, and $iii$) the magnomechanical
backaction noise from the magnon mode ${\cal M}(\omega)$, with

$\begin{split}{\cal
C}(\omega):=\,&G_{c}\left[\left(\frac{\kappa_{c}}{2}-i\omega\right)X_{c}^{\rm
in}(\omega)+\tilde{\Delta}_{c}Y_{c}^{\rm in}(\omega)\right],\\\ {\cal
M}(\omega):=\,&{\rm
Re}G_{m}\left[\left(\frac{\kappa_{m}}{2}-i\omega\right)X_{m}^{\rm
in}(\omega)+\tilde{\Delta}_{m}Y_{m}^{\rm in}(\omega)\right]\\\ &+{\rm
Im}G_{m}\left[-\tilde{\Delta}_{m}X_{m}^{\rm
in}(\omega)+\left(\frac{\kappa_{m}}{2}-i\omega\right)Y_{m}^{\rm
in}(\omega)\right].\end{split}$ (14)

We have introduced the natural susceptibility $\chi_{b}(\omega)$ of the
mechanical mode, $\chi_{c}(\omega)$ of the cavity field, and
$\chi_{m}(\omega)$ of the magnon mode, given by

$\begin{split}\chi_{b}(\omega)&=\frac{\omega_{b}}{\omega_{b}^{2}-\omega^{2}-i\gamma_{b}\omega},\\\
\chi_{c}(\omega)&=\frac{1}{\frac{\kappa_{c}}{2}+i(\tilde{\Delta}_{c}-\omega)},\\\
\chi_{m}(\omega)&=\frac{1}{\frac{\kappa_{m}}{2}+i(\tilde{\Delta}_{m}-\omega)}.\end{split}$
(15)

Figure 4: NSD of the phase $S_{Y_{c}}(\omega)$ for $N_{m}=10^{10}$ and $T=293$
K. Other parameters are listed in the text.

Figure 5: (a) Standard deviation $\sqrt{\left\langle\delta
Y_{c}^{2}\right\rangle}$ and (b) SNR (in decibels) versus magnon population
$\left|\left\langle m\right\rangle\right|^{2}$ at different temperatures
$T=293$ K (red solid), $T=4$ K (blue solid), and $T=10$ mK (green solid). The
black dashed line denotes vacuum fluctuations. The other parameters are those
as in Fig. 4, that is, the average $\left|\left\langle
Y_{c}\right\rangle\right|$ used to calculate SNR is the same as in Fig. 2(a).

From Eq. (12), we can define the noise spectral density (NSD) of the phase
quadrature

$\begin{split}S_{Y_{c}}(\omega)=\frac{1}{4\pi}&\int_{-\infty}^{+\infty}\\!d\omega^{\prime}\,e^{-i(\omega+\omega^{\prime})t}\,\,\times\\\
&\left\langle\delta Y_{c}(\omega)\delta Y_{c}(\omega^{\prime})+\delta
Y_{c}(\omega^{\prime})\delta Y_{c}(\omega)\right\rangle.\end{split}$ (16)

Here $S_{Y_{c}}(\omega)$ can be achieved using the input noise correlations in
the frequency domain, which is shown in Fig. 4. Clearly, a resonant peak
appears at the mechanical frequency $\omega=\omega_{b}$. The parameters are
those as in Fig. 2 and the others are $\omega_{m}/2\pi=10$ GHz,
$\kappa_{m}/2\pi=1$ MHz, $\tilde{\Delta}_{m}=0.5\kappa_{m}$,
$\gamma_{b}/2\pi=100$ Hz, $N_{m}=10^{10}$, and $T=293$ K.

The variance of the phase can be obtained by integrating $S_{Y_{c}}(\omega)$
over the frequency in the whole regime, i.e.,

$\left\langle\delta
Y_{c}(t)^{2}\right\rangle=\frac{1}{2\pi}\\!\int_{-\infty}^{+\infty}\\!\\!d\omega\,S_{Y_{c}}(\omega),$
(17)

which reflects the noise/imprecision of the phase. In our notation,
$\left\langle\delta Y_{c}(t)^{2}\right\rangle=\frac{1}{2}$ denotes vacuum
fluctuations, and the corresponding standard deviation
$\sqrt{\left\langle\delta Y_{c}(t)^{2}\right\rangle}\simeq 0.707$.

The classical average and the standard deviation of the phase allow us to
define the SNR

$\textup{SNR}=\frac{\left|\left\langle
Y_{c}\right\rangle\right|}{\sqrt{\left\langle\delta Y_{c}^{2}\right\rangle}},$
(18)

which is a key parameter and determines the resolution of magnons in our
approach. The resolution can be defined as the value of the magnon population
that can be detected with a unit SNR. It can be improved by enlarging the
signal $\left|\left\langle Y_{c}\right\rangle\right|$, i.e., increasing
$g_{mb}/\\!\sqrt{\kappa_{c}\omega_{b}}$ from Eq. (8), or reducing the noise
$\left\langle\delta Y_{c}^{2}\right\rangle$. The latter can be realized by
placing the system at cryogenic temperatures, or injecting a squeezed light to
reduce the phase noise. In Fig. 5(a), we show the phase noise (standard
deviation) at $T=10$ mK, 4 K, and 293 K. Obviously, the phase noise is reduced
by lowering the environmental temperature, which accordingly results in an
improved SNR as displayed in Fig. 5(b). With the parameters employed in Fig.
5, the SNR ranges from 24 to 35 dB at $T\leq 4$ K (from 16 to 26 dB at room
temperature) for $N_{m}\in[1.7\times 10^{9},1.5\times 10^{10}]$. The
resolution of magnons is in the order of magnitude of $10^{6}$. The SNR can be
further improved by using a squeezed light to break the standard quantum
limit. A 15 dB squeezing 15dB  can improve the SNR by 7.5 dB, leading to the
resolution in the order of magnitude of $10^{5}$. A potentially significant
improvement is possible if design the system to have a much larger value of
$g_{mb}/\\!\sqrt{\kappa_{c}\omega_{b}}$, such as by increasing the cavity
length $L$, since $\kappa_{c}\propto L^{-1}$, and reducing the sphere radius
$R$ if using a YIG sphere, since $g_{mb}\propto R^{-2}$ and $\omega_{b}\propto
R^{-1}$ Tang .

Lastly, it would be beneficial to compare our optical approach with the
approaches using a superconducting qubit Nak17 ; Nak20L . The latter
approaches can achieve very high resolutions of magnons. However,
superconducting qubits require ultra-cold temperatures (tens of mK) to
operate. Instead, our approach is hard to achieve as high resolution as in
Refs. Nak17 ; Nak20L  (however, depending on the system parameters), but our
approach can work at room temperature and with a very high SNR. This would be
highly appreciated for room-temperature experiments in cavity magnonics.
Moreover, the methods of Refs. Nak17 ; Nak20L  typically measure a very low
magnon population. Our method can, however, measure a much wider range of the
magnon population.

##  IV Conclusion

We propose an optical approach for measuring the steady-state magnon
population in a ferromagnet or ferrimagnet. It utilizes the magnetoelasticity
of the ferromagnet and the optomechanical coupling between the deformation
displacement and an optical cavity. The linear dependence between the optical
phase and the magnon population (under appropriate conditions) makes our
protocol a good meter for measuring the steady-state magnon population. The
study of the phase noise confirms that the protocol still works in the
presence of thermal noises even at room temperature, reflected in its high
SNR. We also provide strategies on how to improve the SNR and the resolution
of magnon excitation numbers. We expect that our work can offer an alternative
method for measuring the magnon population, as a complement to the existing
approaches prb21 ; Melkov ; invSH1 ; invSH2 ; invSH3 ; BLS1 ; BLS2 ; Nak17 ;
Nak20L .

## Acknowledgments

This work has been supported by Zhejiang Province Program for Science and
Technology (Grant No. 2020C01019), the National Natural Science Foundation of
China (Grants Nos. U1801661, 11874249, 11934010, 12174329), and the
Fundamental Research Funds for the Central Universities (No. 2021FZZX001-02).

## Appendix

Here we provide the details on the quantization of the magnetization, the
elastic strain, and the corresponding magnetoelastic energy, which lead to the
Hamiltonian we use and a dominant dispersive coupling between magnons and
phonons under appropriate assumptions.

The interaction between the magnetization and the elastic strain is described
by the magnetoelastic coupling Becker ; Kittle . In general, the
magnetoelastic energy density is given by

$\begin{split}f_{me}=&\
\frac{b_{1}}{M_{S}^{2}}\left(M_{x}^{2}\epsilon_{xx}+M_{y}^{2}\epsilon_{yy}+M_{z}^{2}\epsilon_{zz}\right)\\\
&+\frac{2b_{2}}{M_{S}^{2}}\left(M_{x}M_{y}\epsilon_{xy}+M_{x}M_{z}\epsilon_{xz}+M_{y}M_{z}\epsilon_{yz}\right),\end{split}$
(A1)

where $b_{1}$ and $b_{2}$ are the magnetoelastic coupling coefficients (not
specified), $M_{S}$ is the saturation magnetization, and $M_{x,y,z}$ are the
corresponding magnetization components. The magnetoelastic energy density is
related to the strain tensor $\epsilon_{ij}$,
$\epsilon_{ij}=\frac{1}{2}\left(\frac{\partial u_{i}}{\partial
l_{j}}+\frac{\partial u_{j}}{\partial l_{i}}\right)$, where $u_{i}$ are the
components of the displacement vector $\vec{u}$.

The magnetization can be quantized (i.e. magnons) as

$\hat{m}=\sqrt{\frac{V}{2\hbar\gamma M_{S}}}\left(M_{x}-iM_{y}\right),$ (A2)

where $V$ is the volume of the ferromagnet. We therefore obtain

$\begin{split}M_{x}=&\ \sqrt{\frac{\hbar\gamma
M_{S}}{2V}}\left(\hat{m}+\hat{m}^{\dagger}\right)\\\ M_{y}=&\
i\sqrt{\frac{\hbar\gamma
M_{S}}{2V}}\left(\hat{m}-\hat{m}^{\dagger}\right),\end{split}$ (A3)

and consequently,

$M_{z}=\left(M_{S}^{2}-M_{x}^{2}-M_{y}^{2}\right)^{\frac{1}{2}}\approx
M_{S}-\frac{\hbar\gamma}{V}\hat{m}^{\dagger}\hat{m}.$ (A4)

By replacing $M_{x,y,z}$ with Eqs. (A3) and (A4) in the magnetoelastic energy
density $f_{me}$, and integrating over the whole volume of the ferromagnet,
the first term in Eq. (A1) yields the Hamiltonian (neglecting nonresonant
fast-oscillating terms)

$H_{1}\simeq\frac{b_{1}}{M_{S}}\frac{\hbar\gamma}{V}\hat{m}^{\dagger}\hat{m}\int
dl^{3}\left(\epsilon_{xx}+\epsilon_{yy}-2\epsilon_{zz}\right),$ (A5)

which accounts for the dispersive interaction between magnons and phonons, and
the second term leads to the Hamiltonian

$\begin{split}H_{2}=&\,i\frac{b_{2}}{M_{S}}\frac{\hbar\gamma}{V}\left(\hat{m}^{2}-\hat{m}^{\dagger
2}\right)\int dl^{3}\epsilon_{xy}\\\
+&\frac{2b_{2}}{M_{S}^{2}}\sqrt{\frac{\hbar\gamma
M_{S}}{2V}}\left(M_{S}-\frac{\hbar\gamma}{V}\hat{m}^{\dagger}\hat{m}\right)\left[\hat{m}\int
dl^{3}\left(\epsilon_{xz}+i\epsilon_{yz}\right)+h.c.\right],\end{split}$ (A6)

which describes the parametric magnon generation when the phonon frequency is
twice the magnon frequency, or the linear magnon-phonon coupling when they are
(nearly) resonant.

The present work utilizes the dispersive magnon-phonon coupling, corresponding
to the situation where the phonon frequency is much lower than the magnon’s.
This typically occurs for a large-sized ferromagnet Tang ; Davis .

The magnetoelastic displacement can be expressed as a superposition

$\vec{u}=\sum_{n,m,k}d^{\left(n,m,k\right)}\vec{\chi}^{(n,m,k)}\left(x,y,z\right),$
(A7)

with $\vec{\chi}^{(n,m,k)}\left(x,y,z\right)$ the (normalized) displacement
eigenmode and $d^{\left(n,m,k\right)}$ the corresponding amplitude. We can
then quantize the mechanical motion as

$\begin{split}d^{\left(n,m,k\right)}=d_{\text{zpm}}^{\left(n,m,k\right)}\left(\hat{a}_{n,m,k}+\hat{a}^{\dagger}_{n,m,k}\right),\end{split}$
(A8)

where $d_{\text{zpm}}^{\left(n,m,k\right)}$ is the amplitude of the zero-point
motion, and $\hat{a}_{n,m,k}$ and $\hat{a}^{\dagger}_{n,m,k}$ are boson
operators for each mode with specified mode indices $\left(n,m,k\right)$.

Substituting Eqs. (A7) and (A8) into the dispersive interaction Hamiltonian
$H_{1}$, we obtain

$H_{1}=\sum_{n,m,k}\hbar
g_{mb}^{\left(n,m,k\right)}\hat{m}^{\dagger}\hat{m}\left(\hat{a}_{n,m,k}+\hat{a}^{\dagger}_{n,m,k}\right),$
(A9)

where $g_{mb}^{\left(n,m,k\right)}$ is the magnon-phonon coupling strength,
given by

$g_{mb}^{\left(n,m,k\right)}=\frac{b_{1}}{M_{S}}\frac{\gamma}{V}\int
dl^{3}d_{\text{zpm}}^{\left(n,m,k\right)}\left(\frac{\partial\chi_{x}^{(n,m,k)}}{\partial
x}+\frac{\partial\chi_{y}^{(n,m,k)}}{\partial
y}-2\frac{\partial\chi_{z}^{(n,m,k)}}{\partial z}\right).$ (A10)

If considering a specific mechanical mode and its motion in only one
direction, we have the following Hamiltonian

$H_{1}^{\rm sing}=\hbar
g_{mb}\hat{m}^{\dagger}\hat{m}\left(\hat{a}+\hat{a}^{\dagger}\right),$ (A11)

that is the one we use in the Hamiltonian Eq. (2) for the dispersive
magnomechanical interaction (since
$\hat{q}=\frac{1}{\sqrt{2}}(\hat{a}+\hat{a}^{\dagger})$). For simplicity, we
have omitted hat signs for the operators in the main text.

## References

  * (1) 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).
  * (2) Y. Tabuchi, S. Ishino, T. Ishikawa, R. Yamazaki, K. Usami, and Y. Nakamura, Phys. Rev. Lett. 113, 083603 (2014).
  * (3) X. Zhang, C. L. Zou, L. Jiang, and H. X. Tang, Phys. Rev. Lett. 113, 156401 (2014).
  * (4) D. Lachance-Quirion, Y. Tabuchi, A. Gloppe, K. Usami, and Y. Nakamura, Appl. Phys. Express 12, 070101 (2019).
  * (5) P. Pirro, V. I. Vasyuchka, A. A. Serga, and B. Hillebrands, Nat. Rev. Mater. 6, 1114 (2021).
  * (6) B. Z. Rameshti et al., arXiv:2106.09312.
  * (7) D. Lachance-Quirion, Y. Tabuchi, S. Ishino, A. Noguchi, T. Ishikawa, R. Yamazaki and Y. Nakamura, Sci. Adv. 3, e1603150 (2017).
  * (8) D. Lachance-Quirion, S. P. Wolski, Y. Tabuchi, S. Kono, K. Usami, and Y. Nakamura, Science 367, 425 (2020).
  * (9) J. Li, S.-Y. Zhu, and G. S. Agarwal, Phys. Rev. Lett. 121, 203601 (2018).
  * (10) J. Li, S.-Y. Zhu, and G. S. Agarwal, Phys. Rev. A 99, 021801(R) (2019).
  * (11) T. Hioki, H. Shimizu, T. Makiuchi, and E. Saitoh, Phys. Rev. B 104, L100419 (2021).
  * (12) A. G. Gurevich and G. A. Melkov, Magnetization Oscillations and Waves (CRC press, Boca Raton, 1996).
  * (13) E. Saitoh, M. Ueda, H. Miyajima, and G. Tatara, Appl. Phys. Lett. 88, 182509 (2006).
  * (14) Y. Kajiwara et al., Nature (London) 464, 262 (2010).
  * (15) A. V. Chumak, A. A. Serga, M. B. Jungfleisch, R. Neb, D. A. Bozhko, V. S. Tiberkevich, and B. Hillebrands, Appl. Phys. Lett. 100, 082405 (2012).
  * (16) S. O. Demokritov, B. Hillebrands, and A. N. Slavin, Phys. Rep. 348, 441 (2001).
  * (17) T. Sebastian, K. Schultheiss, B. Obry, B. Hillebrands, and H. Schultheiss, Front. Phys. 3, 35 (2015).
  * (18) 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).
  * (19) Y.-P. Wang, G.-Q. Zhang, D. Zhang, T.-F. Li, C.-M. Hu, and J. Q. You, Phys. Rev. Lett. 120, 057202 (2018).
  * (20) R.-C. Shen, Y.-P. Wang, J. Li, S.-Y. Zhu, G. S. Agarwal, and J. Q. You, Phys. Rev. Lett. 127, 183202 (2021).
  * (21) S. P. Wolski, D. Lachance-Quirion, Y. Tabuchi, S. Kono, A. Noguchi, K. Usami, and Y. Nakamura, Phys. Rev. Lett. 125, 117701 (2020).
  * (22) C. Kittel, Phys. Rev. 110, 836 (1958).
  * (23) X. Zhang, C.-L. Zou, L. Jiang, and H. X. Tang, Sci. Adv. 2, e1501286 (2016).
  * (24) C. A. Potts, E. Varga, V. Bittencourt, S. V. Kusminskiy, and J. P. Davis, Phys. Rev. X 11, 031053 (2021).
  * (25) M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, Rev. Mod. Phys. 86, 1391 (2014).
  * (26) The model is, however, not limited to a spherical structure and YIG, but valid for any ferrimagnet or ferromagnet with a dispersive magnon-phonon coupling and a detectable deformation displacement.
  * (27) C. Kittel, Phys. Rev. 73, 155 (1948).
  * (28) A. Osada et al., Phys. Rev. Lett. 116, 223601 (2016).
  * (29) J. Li and S.-Y. Zhu, New J. Phys. 21, 085001 (2019).
  * (30) H. Tan, Phys. Rev. Res. 1, 033161 (2019).
  * (31) C. Kong, B. Wang, Z.-X. Liu, H. Xiong, and Y. Wu, Opt. Express 27, 5544 (2019).
  * (32) M.-S. Ding, L. Zheng, and C. Li, J. Opt. Soc. Am. B 37, 627 (2020).
  * (33) M. Yu, H. Shen, and J. Li, Phys. Rev. Lett. 124, 213604 (2020).
  * (34) C. A. Potts, V. A. S. V. Bittencourt, S. V. Kusminskiy, and J. P. Davis, Phys. Rev. Applied 13, 064001 (2020).
  * (35) Z.-B. Yang, J.-S. Liu, A.-D. Zhu, H.-Y. Liu, and R.-C. Yang, Ann. Phys. (Amsterdam) 532, 2000196 (2020).
  * (36) J. Li and S. Gröblacher, Quantum Sci. Technol. 6, 024005 (2021).
  * (37) J. Li, Y.-P. Wang, J. Q. You, and S.-Y. Zhu, arXiv:2101.02796. Nat. Sci. Rev. (in press).
  * (38) W. Zhang, D.-Y. Wang, C.-H. Bai, T. Wang, S. Zhang, and H.-F. Wang, Opt. Express 29, 11773 (2021).
  * (39) B. Sarma, T. Busch, and J. Twamley, New J. Phys. 23, 043041 (2021).
  * (40) S.-F. Qi and J. Jing, Phys. Rev. A 103, 043704 (2021).
  * (41) Y.-T. Chen, L. Du, Y. Zhang, and J.-H. Wu, Phys. Rev. A 103, 053712 (2021).
  * (42) T.-X. Lu, H. Zhang, Q. Zhang, and H. Jing, Phys. Rev. A 103, 063708 (2021).
  * (43) M.-S. Ding, X.-X. Xin, S.-Y. Qin, and C. Li, Opt. Commun. 490, 126903 (2021).
  * (44) S. Gröblacher, K. Hammerer, M. R. Vanner, and M. Aspelmeyer, Nature 460, 724 (2009).
  * (45) C. W. Gardiner and M. J. Collett, Phys. Rev. A 31, 3761 (1985).
  * (46) R. Benguria and M. Kac, Phys. Rev. Lett. 46, 1 (1981); V. Giovannetti and D.Vitali, Phys. Rev.A 63, 023812 (2001).
  * (47) M. Aspelmeyer, S. Gröblacher, K. Hammerer, and N. Kiesel, J. Opt. Soc. Am. B 27, A189 (2010).
  * (48) P. Grangier, J. A. Levenson, and J.-P. Poizat, Nature 396, 537 (1998).
  * (49) I. S. Gradshteyn and I. M. Ryzhik, Table of Integrals, Series and Products (Academic, Orlando, 1980), Page 1119.
  * (50) H. Vahlbruch, M. Mehmet, K. Danzmann, and R. Schnabel, Phys. Rev. Lett. 117, 110801 (2016).
  * (51) R. Becker and W. Döring, Ferromagnetismus (Verlag Julius Springer, Berlin, 1939).


# A note on heterogeneity, trade integration and spatial inequality††thanks:
This work is dedicated to Carlos Hervés-Beloso. I am deeply indebted to Sofia
B. S. D. Castro and João Correia da Silva for their invaluable comments and
suggestions and also for their help in revising previous and related versions
of this work. I am also very thankful to the audience at the "Porto Workshop
on Economic Theory and Applications - in honor of Carlos Hervés-Beloso" for
the fruitful discussions. This research has been financed by Portuguese public
funds through FCT - Fundação para a Ciência e a Tecnologia, I.P., in the
framework of the projects with references UIDB/04105/2020 and UIDB/00731/2020.
Part of this research was developed while José M. Gaspar was a researcher at
the Research Centre in Management and Economics, Católica Porto Business
School, Universidade Católica Portuguesa, through the grant
CEECIND/02741/2017.

José M. Gaspar School of Economics and Management and CEF.UP, University of
Porto. Email<EMAIL_ADDRESS>
###### Abstract

We study the impact of economic integration on spatial development in a model
where all consumers are inter-regionally mobile and have heterogeneous
preferences regarding their residential location choices. This heterogeneity
is the unique dispersion force in the model. We show that, under reasonable
values for the elasticity of substitution among varieties of consumption
goods, a higher trade integration always promotes more symmetric patterns,
irrespective of the functional form of the dispersion force. We also show that
an increase in the degree of heterogeneity in preferences for location leads
to less spatial inequality.

Keywords: heterogeneous location preferences; economic integration; economic
geography; spatial development.

JEL codes: R10, R12, R23

##  1 Introduction

In this paper, we study the possible spatial distributions in the $2$-region
Core-Periphery (CP) model by Murata (2003), where all labour is free to
migrate between regions and workers are heterogeneous regarding their
preferences for location. Our aim is to investigate how spatial patterns
evolve as regions become more integrated and whether the degree of
heterogeneity has qualitative impact on this relationship. Heterogeneity in
preferences for location generates a local dispersion force (akin to
congestion), which we model as a utility penalty, following the recent work of
Castro et al. (2022). This framework can been shown to encompass various
functional forms used in discrete choice models of probabilistic migration,
including the well-known Logit model (used e.g. by Murata, 2003; Tabuchi and
Thisse, 2002).

We show that, irrespective of the degree of heterogeneity, a higher trade
integration always promotes more symmetric spatial patterns. Since
heterogeneity does not co-vary with transportation costs, the latter only
affect the strength of agglomeration forces due to increasing returns. This is
a reversion of the usual prediction that lower transport costs lead to
agglomeration, and contrasts the findings that consumer heterogeneity leads to
a bell-shaped relationship between decreasing transport costs and spatial
inequality (Murata, 2003; Tabuchi and Thisse, 2002). Since there is no
immobile workforce (no fixed regional internal demand), there is a lower
incentive for firms to relocate to smaller markets in order to capture a
higher share of local demand and avoid competition when transport costs are
higher. Krugman and Elizondo (1996), Helpman (1998) and Murata and Thisse
(2005), have also obtained similar results concerning the relationship between
transportation costs and spatial inequality. The first includes a congestion
cost in the core, the second considers a non-tradable housing sector and the
treatment is only numerical, while the latter’s prediction cannot be
disassociated from the interplay between inter-regional transportation costs
and intra-regional commuting costs. Recently, Tabuchi et al. (2018) have
reached similar conclusions to ours arguing that falling transport costs
increase the incentives for firms in peripheral regions to increase production
since they have a better access to the core, thus contributing to the
dispersion of economic activities. Allen and Arkolakis (2014) estimated that
the removal of the US Interstate Highway System, which, by limiting accesses,
can be interpreted as an increase in transportation costs, would cause a
redistribution of the population from more economically remote regions to less
remote regions in the US. This adds empirical validity to our findings.

The rest of the paper is organized as follows. In section 2 we describe the
model and the short-run equilibrium. In section 3 we study the qualitative
properties of the long-run equilibria. Finally, section 4 is left for some
concluding remarks.

##  2 The model

There are two regions $L$ and $R$ with a total population of mass $1$.
Consumers are assumed to have heterogeneous preferences with respect to the
region in which they reside. We follow Castro et al. (2020) and assume that
these preferences are described by a parameter $x,$ uniformly distributed
along the interval $[0,1]$. The consumer with preference described by $x=0$
(resp. $x=1$) has the highest preference for residing in region $L$ (resp.
region $R$). An agent whose preference corresponds to $x=1/2$ is indifferent
between either region. We can thus refer to each consumer with respect to his
preference towards living in region $L$ as $x\in[0,1]$.

For a consumer with preference $x$, the utilities from living in region $L$
and $R$ are given, respectively, by:

$\displaystyle U_{L}(x)=$ $\displaystyle
U\left(u\left(C_{L}\right),t(x)\right)$ $\displaystyle U_{R}(x)=$
$\displaystyle U\left(u\left(C_{R}\right),t(1-x)\right),$ (1)

where $C_{i}$ denotes the level of consumption of a consumer living in region
$i\in\left\\{L,R\right\\}$, and $t(x)$ is the utility penalty associated with
living in $L$, while $t(1-x)$ is the utility penalty a consumer faces from
living in region $R.$ Regions are symmetric in all respects. We assume that
$t^{\prime}(.)>0$. There are two regions $L$ and $R$ with a total population
of mass $1$.

The short-run equilibrium is equivalent to Murata (2003) and Tabuchi et al.
(2018), which we briefly describe here. The consumption aggregate is a CES
composite given by:

$C=\left[\int_{s\in
S}c(s)^{\tfrac{\sigma-1}{\sigma}}ds\right]^{\tfrac{\sigma}{\sigma-1}},$

where $s$ stands for the variety produced by each monopolistically competitive
firm and $\sigma$ is the constant elasticity of substitution between any two
varieties. The consumer is subject to the budget constraint
$P_{i}C_{i}=w_{i}$, where $P_{i}$ is the price index and $w_{i}$ is the
nominal wage in region $i$. Utility maximization by a consumer in region $i$
yields the following demand for each manufactured variety produced in $j$ and
consumed in $i$:

$c_{ij}=\dfrac{p_{ij}^{-\sigma}}{P_{i}^{1-\sigma}}w_{i},$ (2)

where:

$P_{i}=\left[\int_{s\in
S}p_{i}(s)^{1-\sigma}d(s)\right]^{\tfrac{1}{1-\sigma}},$ (3)

is the manufacturing price index for region $i$. A manufacturing firm faces
the following cost function:

$TC(q)=w\left(\beta q+\alpha\right),$

where $q$ corresponds to a firm’s total production of manufacturing goods,
$\beta$ is the input requirement (per unit of output) and $\alpha$ is the
fixed input requirement. The manufacturing good is subject to trade barriers
in the form of iceberg costs, $\tau\in(1,+\infty)$: a firm ships $\tau$ units
of a good to a foreign region for each unit that arrives there. Assuming free
entry in the manufacturing sector, at equilibrium firms will earn zero
profits, which translates into the following condition:

$\displaystyle\pi\equiv\left(p-\beta w\right)q-\alpha w$ $\displaystyle=0,$

which gives the firm’s total equilibrium output, symmetric across regions:

$q_{L}=q_{R}=\frac{\alpha(\sigma-1)}{\beta},$ (4)

and the following profit maximizing prices:

$p_{LL}=\dfrac{\beta\sigma}{\sigma-1}w_{L}\text{ \ and\
}p_{RL}=\dfrac{\beta\sigma\tau}{\sigma-1}w_{L}.$ (5)

Labour-market clearing implies that the number of agents in region $L$ equals
labour employed by a firm times the number of varieties produced in a region:

$h=n_{L}\left(\alpha+\beta q_{L}\right),$

from where, using (4), we get:

$n_{L}=\dfrac{h}{\sigma\alpha}\text{ \ and \
}n_{R}=\dfrac{1-h}{\sigma\alpha}.$ (6)

Choosing labour in region $R$ as the numeraire, we can normalize $w_{R}$ to
$1$. The price indices in $L$ and $R$ using (3) are given, after (6), by:

$\displaystyle P_{L}$
$\displaystyle=\left[\dfrac{h}{\sigma\alpha}\left(\dfrac{\beta\sigma}{\sigma-1}w\right)^{1-\sigma}+\dfrac{1-h}{\sigma\alpha}\left(\dfrac{\beta\sigma\tau}{\sigma-1}\right)^{1-\sigma}\right]^{\tfrac{1}{1-\sigma}},$
$\displaystyle P_{R}$
=$\left[\dfrac{1-h}{\sigma\alpha}\left(\dfrac{\beta\sigma}{\sigma-1}\right)^{1-\sigma}+\dfrac{h}{\sigma\alpha}\left(\dfrac{\beta\sigma\tau}{\sigma-1}w\right)^{1-\sigma}\right]^{\tfrac{1}{1-\sigma}}$.
(7)

Rewriting a firm’s profit in region $L$ as:

$\pi_{L}=\left(p_{LL}-\beta w\right)\left[hc_{LL}+(1-h)\tau
c_{RL}\right]-\alpha w,$

the zero profit condition yields, after using (2), (5) and (7), the following
wage equation:

$h=\dfrac{w^{\sigma}-\phi}{w^{\sigma}-\phi+w(w^{-\sigma}-\phi)},$ (8)

where $\phi\equiv\tau^{1-\sigma}\in(0,1)$ is the _freeness of trade_. The
nominal wage $w$ can be implicitly defined as a function of the spatial
distribution in region $L$, $h\in[0,1]$.111The conditions for application of
the Implicit Function Theorem are shown to be satisfied in Appendix A. We have
$w(0)=\phi^{1/\sigma}<1$, $w(1/2)=1$ and $w(1)=\phi^{-1/\sigma}>1$. We say
that $L$ is larger than $R$ when $h>1/2,$ and smaller otherwise. Figure 1
illustrates $w(h)$ for $h\in[0,1]$ with parameter values $\sigma=2$ and
$\phi=0.5$.222Our choice of parameters is made for graphical convenience and
extends to the whole parameter range, unless stated otherwise. Therefore, the
conclusions under our choice are consistent with more realistic empirical
estimates.

Figure 1: Short-run equilibrium relative wage as a function of the consumers
in region $L$. We set $\sigma=2$ and $\phi=0.5$.

The nominal wage is strictly increasing in the spatial distribution (see
Appendix A). Moreover, see Appendix A, when $h>1/2$, the relative wage $w$
increases when trade barriers are higher. Conversely, if $h<1/2$, higher trade
barriers implies lower wages in $L$. Therefore higher trade barriers increase
the wage divergence between the regions.

The intuition is as follows. When all workers are mobile, they can move to the
region that offers them a relatively better access to varieties. This
advantage of the larger region becomes higher as transport costs increase
because markets become more focused on local demand. This increases
expenditures in the larger region relative to the smaller, which in turn
pushes the nominal relative wage upwards.

We consider a general isoelastic sub-utility for consumption goods:

$\begin{cases}u_{i}=\dfrac{C_{i}^{1-\theta}-1}{1-\theta},&\text{ if
}\theta\in[0,1)\cup(1,+\infty)\\\ u_{i}=\ln C_{i},&\text{ if
}\theta=1,\end{cases},$ (9)

where for $\theta=1$ we take the limit value of the upper expression of
$u_{i}$. The parameter $\theta$ is a positive agglomeration externality. It
influences how a change in the regional consumption differential,
$C_{L}-C_{R}$, impacts the regional utility differential, $u_{L}-u_{R}$, i.e.,
it influences the strength of the self-reinforcing agglomeration mechanism
when one region is more populated than the other. We shall detail further
below why this specification is useful.

The utility differential from consumption goods, $\Delta u\equiv u_{L}-u_{R},$
is given by:

$\Delta
u=\begin{cases}\dfrac{\eta}{1-\theta}\left\\{w^{1-\theta}\left[(1-h)\phi+hw^{1-\sigma}\right]^{\tfrac{1-\theta}{\sigma-1}}-\left[1-h+h\phi
w^{1-\sigma}\right]^{\tfrac{1-\theta}{\sigma-1}}\right\\},&\ \text{if
}\mbox{$\theta\neq$1}\\\ \ln
w+\dfrac{1}{\sigma-1}\ln\left[\frac{hw^{1-\sigma}+(1-h)\phi}{h\phi
w^{1-\sigma}+(1-h)}\right],&\ \text{if $\theta$=1.}\end{cases}$ (10)

We adopt the normalizations by Fujita _et al_. (1999), i.e., $\alpha\sigma=1$
so that the number of consumers in a region equals its number of firms.
Moreover, we assume that $(\sigma-1)/(\sigma\beta)=1$ so that the price of
each manufactured variety in a region equals its workers’ nominal wage. These
imply that $\eta=1$.

Figure 2 depicts $\Delta u$ in (10) for different levels of $\theta$.

Figure 2: Utility differential $\Delta u=u_{L}-u_{R}$ for different levels of
$\theta$. Parameter values are $\sigma=2$ and $\phi=0.5$.

We observe that a higher $\theta$ increases the utility differential $\Delta
u$ for any spatial distribution $h>1/2$. Therefore, if $L$ is relatively more
industrialized, a higher $\theta$ increases the attractiveness of $L$ relative
to $R$ for consumers. Thus, it strengthens the agglomeration forces towards
region $L$.

The purpose of the general specification in (9) is explained as follows.
Firstly, it allows us to derive general stability conditions for long-run
equilibria for any $t(x)$. Secondly, it encompasses the Murata (2003) model as
a particular case when $\theta=0$, which is of great interest for comparison
purposes.

##  3 Long-run equilibria

In a long-run spatial equilibrium, each worker chooses to live in the region
that provides a higher utility. Hence, as in Castro et al. (2022), a long-run
equilibrium must be such that workers with $x\in[0,h)$ live in region $L$ and
workers with $x\in(h,1]$ live in region $R$. Knowing from utility maximization
that $C=W/P$ and assuming that $t(x)$ enters additively in overall utility
$U_{i}(x)$, we rewrite the indirect utilities as:

$\displaystyle V_{L}(h)$
$\displaystyle=\dfrac{C_{L}^{1-\theta}-1}{1-\theta}-t(h)$ $\displaystyle
V_{R}(h)$ $\displaystyle=\dfrac{C_{R}^{1-\theta}-1}{1-\theta}-t(1-h),$ (11)

where $h$ satisfies the short-run equilibrium in (8).

###  3.1 Interior equilibria

Due to symmetry, we focus only on the case where $L$ is larger or the same
size as $R$, i.e., $h\geq 1/2$. We define an interior equilibrium
$h^{*}\in[1/2,1)$ as a spatial distribution $h$ that satisfies both (8) and
$V_{L}=V_{R}$. Such an equilibrium is said to be stable if
$d\left(V_{L}-V_{R}\right)/dh<0$. This leads to the following Lemma.

######  Lemma 1.

An interior equilibrium $h^{*}\in[1/2,1)$ is stable if:

$\zeta\left\\{\varphi\left[\dfrac{w^{\sigma}+w^{1-\sigma}-(w+1)\phi}{w\left(1-\phi^{2}\right)}\right]^{\frac{\theta-1}{\sigma-1}}-\psi\left[\dfrac{w^{2\sigma}-(w+1)\phi
w^{\sigma}+w}{w\left(1-\phi^{2}\right)}\right]^{\frac{\theta+\sigma-2}{\sigma-1}}\right\\}<t^{\prime}\left(h^{*}\right)+t^{\prime}\left(1-h^{*}\right),$

where $h^{*}$ satisfies the short-run equilibrium condition in (8), and:

$\displaystyle\zeta$
$\displaystyle=\dfrac{w^{-(\theta+\sigma+1)}}{(\sigma-1)\left[(\sigma-1)\phi+(\sigma-1)\phi
w^{2\sigma}+\left(-2\sigma+\phi^{2}+1\right)w^{\sigma}\right]};$
$\displaystyle\varphi$
$\displaystyle=w\left\\{\phi\left[\sigma+(\sigma-1)w\right]w^{\sigma}-2\sigma
w+w\right\\}\left[w^{2\sigma}-(w+1)\phi w^{\sigma}+w\right];$
$\displaystyle\psi$
$\displaystyle=\left(1-\phi^{2}\right)w^{\theta+\sigma+1}\left[(2\sigma-1)w^{\sigma}-\sigma(w+1)\phi+\phi\right].$

###### Proof.

See Appendix B. ∎

An interior equilibrium is called _symmetric dispersion_ if $h^{*}=1/2$ and
_partial agglomeration_ otherwise. While the former always exists, the latter
depends on the form of $t(x)$. Knowing that $w=1$ for $h=1/2$ the stability
condition in Lemma 1 for symmetric dispersion simplifies to:

$\dfrac{d\Delta
u}{dh}\left(\tfrac{1}{2}\right)\equiv\frac{2(2\sigma-1)(1-\phi)\left(\frac{1+\phi}{2}\right)^{\frac{1-\theta}{\sigma-1}}}{(\sigma-1)(2\sigma+\phi-1)}<t^{\prime}\left(\tfrac{1}{2}\right).$
(12)

A careful inspection of (12) allows us to conclude that the LHS is decreasing
in $\phi$ if $\theta\geq 1$. For $\theta<1$, it is decreasing in $\phi$ if
$\sigma\geq 1+\sqrt{2}/2\approx 1.71$. Symmetric dispersion then becomes
easier to sustain under lower transportation costs if $\sigma>1+\sqrt{2}/2$,
which, according to recent empirical estimations for $\sigma$, is more than
reasonable.333Estimations evidence that $\sigma$ should be significantly
larger than unity (Crozet, 2004; Head and Mayer, 2004; Niebuhr, 2006; Bosker
_et al.,_ 2010). Anderson and Wincoop (2004), for instance, find that it is
likely to range from 5 to 10. Since $t(h)$ does not depend on $\phi$, this
result does _not_ depend on consumer heterogeneity. Since, at symmetric
dispersion, consumers in each region have access to the same amount of
manufactures, an exogenous migration will induce a lower (higher) benefit from
local consumption goods in the larger market if transport costs are lower
(higher). This is captured by the fact that the relative decrease in prices
and increase in wages (at the symmetric equilibrium) is more pronounced when
transport costs are higher.

Let us now assume that a partial agglomeration $h^{*}\in(1/2,1)$ equilibrium
exists and is unique. Regarding the influence of trade integration on partial
agglomeration equilibria, we establish the following result.

######  Proposition 1.

If a stable partial agglomeration $h^{*}\in(1/2,1)$ exists and is unique, it
becomes more symmetric as trade barriers decrease.

###### Proof.

See Appendix C. ∎

This proposition states that if most consumers reside in region $L$,
increasing the freeness of trade will lead to a smooth exodus from region $L$
to region $R$, irrespective of the value of $t(h^{*})$. More integration thus
increases the incentives for consumers to distribute more equally among the
two regions.

###  3.2 Logit preferences

For the two region case, and following Tabuchi and Thisse (2002) and Murata
(2003), the probability that a consumer will choose to reside in region $L$ is
given by the logit model written as:

$P_{L}(h)=\dfrac{e^{u_{L}(h)/\mu}}{e^{u_{L}(h)/\mu}+e^{u_{R}(h)/\mu}},$ (13)

where $\mu\geq 0$ is a scale parameter which measures the dispersion of
consumer preferences. If $\mu=0$, consumers do not care about their location
preferences but rather solely about relative wages. The number of agents $h$
is the same as the consumer $x$ who is indifferent between living in region
$L$ or in region $R$. Therefore, the probability $P_{L}(h)$ is tantamount to
the indifferent consumer $x=h$ (Castro et al. 2022). Thus, using (13), we can
write:

$h=\dfrac{e^{u_{L}(h)/\mu}}{e^{u_{L}(h)/\mu}+e^{u_{R}(h)/\mu}}.$ (14)

Manipulating (14) yields:

$u_{L}-\mu\ln h=u_{R}-\mu\ln(1-h),$ (15)

which rearranged yields the long-run equilibrium condition, $\Delta
u(h)=\Delta t(h),$ where $\Delta t(h)=\mu\left[\ln h-\ln(1-h)\right]$.

Notice from the (15) that $t(h)=\mu\ln h<0$ for $h\in[0,1]$, so that the
overall utility here is interpreted as the utility from consumption plus a
benefit from living in the most preferred region. The overall utility of a
consumer in region $i$ thus lies on the interval $[U_{i},+\infty)$. For a
strictly positive $\mu$, the consumer who likes region $R$ the most (at $x=1$)
will never want to live in $L$, because, at $h=1$, the overall utility gain
from moving from region $L$ to region $R$ is infinite. This means that the
Logit model penalizes (benefits) the consumers who are less (more) willing to
leave a region very strongly. Therefore, agglomeration in any region is not
possible.

This discussion allows us to formalize the following result, which summarizes
the possible spatial outcomes under Logit type preferences, depending on the
degree of heterogeneity $\mu.$

######  Proposition 2.

Under Logit type preferences, the spatial distribution depends on the level of
consumer heterogeneity as follows:

  * •

_Symmetric dispersion is the only stable equilibrium if:_

$\mu>\mu_{d}\equiv\dfrac{(2\sigma-1)(1-\phi)}{(\sigma-1)(2\sigma+\phi-1)};$
(16)

  * •

_Partial agglomeration $h^{*}\in\left(\tfrac{1}{2},1\right)$ is the only
stable equilibrium if $\mu<\mu_{d}$._

###### Proof.

See Appendix E. ∎

When consumer heterogeneity is low there is a single stable partial
agglomeration equilibrium that is very asymmetric (close to agglomeration). As
consumer heterogeneity increases, partial agglomeration corresponds to a more
even distribution. Finally, if consumer heterogeneity is high enough
($\mu>\mu_{d}$ in (16)), consumers disperse symmetrically across the two
regions.

In Figure 3 we show how the spatial distribution of industry changes as
regions become more integrated. We set $\sigma=2.5$ and $\mu=0.2$.

Figure 3: Utility differential $\Delta u$ (dashed lines) and penalty
differential $t(h)-t(1-h)$ (thick line) for increasing levels of $\phi$. For
$\phi=0.3$ (upper dashed line), dispersion is unstable and partial
agglomeration is stable. For$\phi=0.5$ (medium dashed line) dispersion is
unstable, partial agglomeration is stable and is less asymmetric. For
$\phi=0.9$ (lower dashed line) dispersion is stable and is the only
equilibrium. The values for the parameters are $\sigma=2.5$ and $\mu=0.2$.

For low levels of the freeness of trade, a unique stable partial agglomeration
equilibrium exists where most consumers reside in region $L$. This means that
there is just a small fraction of consumers in $R$ that are not willing to
forego their preferred region because the gain in consumption goods from doing
so is not high enough. As regions become more integrated, the home-market
effect becomes weaker and partial agglomeration becomes more symmetric.
Finally, for a high level of inter-regional integration, agglomeration forces
are so weak that no consumer is willing to leave his most preferred region and
only symmetric dispersion is stable.

The preceding analysis is summarized in the two bifurcation diagrams in Figure
4, along a smooth parameter path where $\mu$ and $\phi$ increase. We
numerically uncover two supercritical pitchfork bifurcations in $\mu$ and
$\phi$. The picture to the right in Figure 4 (increasing $\phi$) differs from
other supercritical pitchforks found in the literature (e.g., Pflüger, 2004;
Ikeda et al., 2022) crucially in the sense that the direction of change in
stability as $\phi$ increases is reversed, i.e., lower trade barriers leads to
more symmetric spatial distributions.

Figure 4: Bifurcation diagrams. To the left, the bifurcation parameter is
$\mu\in[0,1]$ with $\phi=0.4$. To the right, we have $\phi\in(0,1)$ as the
bifurcation parameter and set $\mu=0.2$. For both scenarios we use $\sigma=2$.

##  4 Discussion and concluding remarks

We have seen that heterogeneity in preferences for location alone bears no
impact on the relationship between trade integration and spatial inequality,
which is a monotonic decreasing one. This contrasts the findings in other
works with heterogeneity in consumer preferences, such as Tabuchi and Thisse
(2002) and Murata (2003), who show evidence of a bell-shaped relationship
between trade integration and spatial inequality. The former’s setting differs
from ours because the authors consider an inter-regionally immobile workforce
whose role as a dispersive force is enhanced by higher transportation costs.

However, it is particularly worthwhile to discuss the results of Murata
(2003). Murata (2003) found that the relationship between trade integration
and spatial inequality need not be monotonic and depends on the degree of
consumer heterogeneity, which is at odds with our findings. For instance, for
an intermediate degree of consumer heterogeneity, Murata finds that increasing
trade integration initially fosters agglomeration and later leads to re-
dispersion of industry. However, these conclusions can be shown to stem from
the author’s particular choice of the value for the elasticity of
substitution, $\sigma=1.25$. As we have argued in Section 3, such a low value
is empirically implausible. For exceedingly low values $(\sigma<1.71)$,
increasing returns at the firm level are too strong. Strong enough that the
utility gain at dispersion becomes increasing in $\phi$, instead of
decreasing. This would justify an initial concentration of industry as a
result of an increase in $\phi.$ However, we have seen that for a plausible
range of $\sigma$ the utility gain at dispersion always decreases with $\phi$.
Moreover, if $\theta>0$, this holds even for lower values of $\sigma$.444If
$\theta\geq 1$, the result holds for $\sigma>1$. Therefore, a higher $\phi$
always promotes more equitable distributions as opposed to asymmetric ones.

Hence, when workers are completely mobile, more trade integration ubiquitously
reduces the spatial inequalities between the two regions, irrespective of the
degree of heterogeneity in location preferences. Therefore, a _de facto_ lower
inter-regional labour mobility induced by consumer heterogeneity alone cannot
account for the predictions that a higher inter-regional integration will lead
to more unequal spatial development or an otherwise bell-shaped relationship
between the two.

By considering that all consumers are allowed to migrate if they so desire, we
have shown that a higher inter-regional trade integration always leads to less
spatial inequality. This result is independent of the level and impact of
consumer heterogeneity. This conclusion may be of potential use for policy
makers. Namely, the predictions that globalization is likely to lead to more
spatial inequality (World Bank, 2009) in the future may be reversed if
policies are undertaken to promote inter-regional mobility.

## References

  * (1) Allen, T., Arkolakis, C. (2014), “Trade and the Topography of the Spatial Economy”, _The Quarterly Journal of Economics_ , 129 (3), 1085-1140.
  * Anderson, J., van Wincoop, E. (2004) Anderson, J., van Wincoop, E. (2004), “Trade Costs”, _Journal of Economic Literature_ , 42, 691-751.
  * Bosker, M., Brakman, S., Garretsen, H., & Schramm, M. (2010) Bosker, M., Brakman, S., Garretsen, H., & Schramm, M. (2010), “Adding geography to the new economic geography: bridging the gap between theory and empirics”, _Journal of Economic Geography_ , 10 (6), 793-823.
  * Castro, S.B.S.D., Correia-da-Silva, J., & Gaspar, J.M.. (2022) Castro, S.B.S.D., Correia-da-Silva, J. & Gaspar, J.M.. (2022), “Economic geography meets Hotelling: the home-sweet-home effect”, _Economic Theory_ , 73(1), 183–209.
  * Crozet, M. (2004) Crozet, M. (2004), “Do migrants follow market potentials? An estimation of a new economic geography model”, _Journal of Economic Geography_ , 4(4), 439-458.
  * Forslid, R., Ottaviano, G. (2003) Forslid, R., Ottaviano, G. (2003), “An Analytically Solvable Core-Periphery Model”, Journal of Economic Geography, 3, 229-240.
  * Fujita, M., Krugman, K., Venables, A. (1999) Fujita, M., Krugman, K., Venables, A. (1999), _The Spatial Economy_ , MIT Press.
  * Head, K., Mayer, T. (2004) Head, K., Mayer, T. (2004), “The empirics of agglomeration and trade”, _Handbook of Regional and Urban Economics_ , 4, 2609-2669.
  * Helpman, E. (1998) Helpman, E. (1998), “The size of regions”, in: Pines, D., Zilcha, I., _Topics in Public Economics_ , Cambridge University Press.
  * Ikeda, K., Takayama, Y., Gaspar, J. M., & Osawa, M. (2022) Ikeda, K., Takayama, Y., Gaspar, J. M., & Osawa, M. (2022), “Perturbed cusp catastrophe in a population game: Spatial economics with locational asymmetries”, _Journal of Regional Science_ , 62(4), 961-980.
  * Krugman, P. (1991) Krugman, P. (1991), “Increasing Returns and Economic Geography”, _Journal of Political Economy_ , 99 (3), 483-499.
  * Krugman, P., Elizondo, R. L. (1996) Krugman, P., Elizondo, R. L. (1996), “Trade policy and the third world metropolis”, _Journal of Development Economics_ , 49 (1), 137-150.
  * Murata, Y. (2003) Murata, Y. (2003), “Product diversity, taste heterogeneity, and geographic distribution of economic activities: market vs. non-market interactions”, _Journal of Urban Economics_ , 53, 126-144.
  * Mossay, (2003) Mossay, P. (2003), “Increasing returns and heterogeneity in a spatial economy”, _Regional Science and Urban Economics_ , 33 (4), 419-444.
  * Murata, Y., Thisse, J.-F. (2005) Murata, Y., Thisse, J.-F. (2005), “A simple model of economic geography à la Helpman–Tabuchi”, _Journal of Urban Economics_ , 58 (1), 137-155.
  * Niebuhr, A. (2006) Niebuhr, A. (2006), “Market access and regional disparities”, _The Annals of Regional Science_ , 40(2), 313-334.
  * Pflüger, M. (2004) Pflüger, M. (2004), “A simple, analytically solvable, Chamberlinian agglomeration model”, _Regional Science and Urban Economics_ , 34 (5), 565-573.
  * Tabuchi, T., Thisse, J.-F. (2002) Tabuchi, T., Thisse, J.-F. (2002), “Taste Heterogeneity, Labor Mobility and Economic Geography, _Journal of Development Economics_ , 69, 155-177.
  * Tabuchi, T., Thisse, J. F., & Zhu, X. (2018) Tabuchi, T., Thisse, J. F., & Zhu, X. (2018), “Does technological progress magnify regional disparities?”, _International Economic Review_ , 59(2), 647-663.
  * World Bank (2009) World Bank (2009), “World Development Report 2009 : Reshaping Economic Geography”. World Bank. © World Bank. https://openknowledge.worldbank.org/handle/10986/5991 License: CC BY 3.0 IGO.”

## Appendix A - Wages and freeness of trade

Consider $h=f(w)$ in (8). Let us now define $F(h,w)=f(w)-h=0.$ Differentiation
of $F(h,w)$ yields:

$\dfrac{\partial F(h,w)}{\partial
w}=\dfrac{w^{\sigma}G(w^{\sigma})}{\left[w^{2\sigma}-(w+1)\phi
w^{\sigma}+w\right]^{2}},$ (17)

where:

$G(w^{\sigma})=-\left[\phi(\sigma-1)+(\sigma-1)\phi
w^{2\sigma}-\left(2\sigma-\phi^{2}-1\right)w^{\sigma}\right].$

The derivative in (17) is zero if $G(w^{\sigma})=0$. One can observe that
$G(w^{\sigma})$ has either two (real) zeros given by some
$\left\\{w^{-},w^{+}\right\\},$ or none. Moreover, $G(w^{\sigma})$ is concave
in $w^{\sigma}.$ Since
$G\left(w^{\sigma}=\phi\right)=\sigma\left[\phi(1-\phi^{2})\right]>0$ and
$G\left(w^{\sigma}=\phi^{-1}\right)=-\left[\sigma\phi\left(1-\phi^{-2}\right)\right]>0$,
it must be that $G(w^{\sigma})>0$ for
$w^{\sigma}\in[\phi,\phi^{-1}]\subset(w^{-},w^{+})$.

We now proceed to show that $w\in\left(0,\phi^{1/\sigma}\right)$ and
$w\in\left(\phi^{-1/\sigma},+\infty\right)$ are not defined in $h\in[0,1]$.
Using (9), we have the following:

$\lim_{w\rightarrow 0}h(w)=0;\quad\lim_{w\rightarrow+\infty}h(w)=1;\quad
h(\phi^{1/\sigma})=0;\quad h(\phi^{-1/\sigma})=1;$

Since $dh/dw=\partial F/\partial w$, we have that $h(w)$ is increasing for
$w^{\sigma}\in(w^{-},w^{+})$. The limits above, together with the knowledge
that the zeros of $dh/dw$ lie to the left and right of $[\phi,\phi^{-1}]$,
ensure that $h(w)<0$ for $w\in\left(0,\phi^{1/\sigma}\right)$ and $h(w)>1$ for
$w\in\left(\phi^{-1/\sigma},+\infty\right)$.

Knowing that $\partial F/\partial h=-1$ and $\partial F/\partial w$ are
continuous, by the IFT we can write
$w:h\in[0,1]\subset\mathbb{R}\mapsto\mathbb{R}$ such that $dw/dh$ exists and
$F(h,w(h))=0$.

Using implicit differentiation on (8), we get:

$\dfrac{dw}{dh}=\frac{\left[w^{2\sigma}-(w+1)\phi
w^{\sigma}+w\right]^{2}}{w^{\sigma}G(w^{\sigma})}.$ (18)

This derivative is positive for $w\in[\phi^{1/\sigma},\phi^{-1/\sigma}]$,
which implies that $w(h)$ is increasing in $[0,1]$.

From $F(h,w)=0$ defined above, differentiating with respect to $\phi,$ we get:

$\displaystyle\dfrac{\partial F}{\partial\phi}+\dfrac{\partial F}{\partial
w}\dfrac{dw}{d\phi}$ $\displaystyle=0\iff$ $\displaystyle\dfrac{dw}{d\phi}$
$\displaystyle=-\left.\dfrac{\partial F}{\partial\phi}\right/\dfrac{\partial
F}{\partial w}.$

Using (8) we get:

$\dfrac{\partial
F}{\partial\phi}=\frac{w^{\sigma+1}\left(w^{2\sigma}-1\right)}{\left[w^{2\sigma}-(w+1)\phi
w^{\sigma}+w\right]^{2}},$

which is positive if $h>1/2$, because the latter implies $w>1$ (because $w$ is
increasing in $h$). As a result, we have $dw/d\phi<0$. Therefore, the nominal
wage $w$ is decreasing in the freeness of trade when $L$ is the largest region
($h>1/2$).

## Appendix B - Interior Equilibria

Proof of Lemma 1

###### Proof.

Taking the derivative of (10) with respect to $h$ yields:

$\dfrac{d\Delta
u}{dh}=-\frac{w^{-(\theta+\sigma)}\left\\{w^{\theta}\left(h\phi
w^{1-\sigma}-h+1\right)^{-\frac{\theta+\sigma-2}{\sigma-1}}\left[-h(\sigma-1)\phi
w^{\prime}-w^{\sigma}+\phi
w\right]+\left(hw^{1-\sigma}-h\phi+\phi\right)^{-\left(1+\frac{\theta}{\sigma-1}\right)}\varOmega\right\\}}{(\sigma-1)},$
(19)

where:

$\Omega=-w\left(hw^{1-\sigma}-h\phi+\phi\right)^{\frac{1}{\sigma-1}}\left[-h(\sigma-1)w^{\prime}-\phi
w^{\sigma}+w\right]-(\sigma-1)w^{\sigma}w^{\prime}\left(hw^{1-\sigma}-h\phi+\phi\right)^{\frac{\sigma}{\sigma-1}},$

Substituting for $w^{\prime}$ using (18) from Appendix A, and given that any
partial equilibrium $h^{*}$ must satisfy the short-run equilibrium condition
given by (8), we reach:

$\left.\dfrac{d\Delta
u}{dh}\right|_{h=h^{*}}=\zeta\left\\{\varphi\left[\dfrac{w^{\sigma-1}+w^{-\sigma}-(1+w^{-1})\phi}{1-\phi^{2}}\right]^{\frac{\theta-1}{\sigma-1}}-\psi\left[\dfrac{w^{2\sigma}-(w+1)\phi
w^{\sigma}+w}{w\left(1-\phi^{2}\right)}\right]^{\frac{\theta+\sigma-2}{\sigma-1}}\right\\},$
(20)

with:

$\displaystyle\zeta$
$\displaystyle=\dfrac{w^{-(\theta+\sigma+1)}}{(\sigma-1)\left[(\sigma-1)\phi+(\sigma-1)\phi
w^{2\sigma}+\left(-2\sigma+\phi^{2}+1\right)w^{\sigma}\right]};$
$\displaystyle\varphi$
$\displaystyle=w\left\\{\phi\left[\sigma+(\sigma-1)w\right]w^{\sigma}-2\sigma
w+w\right\\}\left[w^{2\sigma}-(w+1)\phi w^{\sigma}+w\right];$
$\displaystyle\psi$
$\displaystyle=\left(1-\phi^{2}\right)w^{\theta+\sigma+1}\left[(2\sigma-1)w^{\sigma}-\sigma(w+1)\phi+\phi\right].$

Since an interior equilibrium is stable if $d\Delta
u/dh<d\left[t(h)-t(1-h)\right]/dh$ at $h=h^{*}$, it is stable if (20) is
negative, which concludes the proof. ∎

Proof of Proposition 1:

Assume that partial agglomeration $h^{*}=h\in\left(\frac{1}{2},1\right)$
exists and is stable. Then, according to Castro et al. (2020, Prop. 9, pp.
197), it becomes more symmetric if $\frac{d\Delta u}{d\phi}<0$. We have:

$\dfrac{d\Delta u}{d\phi}=\frac{\partial\Delta
u}{\partial\phi}+\frac{\partial\Delta u}{\partial w}\frac{dw}{d\phi},$

which, by using (8) and (10), becomes:

$-\frac{w}{a_{3}}\left\\{a_{1}\left[-\frac{\left(\phi^{2}-1\right)w^{\sigma+1}}{w^{2\sigma}-(w+1)\phi
w^{\sigma}+w}\right]^{-\frac{\theta+\sigma-2}{\sigma-1}}+a_{2}\left[\frac{w-w\phi^{2}}{w^{2\sigma}-(w+1)\phi
w^{\sigma}+w}\right]^{-\frac{\theta+\sigma-2}{\sigma-1}}\right\\}$ (21)

where:

$\displaystyle a_{1}$ $\displaystyle=w^{1-\theta}\left(\phi
w^{\sigma}-1\right)\left[-2\sigma+2(\sigma-1)\phi
w^{\sigma}+\phi^{2}+1\right]$ $\displaystyle a_{2}$
$\displaystyle=w^{-\sigma}\left(w^{\sigma}-\phi\right)\left[2(\sigma-1)\phi+\left(-2\sigma+\phi^{2}+1\right)w^{\sigma}\right]$
$\displaystyle a_{3}$
$\displaystyle=(\sigma-1)\left[w^{1-\sigma}+w^{\sigma}-(w+1)\phi\right]\left[(\sigma-1)\phi+(\sigma-1)\phi
w^{2\sigma}+\left(-2\sigma+\phi^{2}+1\right)w^{\sigma}\right]$

Notice that (21) is equivalent to

$-\frac{w\left[\frac{w-w\phi^{2}}{w^{2\sigma}-(w+1)\phi
w^{\sigma}+w}\right]^{-\frac{\theta+\sigma-2}{\sigma-1}}\left(a_{1}w^{-\sigma\frac{\theta+\sigma-2}{\sigma-1}}+a_{2}\right)}{a_{3}}.$

Thus, the sign of (21) is given by the sign of:

$-\frac{a_{1}w^{-\sigma\frac{\theta+\sigma-2}{\sigma-1}}+a_{2}}{a_{3}}.$

Since $w>1$ for $h\in(\frac{1}{2},1]$, we have $a_{3}<0$ because
$\left(w^{1-\sigma}+w^{\sigma}-(w+1)\phi\right)>0$ and

$\left[(\sigma-1)\phi+(\sigma-1)\phi
w^{2\sigma}+\left(-2\sigma+\phi^{2}+1\right)w^{\sigma}\right]<0.$

Next, evaluating the numerator at $\theta=1$, we get:

$w^{-\sigma}\left(\phi w^{\sigma}-1\right)\left[-2\sigma+2(\sigma-1)\phi
w^{\sigma}+\phi^{2}+1\right]+w^{-\sigma}\left(w^{\sigma}-\phi\right)\left[2(\sigma-1)\phi+\left(-2\sigma+\phi^{2}+1\right)w^{\sigma}\right]<0.$

Differentiating the numerator with respect to $\theta$ yields:

$\frac{\ln(w)w^{\frac{\theta-1}{\sigma-1}+\sigma}\left(\phi
w^{\sigma}-1\right)\left[-2\sigma+2(\sigma-1)\phi
w^{\sigma}+\phi^{2}+1\right]}{\sigma-1},$

which is positive because $1<w^{\sigma}<1/\phi$ and
$\left[-2\sigma+2(\sigma-1)\phi w^{\sigma}+\phi^{2}+1\right]<0$. Thus, the
numerator of (21) remains negative for all $\theta<1$. We conclude that (21)
is negative and thus $\frac{d\Delta u}{d\phi}<0$ for any $\theta\in[0,1]$.
This concludes the proof. $\square$

## Appendix C - Logit heterogeneity

Proof of Proposition 2:

When $t(x)=\mu\ln x$, the RHS of expression (20) is given by:

$t^{\prime}\left(h^{*}\right)-t^{\prime}\left(1-h^{*}\right)=\mu\left[\dfrac{1}{h(1-h)}\right].$

Replacing in (20) and rearranging, partial agglomeration is stable if:

$-\frac{w^{-\sigma-1}\left[w^{2\sigma}-(w+1)\phi
w^{\sigma}+w\right]^{2}\Gamma}{(\sigma-1)\left(w^{\sigma}-\phi\right)\left(1-\phi
w^{\sigma}\right)\left[(\sigma-1)\phi+(\sigma-1)\phi
w^{2\sigma}+\left(-2\sigma+\phi^{2}+1\right)w^{\sigma}\right]}<0,$ (22)

where

$\displaystyle\Gamma$
$\displaystyle=\phi\left[\mu(\sigma-1)^{2}-2\sigma+1\right]+\phi\left[\mu(\sigma-1)^{2}-2\sigma+1\right]w^{2\sigma}-$
$\displaystyle-w^{\sigma}\left\\{\phi^{2}\left[-(\mu+2)\sigma+\mu+1\right]+(2\sigma-1)\left[\mu(\sigma-1)-1\right]\right\\}.$

The numerator of (22) except $\Gamma$ is positive, and so is
$\left(w^{\sigma}-\phi\right)\left(1-\phi w^{\sigma}\right)$. We have:

$(\sigma-1)\phi+(\sigma-1)\phi
w^{2\sigma}+\left(-2\sigma+\phi^{2}+1\right)w^{\sigma}<0.$

Therefore, (20) holds if and only if $\Gamma<0$, which gives:

$\mu>\mu_{p}\equiv\frac{(2\sigma-1)(w^{\sigma}-\phi)(1-w^{\sigma}\phi)}{(\sigma-1)\left[\left(2\sigma-\phi^{2}-1\right)w^{\sigma}-(\sigma-1)\phi
w^{2\sigma}-(\sigma-1)\phi\right]}.$

The condition for partial agglomeration in (22) holds for any interior
equilibrium including symmetric dispersion $h^{*}=1/2$. At symmetric
dispersion we have $w=1$ and the condition above simplifies to:

$\mu>\mu_{d}\equiv\dfrac{(2\sigma-1)(1-\phi)}{(\sigma-1)(2\sigma+\phi-1)}.$

The derivative of $\mu_{p}$ with respect to $X\equiv w^{\sigma}$ equals:

$\dfrac{\partial\mu_{p}}{\partial
X}=\frac{\sigma(2\sigma-1)\left(X^{2}-1\right)\phi\left(\phi^{2}-1\right)}{(\sigma-1)\left[(\sigma-1)\left(X^{2}+1\right)\phi-2\sigma
X+X\phi^{2}+X\right]^{2}}<0.$

We can also see that $\mu_{p}$ approaches zero as $w^{\sigma}$ approaches
$\phi^{-1}$. Therefore, we conclude that $0<\mu_{p}<\mu_{d}$. Under the
assumption that only one partial agglomeration equilibrium $h^{*}\in(1/2,1)$
exists, we have two possibilities: (i) if $\mu\in(\mu_{p},u_{d})$, partial
agglomeration is the only stable equilibrium; and (ii) if $\mu>\mu_{d},$
symmetric dispersion is the only stable equilibrium. Assume now, by way of
contradiction, that $\mu\in(0,\mu_{p})$. Then both dispersion and partial
agglomeration are unstable. Since the state space is one dimensional, this
would require agglomeration to be stable. However, we know that agglomeration
is always unstable, which implies that either partial agglomeration or
dispersion are stable. Hence, $\mu\notin(0,\mu_{p})$. Therefore, symmetric
dispersion is the only stable equilibrium if $\mu>\mu_{d}$; otherwise, partial
agglomeration is the only stable equilibrium. $\square$