# Supersolid Devil’s Staircases of Spin-Orbit-Coupled Bosons in Optical
Lattices

Daisuke Yamamoto1<EMAIL_ADDRESS>Kotaro Bannai2 Nobuo
Furukawa2 Carlos A. R. Sá de Melo3 1Department of Physics, College of
Humanities and Sciences, Nihon University, Sakurajosui, Setagaya, Tokyo
156-8550, Japan  2Department of Physics and Mathematics, Aoyama Gakuin
University, Sagamihara, Kanagawa 252-5258, Japan  3School of Physics, Georgia
Institute of Technology, Atlanta, Georgia 30332, USA

###### Abstract

We study the emergence of supersolid Devil’s staircases of spin-orbit coupled
bosons loaded in optical lattices. We consider two- and three-dimensional
systems of pseudo-spin-$1/2$ bosons interacting via local spin-dependent
interactions. These interactions together with spin-orbit coupling produce
length scales that are commensurate to the lattice spacing. This
commensurability leads to Devil’s staircases of supersolids, with fractal
Hausdorff dimensions, which arise from uniform superfluid phases. We show that
umklapp processes are essential for the existence of commensurate supersolids,
and that without them the Devil’s staircase does not exist. Lastly, we
emphasize the generality of our results, suggest experiments that can unveil
these unusual predictions, and discuss potential applications to the case of
87Rb.

Introduction: For the solid phase of quantum fluids the question - Can solids
be superfluids? - has be asked many years ago leggett-1970  in the context of
solid ${}^{4}{\rm He}$, where it was investigated experimentally, but has so
far yielded a negative answer chan-2012 ; chan-2013 ; reppy-2016 . This
question essentially provides the origin of the name supersolid: the solid
state of a quantum fluid possessing correlated defects that exhibit
superfluidity chester-1970 . The existence of supersolidity is a very
important issue with ramifications in the areas of low temperature condensed
matter physics $(^{4}{\rm He})$ beamish-2007 ; maris-2012 , astrophysics
(neutron star cores) rajagopal-2006 ; samuelsson-2007 , and ultracold atoms
and molecules (large spin atoms, dipolar molecules, and spin-orbit coupled
systems) pfau-2019 ; modugno-2019 ; ferlaino-2019 ; ferlaino-2021 ;
modugno-2021 ; lewenstein-2002 ; mitra-2009 ; ketterle-2017 ; ho-2011 ;
stringari-2012 .

In condensed matter physics, studies of supersolidity focus on the emergence
of superfluid properties from a solid with crystalline order. However, in
ultracold atoms and molecules, investigations of supersolidity focus on the
emergence of crystalline order (solid) from a superfluid. Since quantum gases
are low density systems, standard solid phases, like those found in
${}^{4}{\rm He}$, are out of reach. Thus, in the context of ultracold atoms
and molecules, the reverse question needs to be asked: Can a superfluid be a
solid with crystalline order? This implies that the original definition of
supersolidity, described in the opening paragraph, needs to be revised to
represent a state of matter where both superfluid and solid-like crystalline
order coexist. Using the latter concept, some experimental groups have
recently reported the existence of supersolids in ultracold dipolar bosons
with internal magnetic moments ferlaino-2021 ; modugno-2021 , following
earlier experimental indications of at least metastable supersolidity
pfau-2019 ; modugno-2019 ; ferlaino-2019 , which had been theoretically
suggested as a compromise phase between superfluidity and Wigner-
crystalization in dipolar Bose gases mitra-2009 . These recent experiments
pfau-2019 ; modugno-2019 ; ferlaino-2019 ; ferlaino-2021 ; modugno-2021  have
estimulated a flurry of recent theoretical work on supersolid phases of
dipolar bosons pohl-2019 ; blakie-2020 ; recati-2020 ; stringari-2020 ;
reatto-2021 ; pfau-2021a ; pfau-2021b .

While recent experimental and theoretical work about supersolids in ultracold
quantum gases focused on continuum and trapped (harmonic and boxed) dipolar
bosons pfau-2019 ; modugno-2019 ; ferlaino-2019 ; ferlaino-2021 ; modugno-2021
; pohl-2019 ; blakie-2020 ; recati-2020 ; stringari-2020 ; reatto-2021 ;
pfau-2021a ; pfau-2021b , experimental investigation of supersolidity of
dipolar bosons in optical lattices is still lacking, albeit the existence of
early theoretical work lewenstein-2007 ; sun-2007 ; danshita-2009 ;
pupillo-2010 ; yamamoto-2012  describing the existence of insulating,
superfluid and supersolid phases. In contrast to dipolar systems, we propose
experiments to create and detect supersolid phases of spin-orbit coupled (SOC)
bosons loaded in optical lattices, similar to harmonically trapped bosons
ketterle-2017 . Furthermore, we theoretically study the emergence of
supersolid Devil’s staircases for pseudo-spin-1/2 bosons via local spin-
dependent interactions, like 87Rb.

Experimental Proposal: We propose two experimental setups that could be used
to create a Devil’s staircase of supersolid phases of spin-1/2 bosons with SOC
in optical lattices. The simplest case is the creation of either a two-
dimensional (2D) square or a three-dimensional (3D) cubic optical lattice with
the application of two counter-propagating Raman beams spielman-2011  parallel
to the optical lattice $xy$ plane, but making angle $\theta$ with respect to
the $x$ axis. The second experimental setup involves the utilization of radio-
frequency chips spielman-2013  or monolithic microwave integrated circuits
(MMICs) sademelo-2020 , where the axis of the spin-dependent momentum transfer
can be changed from the $x$ direction, through a relative rotation of the
device with respect to the optical lattice.

Hamiltonian: To investigate the supersolid phases of spin-orbit coupled bosons
in optical lattices, we consider the Hamiltonian for a 2D square lattice

$\displaystyle\hat{\mathcal{H}}$ $\displaystyle=$ $\displaystyle\sum_{\langle
i,j\rangle}\left(\hat{\bm{b}}_{i}^{\dagger}\mathcal{T}_{ij}\hat{\bm{b}}_{j}+{H.c.}\right)+\sum_{i}\hat{\bm{b}}_{i}^{\dagger}\mathcal{M}\hat{\bm{b}}_{i}$
(1)
$\displaystyle+\sum_{is}\frac{U_{ss}}{2}\hat{n}_{is}(\hat{n}_{is}-1)+U_{\uparrow\downarrow}\sum_{i}\hat{n}_{i\uparrow}\hat{n}_{i\downarrow},$

where $\hat{\bm{b}}_{i}=(\hat{b}_{i\uparrow}~{}\hat{b}_{i\downarrow})^{T}$
denotes the annihilation operators of bosons with internal state (pseudo-spin)
$s=\uparrow,\downarrow$ at site $i$ of a square optical lattice with lattice
vectors $\bm{a}_{1}=(a,0)$ and $\bm{a}_{2}=(0,a)$, and
$\hat{n}_{is}=\hat{b}_{is}^{\dagger}\hat{b}_{is}$ counts the local number of
$s$ bosons. The $2\times 2$ matrices are
$\mathcal{T}_{ij}=-t\exp[-i{\bm{\sigma}}_{z}\bm{k}_{T}\cdot(\bm{r}_{i}-\bm{r}_{j})]$
and $\mathcal{M}=-\mu{\bf
1}+\frac{\hbar\Omega}{2}{\bm{\sigma}}_{x}+\frac{\hbar\delta}{2}{\bm{\sigma}}_{z}$
with ${\bm{\sigma}}$ being the Pauli matrices. We consider the equal-Rashba-
Dresselhaus SOC spielman-2011  with momentum transfer
$\bm{p}_{T}=\hbar\bm{k}_{T}=\hbar
k_{T}(\cos\theta\bm{e}_{x}+\sin\theta\bm{e}_{y})$ along the direction tilted
from the lattice $x$-axis by angle $\theta$ in the $xy$ plane. Here,
$\bm{e}_{x}$ ($\bm{e}_{y}$) denotes the unit vector in the $x$ ($y$) direction
of the lattice. The Hamiltonian also includes the Rabi coupling $\Omega$ and
detuning $\delta$, as well as the standard Bose-Hubbard parameters: nearest-
neighbor hopping $t$, chemical potential $\mu$, and intraspin ($s=s^{\prime}$)
and interspin ($s\neq s^{\prime}$) onsite repulsions $U_{ss^{\prime}}>0$ with
$U_{\uparrow\downarrow}^{2}<U_{\uparrow\uparrow}U_{\downarrow\downarrow}$ to
prevent phase separation. We focus on the 2D case, but the 3D case is
analogous, in particular, if one uses the same tilt angle $\theta$ in the $xy$
plane. We explore the model of Eq. (1) in the regime dominated by $t$ in
comparison to $U_{ss^{\prime}}$, and investigate the emergence of supersolids
from superfluid phases.

Figure 1:  (a) Excitation spectra of non-interacting particles
$E_{\bm{k}_{\pm}}$ measured from $\mu$ and (b) the distributions of the BEC
components for $\hbar\Omega=3t$, $k_{T}=0.57\pi/a$, $\theta=0.15\pi$,
$\delta=0$, $U_{\uparrow\uparrow}=U_{\downarrow\downarrow}=U=10t/\rho$, and
$U_{\uparrow\downarrow}=0.9U$. The blue and orange bars in (b) represent the
BEC amplitude $|\psi_{n,s}|$ for $s=\uparrow$ and $s=\downarrow$,
respectively. (c) The stripe pattern of the particle density
$\sum_{s}|\Psi_{is}|^{2}$, resulting from the multiple BECs of (b), in the
real space. The density modulation is $\sim\pm 7.1$ percent of the average
filling factor $\rho\gg 1$ for the lattice sites with the lighter and darker
colors, respectively.

Multiple Condensates: Diagonalizing the Hamiltonian Eq. (1) with
$U_{ss^{\prime}}=0$, we obtain the excitation spectra of non-interacting
particles:

$\displaystyle
E_{\bm{k}\pm}=\frac{\epsilon_{\bm{k}\uparrow}+\epsilon_{\bm{k}\downarrow}}{2}-\mu\pm\sqrt{\left(\frac{\epsilon_{\bm{k}\uparrow}-\epsilon_{\bm{k}\downarrow}}{2}\right)^{2}+\left(\frac{\hbar\Omega}{2}\right)^{2}}$
(2)

with
$\epsilon_{\bm{k}s}=-2t[\cos(k_{x}a+\tau_{s}{k_{T}a\cos\theta})+\cos(k_{y}a+\tau_{s}{k_{T}a\sin\theta})]+\tau_{s}{\hbar\delta/2}$
($\tau_{\uparrow}=1$ and $\tau_{\downarrow}=-1$). We are interested in the
situation where $\hbar|\Omega|/t$ and $\hbar|\delta|/t$ are sufficiently small
to create two minima in the lower branch $E_{\bm{k}-}$ within the first
Brillouin zone (BZ). For instance, this is achieved for
$\hbar|\Omega|/t<4\sum_{j=1,2}|\sin\bm{k}_{T}\cdot\bm{a}_{j}\tan\bm{k}_{T}\cdot\bm{a}_{j}|$,
when $\delta=0$. In this case, when the temperature is sufficiently low, the
particles form Bose-Einstein condensates (BECs) with wave vectors at the two
minima of $E_{\bm{k}-}$, say, $\bm{k}_{1}$ and $\bm{k}_{2}$. However, in the
presence of interactions $U_{ss{{}^{\prime}}}$, the momenta where condensation
occurs are modified to ${\widetilde{\bm{k}}}_{1}$ and
${\widetilde{\bm{k}}}_{2}$. These momenta can be parametrized as
$\widetilde{{\bm{k}}}_{1}=-\bar{\bm{q}}+\delta{\bm{q}}$ and
$\widetilde{{\bm{k}}}_{2}=\bar{\bm{q}}+\delta{\bm{q}}$. The deviation
$\delta\bm{q}$ reflects the parity asymmetry caused by non-zero detuning
$\delta$ and/or by broken $Z_{2}$-symmetry interactions
$U_{\uparrow\uparrow}\neq U_{\downarrow\downarrow}$.

When BECs are formed at $\widetilde{\bm{k}}_{1}$ and $\widetilde{\bm{k}}_{2}$,
the interference of the two fundamental matter waves produces higher harmonics
with wave vectors differing by an integer multiple of $2\bar{\bm{q}}$. This
indicates that the expectation value of $\langle\hat{b}_{is}\rangle$ acquires
a spatial modulation of the form

$\langle\hat{b}_{is}\rangle=\Psi_{is}=\sqrt{\rho}\sum_{n}\psi_{ns}e^{i\bm{q}_{n}\cdot\bm{r}_{i}}$
(3)

with $\bm{q}_{n}=(2n-1)\bar{\bm{q}}+\delta{\bm{q}}$ ($n\in\mathbb{Z}$). Here,
$\rho=\sum_{is}|\Psi_{is}|^{2}/M$ is the particle filling per site, where $M$
is the number of lattice sites. The Fourier amplitudes $\psi_{ns}$ represent
BEC order parameters at each momentum $\bm{q}_{n}$ and are normalized to one,
that is, $\sum_{ns}|\psi_{ns}|^{2}=1$. The subscript $n$ labels the harmonic
components (HCs), for example, $n=0,1$ correspond to the two fundamental
matter waves with wave vectors $\mp\bar{\bm{q}}+\delta{\bm{q}}$, and $n=-1,2$
are the second harmonics with $\mp 3\bar{\bm{q}}+\delta{\bm{q}}$, and so on.

Using Eq. (3), we minimize the variational energy per particle

$\displaystyle\frac{E_{0}}{M\rho}$ $\displaystyle=$
$\displaystyle\sum_{n}\left(\begin{array}[]{cc}\psi_{n\uparrow}^{\ast}&\psi_{n\downarrow}^{\ast}\end{array}\right)\left(\begin{array}[]{cc}\epsilon_{\bm{q}_{n}\uparrow}-\mu&\hbar\Omega/2\\\
\hbar\Omega/2&\epsilon_{\bm{q}_{n}\downarrow}-\mu\end{array}\right)\left(\begin{array}[]{c}\psi_{n\uparrow}\\\
\psi_{n\downarrow}\end{array}\right)$ (10)
$\displaystyle+\sum^{\prime}_{n_{1}+n_{2}=n_{3}+n_{4}}\sum_{ss^{\prime}}\frac{U_{ss^{\prime}}\rho}{2}\psi_{n_{1}s}^{\ast}\psi_{n_{2}s^{\prime}}^{\ast}\psi_{n_{3}s^{\prime}}\psi_{n_{4}s},$

with respect to the order parameters $\psi_{ns}$ and the wave vectors
$\bar{\bm{q}}$ and $\delta{\bm{q}}$ under the condition
$\sum_{ns}|\psi_{ns}|^{2}=1$. The sum in the interaction term is over all
possible subsets of the HCs that satisfy momentum conservation implying the
restriction $n_{1}+n_{2}=n_{3}+n_{4}$. When $2\bar{\bm{q}}$ is commensurate to
the lattice spacing, we also need to consider umklapp scattering processes
with momentum conservation modulo the reciprocal lattice vectors
$\bm{G}_{1}=(2\pi/a,0)$ and $\bm{G}_{2}=(0,2\pi/a)$, as we shall explain
later.

In Fig. 1(b), we show the ground-state distribution of the order parameter
amplitude $|\psi_{ns}|$, in the first BZ, for parameters given in the caption.
Remarkably, we find that the $({\bar{q}}_{x},{\bar{q}}_{y})$ components of
$\bar{\bm{q}}$, in units of $\pi/a$, may be rational numbers even when the SOC
momentum components $(k_{Tx},k_{Ty})$, in units of $\pi/a$, are irrational
numbers, while the components of $\delta\bm{q}$ can have any real value. For
the parameters of Fig. 1(b), where $\delta=0$ and
$U_{\uparrow\uparrow}=U_{\downarrow\downarrow}$, then $\delta\bm{q}=\bm{0}$,
and the fundamental BEC wave vectors are $\pm\bar{\bm{q}}=\pm(1/2,1/4)\pi/a$.
The fundamental wave vectors for the density modulation
$\sum_{s}|\Psi_{is}|^{2}$ are then $\bm{Q}=\pm 2\bar{\bm{q}}$, which are
expressed as $\bm{Q}=\pm\left[(1/2)\bm{G}_{1}+(1/4)\bm{G}_{2}\right]$, showing
a modulation period $\Lambda_{x}=2a$ $(\Lambda_{y}=4a)$ along the $x$ $(y)$
direction, as seen in Fig. 1(c).

When $\bar{q}_{x}$ and $\bar{q}_{y}$ are commensurate with $\pi/a$, the number
of the HCs is finite. For example, in Fig. 1(b), the wave vectors of second-
harmonic components are given as $\bm{q}_{-1,2}=\pm
3\bar{\bm{q}}=\pm(3/2,3/4)\pi/a$, which are equivalent to $\pm(-1/2,3/4)\pi/a$
in the first BZ. In this case, third and higher harmonics are reduced to
either the fundamental or the second harmonic wave vectors due to momentum-
space periodicity, and thus it is sufficient to consider up to the second
components. For a general commensurate wave vector
$\bar{\bm{q}}=(\xi_{1}/\eta_{1},\xi_{2}/\eta_{2})\pi/a$, with relatively prime
integers $\xi_{\ell}$ and $\eta_{\ell}$, the number of independent HCs is
given by $N_{\rm HC}={\rm LCM[\eta_{1},\eta_{2}]}$, where LCM means least
common multiple. The interference of BECs with $N_{\rm HCs}$ components
produces striped interference patterns (supersolids), with wavelengths
$\Lambda_{x}=\eta_{1}a$ and $\Lambda_{y}=\eta_{2}a$, in the density profile
$|\Psi_{is}|^{2}$.

Commensurability and umklapp scattering: Next, we show that the difference
vector $2\bar{\bm{q}}$ between the two fundamental matter waves
$\bm{q}_{0,1}=\pm\bar{\bm{q}}+\delta\bm{q}$ must always be of the form
$\frac{\xi_{1}}{\eta_{1}}\bm{G}_{1}+\frac{\xi_{2}}{\eta_{2}}\bm{G}_{2}$ in the
ground state, as seen from the conditions satisfied by the phases $\phi_{ns}$
of the order parameters $\psi_{ns}=|\psi_{ns}|e^{i\phi_{ns}}$ for the
minimization of the ground-state energy in Eq. (10). First, from the spin-flip
terms in Eq. (10), the relative phase between $\psi_{n\uparrow}$ and
$\psi_{n\downarrow}$ with the same $n$ is determined by the sign of $\Omega$
to be $\phi_{n\uparrow}-\phi_{n\downarrow}=\pi$
$(\phi_{n\uparrow}-\phi_{n\downarrow}=0)$ when $\Omega>0$ ($\Omega<0$).
Second, given the previous condition, the interaction terms involving both
intraspin $U_{ss}$ and interspin $U_{ss^{\prime}}$ (with $s\neq s^{\prime}$)
interactions produce the factor

$\displaystyle{\cal
A}_{ss^{\prime}}\cos[\phi_{n_{1}s}+\phi_{n_{2}s}-\phi_{n_{3}s}-\phi_{n_{4}s}],$
(11)

with positive interaction coefficients $({\cal A}_{ss^{\prime}}>0)$. Here, we
choose the global phase to be
$\phi_{0\uparrow}=-\phi_{1\uparrow}=-\bar{\phi}/2$, without loss of
generality. The second harmonic components $n=2$ and $n=-1$ arise from the
scattering processes of the type $(n_{1},n_{2};n_{3},n_{4})=(2,0;1,1)$ and
$(-1,1;0,0)$, respectively. Thus, their phases must be
$\phi_{2\uparrow}=3\bar{\phi}/2+\pi$ and
$\phi_{-1\uparrow}=-3\bar{\phi}/{2}+\pi$ to minimize the interaction energy in
Eq. (11). Analogously, the momentum conservation $n_{1}+n_{2}=n_{3}+n_{4}$ and
the minimization of the interaction energy lead to the conclusion that the up-
spin phases of the HCs satisfy

$\phi_{n\uparrow}=\frac{2n-1}{2}\bar{\phi}+\arccos[(-1)^{\frac{|2n-1|-1}{2}}]$
(12)

with $n\in\mathbb{Z}$, while the down-spin phases are
$\phi_{n\downarrow}=\phi_{n\uparrow}-\pi$ for $\Omega>0$ or
$\phi_{n\downarrow}=\phi_{n\uparrow}$ for $\Omega<0$. The last degree of
freedom is the relative phase $\bar{\phi}$ between the two fundamental BECs,
discussed next.

When the difference vector between the two fundamental matter waves is
commensurate to the lattice spacing, that is,
$2\bar{\bm{q}}=\frac{\xi_{1}}{\eta_{1}}\bm{G}_{1}+\frac{\xi_{2}}{\eta_{2}}\bm{G}_{2}$,
umklapp scattering processes with the total momentum transfer equal to
reciprocal lattice vectors $n_{1}+n_{2}-n_{3}-n_{4}=\pm{\rm
LCM[\eta_{1},\eta_{2}]}=\pm N_{\rm HC},$ must be considered in the second sum
of Eq. (10). Using Eqs. (11) and (12), it is easy to see that the sum of the
umklapp scattering terms is proportional to ${\cal B}\cos[N_{\rm
HC}\bar{\phi}]$, which shows that the umklapp scattering can always reduce the
energy by taking an appropriate value of $\bar{\phi}$ such that $\cos[N_{\rm
HC}\bar{\phi}]=-1$ $(+1)$ when ${\cal B}>0$ $({\cal B}<0)$. This implies that
the supersolid ground state exhibits $N_{\rm HC}$-fold degeneracy with respect
to the choice of $\bar{\phi}$, and thus breaks $Z_{N}$ symmetry (with
$N=N_{\rm HC}$), in addition to the breaking ${\rm U(1)}$ symmetry associated
with the global phase fixed earlier.

Figure 2:  Contour plots of the values of $\bar{q}_{x}$ and $\bar{q}_{y}$ in
the $k_{T}\cos\theta$ vs. $k_{T}\sin\theta$ plane for (a)-(b)
$\hbar\Omega=1.2t$ and (c)-(d) $\hbar\Omega=3t$ ($\delta=0$,
$U_{\uparrow\uparrow}=U_{\downarrow\downarrow}=U=8t/\rho$, and
$U_{\uparrow\downarrow}=0.9U$ for both). The gray areas represent the regions
where a single-BEC superfluid is the ground state.

Supersolid Devil’s staircase: As discussed above, commensurate ground states
with rational values of $\bar{q}_{x}$ and $\bar{q}_{y}$ (in units of $\pi/a$)
are favored against incommensurate ones when umklapp processes are important.
This indicates that $\bar{q}_{x}$ and $\bar{q}_{y}$ develop plateaus at
rational values for numerous intervals of input variables $k_{T_{x}}$ and
$k_{T_{y}}$, thus forming a Devil’s staircase structure with an infinite
number of steps. In Fig 2, we show the values of $\bar{q}_{x}$ and
$\bar{q}_{y}$ as functions of $k_{T_{x}}=k_{T}\cos\theta$ and
$k_{T_{y}}=k_{T}\sin\theta$ for two sets of parameters given in the caption.
We show only the first quadrant of the first BZ since the function
$\bar{q}_{x}$ ($\bar{q}_{y}$) is odd in $k_{T_{x}}$ ($k_{T_{y}}$) and even in
$k_{T_{y}}$ ($k_{T_{x}}$).

In all panels of Fig. 3, the parameters used are
$U_{\uparrow\uparrow}=U_{\downarrow\downarrow}=U$,
$U_{\uparrow\downarrow}=0.9U$, and $\delta=0$. In Fig. 3(a), we show
supersolid Devil’s staircases for $\bar{q}_{x}$ versus $k_{T}$ at SOC angle
$\theta=0$. A box-counting analysis for the plateau width of the Devil’s
staircases is shown in Fig. 3(b). The function $L(\epsilon)$ is the difference
between the total width of the staircase and the sum of the plateaus widths
larger than $\epsilon>0$. The slope of the log-log plot of
$L(\epsilon)/\epsilon$ versus $1/\epsilon$ in the limit of $\epsilon\to 0$
gives the Hausdorff fractal dimension $D$ of the system salinas-2014 . If
$D<1$, the incommensurate phases form a fractal set of measure zero, meaning
that we have a complete Devil’s staircase of commensurate (supersolid) phases.
In Fig. 3(c), we show $D$ versus $U\rho/t$ to indicate the fractality of the
staircase. Notice that when interactions tend to zero ($U\rho/t\to 0$), then
the Hausdorff dimension $D\to 1$, meaning that the lines in Fig. 3(a) do not
contain a dense set of plateaus. This reinforces that interactions are
essential for the emergence of the supersolid Devil’s staircase and its
fractality.

Figure 3:  (a) Supersolid Devil’s staircases for $\bar{q}_{x}$ versus $k_{T}$,
in units of $\pi/a$, at SOC angle $\theta=0$. The orange, green, and red lines
are vertically shifted by $0.1$, $0.2$, and $0.3$ with respect to the blue
line to avoid overlap. The black triangles at the endpoints indicate the
emergence of single-BEC superfluid phases. (b) Plots of the function
$L(\epsilon)/\epsilon$ vs. $1/\epsilon$ characterizing the plateau widths
$\epsilon$ of the Devil’s staircases. (c) Hausdorff fractal dimension $D$ as
functions of $U\rho/t$.

Experimental Detection: For fixed $t$, $U_{ss^{\prime}}$, $\Omega$, and
$\delta$, there are two experimental parameters that can be adjusted: the
magnitude of the momentum transfer $k_{T}$ and the angle $\theta$ between the
direction of $\bm{k}_{T}$ and the $x$ axis of the optical lattice. However, it
is easier to vary the tilt angle $\theta$, as changing $k_{T}$ requires a
different laser wavelength for the Raman setup, or a different radio-frequency
(microwave) wavelength in the atom-chip (MMIC) configuration. Thus, in Fig. 4,
we show examples of ${\bar{q}}_{x}$ and ${\bar{q}}_{y}$ for fixed $k_{T}$ and
changing $\theta$. In Fig. 4(a), $k_{T}=2.02\pi/a$, but the vector
$\bm{k}_{T}$ is closer to pointing along the $k_{x}$ direction since
$73.5^{\rm o}<\theta<84.5^{\rm o}$ [see Fig. 4(c)]. In this case, the value of
${\bar{q}}_{y}$ is nearly (or exactly) zero, while ${\bar{q}}_{x}$ take
fractional values with the largest steps being at
${\bar{q}}_{x}=\\{1/2,1/3,1/4,1/5\\}\pi/a$. These are supersolid stripes,
where the density $\sum_{s}|\Psi_{is}|^{2}$ is uniform along the $y$ direction
and modulated along the $x$ direction with period $\\{2,3,4,5\\}a$. In Fig.
4(b), $k_{T}=0.72\pi/a$, but the vector $\bm{k}_{T}$ is closer to the diagonal
in the first quadrant of the first BZ [see Fig. 4(d)]. In this case, the most
prominent commensurate (supersolid) phases have the largest ladder steps
characterized by the ordered pairs
$({\bar{q}}_{x},{\bar{q}}_{y})=\\{(1/3,2/3);(2/5,3/5);(1/2,1/2);(2/3,1/3);(3/5,2/5)\\}$
in units of $\pi/a$.

Figure 4:  Supersolid Devil’s staircases for $\bar{q}_{x},\bar{q}_{y}$ vs. SOC
angle $\theta$ for fixed values of $k_{T}$. Panel (a) [(b)] has same
parameters as in Figs. 2(a)-(b) [Figs. 2(c)-(d)] along the arc with
$k_{T}=2.02\pi/a$ $[k_{T}=0.72\pi/a]$, indicated by the arrow in (c) [(d)].
The triangles have the same meaning as in Fig. 3(a).

The simplest experiments to detect supersolid phases and their staircase
structure are momentum space measurements. Both the momentum distribution
$n(\bm{k})$ in time of flight bloch-2002  and the structure factor $S(\bm{q})$
obtained from Bragg spectroscopy hulet-2010  can reveal the fundamental and
higher-order momentum components of the order parameter density in the
supersolid phases. In addition, the real-space periodic modulations of the
supersolids, which are commensurate with the underlying lattice structure, can
be detected, in principle, using quantum gas microscopes greiner-2009 ;
bloch-2016 ; bakr-2018  or magnifiers weitenberg-2021 .

When atoms have anisotropic interactions $U_{\uparrow\uparrow}\neq
U_{\downarrow\downarrow}$, then $Z_{2}$ symmetry is broken. This is the case
for two hyperfine states of ${}^{87}\rm{Rb}$ atoms,
$|\\!\uparrow\rangle=|F=1,m_{F}=0\rangle$ and
$|\\!\downarrow\rangle=|F=1,m_{F}=-1\rangle$, where
$U_{\downarrow\downarrow}\approx U_{\uparrow\downarrow}\approx
0.995U_{\uparrow\uparrow}$ spielman-2011 ; widera-2006 . The broken spin
symmetry can be compensated by the detuning
$\delta=\delta_{0}=-\rho(U_{\uparrow\uparrow}-U_{\downarrow\downarrow})/2$,
since the anisotropy is small
$(|U_{\uparrow\uparrow}-U_{\downarrow\downarrow}|\ll\sum_{ss^{\prime}}U_{ss^{\prime}})$
stringari-2012 . In this case, the shift $\delta\bm{q}$ of the BEC momenta is
negligible and the results for $U_{\uparrow\uparrow}\neq
U_{\downarrow\downarrow}$ with $\delta=\delta_{0}$ are essentially identical
to the results for $U_{\uparrow\uparrow}\to{\overline{U}}$,
$U_{\downarrow\downarrow}\to{\overline{U}}$, and $\delta_{0}\to 0$, where
${\overline{U}}=(U_{\uparrow\uparrow}+U_{\downarrow\downarrow})/2$.

Conclusions: We have investigated the existence of supersolid Devil’s
staircases for spin-orbit coupled bosons in optical lattices, proposed
experimental setups to investigate the phenomenon and suggested detection
techniques for their direct observation. Furthermore, we showed that the
cascade of supersolid phases occurs due to the commensurability of the spatial
modulation of the order parameter with respect to the underlying lattice,
which is induced by spin-orbit coupling and stabilized by interactions
(umklapp processes). This work opens the door for experimental investigations
of supersolidity in ${}^{87}{\rm Rb}$ and other bosonic systems, and, more
generally, provides a fundamental idea for understanding the commensuration of
ordering vectors with possible applications to FFLO superconductivity
Kinnunen-18  and chiral magnets tsuruta-16 .

###### Acknowledgements.

This work was supported by JSPS KAKENHI Grant Nos. 18K03525 and 21H05185
(D.Y.), JST CREST Grant No. JPMJCR1673 (D.Y.), and JST PRESTO Grant No.
JPMJPR2118, Japan (D.Y.).

## References

  * (1) A. J. Leggett, Can a solid be superfluid? Phys. Rev. Lett. 25, 1543 (1970).
  * (2) D. Y. Kim and M. H. W. Chan, Absence of supersolidity in solid Helium in porous Vycor glass, Phys. Rev. Lett. 109, 155301 (2012).
  * (3) M. H. W. Chan, R. Hallock, and L. Reatto, Overview on solid ${}^{4}{\rm He}$ and the issue of supersolidity, J. Low Temp. Phys. 172, 317 (2013).
  * (4) A. Eyal, X. Mi, A. V. Talanov, J. D. Reppy, Multiple mode torsional oscillator studies and evidence for supersolidity in bulk ${}^{4}{\rm He}$, PNAS 113, E3203 (2016).
  * (5) G. V. Chester, Speculations on Bose-Einstein condensation and quantum crystals, Phys. Rev. A 2, 256 (1970).
  * (6) J. Day and J. Beamish, Low-temperature shear modulus changes in solid ${}^{4}{\rm He}$ and connection to supersolidity, Nature 450, 853 (2007).
  * (7) H. J. Maris, Effect of elasticity on torsional oscillator experiments probing the possible supersolidity of Helium, Phys. Rev. B 86, 020502(R) (2012).
  * (8) K. Rajagopal and R. Sharma, Crystallography of three-flavor quark matter, Phys. Rev. D 74, 094019, (2006).
  * (9) B. Haskell, N. Andersson, D. I. Jones, and L. Samuelsson, Are neutron stars with crystalline color-superconducting cores relevant for the LIGO experiment? Phys. Rev. Lett. 99, 231101 (2007).
  * (10) F. Böttcher, J. N. Schmidt, M.Wenzel, J. Hertkorn, M. Guo, T. Langen, and T. Pfau, Transient supersolid properties in an array of dipolar quantum droplets, Phys. Rev. X 9, 011051 (2019).
  * (11) L. Tanzi, E. Lucioni, F. Famà, J. Catani, A. Fioretti, C. Gabbanini, R. N. Bisset, L. Santos, and G. Modugno, Observation of a dipolar quantum gas with metastable supersolid properties, Phys. Rev. Lett. 122, 130405 (2019).
  * (12) L. Chomaz, D. Petter, P. Ilzhöfer, G. Natale, A. Trautmann, C. Politi, G. Durastante, R. M. W. van Bijnen, A. Patscheider, M. Sohmen, M. J. Mark, and F. Ferlaino, Long-lived and transient supersolid behaviors in dipolar quantum gases, Phys. Rev. X 9, 021012 (2019).
  * (13) M. A. Norcia, C. Politi, L. Klaus, E. Poli, M. Sohmen, M. J. Mark, R. Bisset, L. Santos, and F. Ferlaino, Two-dimensional supersolidity in a dipolar quantum gas, Nature 596, 357-361 (2021).
  * (14) L. Tanzi, J. G. Maloberti, G. Biagioni, A. Fioretti, C. Gabbanini, and G. Modugno, Evidence of superfluidity in a dipolar supersolid from non-classical rotational inertia, Science 371, 1162 (2021).
  * (15) K. Góral, L. Santos, and M. Lewenstein, Quantum phases of dipolar bosons in optical lattices, Phys. Rev. Lett. 88, 170406 (2002).
  * (16) K. Mitra, C. J. Williams and C. A. R. Sá de Melo, Hexatic, Wigner crystal, and superfluid phases of dipolar bosons, arXiv:0903.4655v1 (2009).
  * (17) J. R. Li, J. Lee, W. Huang, S. Burchesky, B. Shteynas, F. C. Top, A. O. Jamison, and W. Ketterle, A stripe phase with supersolid properties in spin-orbit-coupled Bose-Einstein condensates, Nature 543, 91 (2017).
  * (18) T. L. Ho and S. Zhang, Bose-Einstein condensates with spin-orbit interaction, Phys. Rev. Lett. 107, 150403 (2011).
  * (19) Y. Li, L. P. Pitaevskii, and S. Stringari, Quantum tricriticality and phase transitions in spin-orbit coupled Bose-Einstein condensates, Phys. Rev. Lett. 108, 225301 (2012).
  * (20) Y. C. Zhang, F. Maucher, and T. Pohl, Supersolidity around a critical point in dipolar Bose–Einstein condensates, Phys. Rev. Lett. 123, 015301 (2019).
  * (21) P. B. Blakie, D. Baillie, L. Chomaz, and F. Ferlaino, Supersolidity in an elongated dipolar condensate, Phys. Rev. Research 2, 043318 (2020).
  * (22) A. Gallemí, S. M. Roccuzzo, S. Stringari, and A. Recati, Quantized vortices in dipolar supersolid Bose–Einstein-condensed gases, Phys. Rev. A 102, 023322 (2020).
  * (23) S. M. Roccuzzo, A. Gallemí, A. Recati, and S. Stringari, Rotating a supersolid dipolar gas, Phys. Rev. Lett. 124, 045702 (2020).
  * (24) F. Ancilotto, M. Barranco, M. Pi, and L. Reatto, Vortex properties in the extended supersolid phase of dipolar Bose–Einstein condensates, Phys. Rev. A 103, 033314 (2021).
  * (25) J. Hertkorn, J. N. Schmidt, M. Guo, F. Böttcher, K. S. H. Ng, S. D. Graham, P. Uerlings, H. P. Büchler, T. Langen, M. Zwierlein, and T. Pfau, Supersolidity in two-dimensional trapped dipolar droplet arrays, Phys. Rev. Lett. 127, 155301 (2021).
  * (26) J. Hertkorn, J. N. Schmidt, M. Guo, F. Böttcher, K. S. H. Ng, S. D. Graham, P. Uerlings, T. Langen, M. Zwierlein, and T. Pfau, Pattern formation in quantum ferrofluids: from supersolids to superglasses, Phys. Rev. Research 3, 033125 (2021).
  * (27) C. Menotti, C. Trefzger, and M. Lewenstein, Metastable states of a gas of dipolar bosons in a 2D optical lattice, Phys. Rev. Lett. 98, 235301 (2007).
  * (28) S. Yi, T. Li, and C. P. Sun, Novel quantum phases of dipolar Bose gases in optical lattices, Phys. Rev. Lett. 98, 260405 (2007).
  * (29) I. Danshita and C A. R. Sá de Melo, Stability of superfluid and supersolid phases of dipolar bosons in optical lattices, Phys. Rev. Lett. 103, 225301 (2009).
  * (30) B. Capogrosso-Sansone, C. Trefzger, M. Lewenstein, P. Zoller, and G. Pupillo, Quantum phases of cold polar molecules in 2D optical lattices, Phys. Rev. Lett. 104, 125301 (2010).
  * (31) D. Yamamoto, I. Danshita, and C. A. R. Sá de Melo, Dipolar bosons in triangular optical lattices: Quantum phase transitions and anomalous hysteresis, Phys. Rev. A 85, 021601(R) (2012).
  * (32) Y. J. Lin, K. Jiménez-García, and I. B. Spielman, Spin–orbit-coupled Bose–Einstein condensates, Nature 471, 83 (2011).
  * (33) B. M. Anderson, I. B. Spielman, and G. Juzeliunas, Magnetically generated spin-orbit coupling for ultracold atoms, Phys. Rev. Lett. 111, 125301 (2013).
  * (34) D. Lao, C. Raman and C. A. R. Sá de Melo, Nematic-orbit coupling and nematic density waves in Spin-1 condensates, Phys. Rev. Lett 124, 173203 (2020).
  * (35) E. S. Nascimento, J. P. de Lima, and S. R. Salina, Modulated phases and Devil’s staircases in a layered mean-field version of the ANNNI model, Physica A 409, 78 (2014).
  * (36) M. Greiner, O. Mandel, T. Esslinger, T.W. Hänsch, and I. Bloch, Quantum phase transition from a superfluid to a Mott insulator in a gas of ultracold atoms, Nature 415, 39 (2002).
  * (37) T. A. Corcovilos, S. K. Baur, J. M. Hitchcock, E. J. Mueller, and R. G. Hulet, Detecting antiferromagnetism of atoms in an optical lattice via optical Bragg scattering, Physical Review A 81, 013415 (2010).
  * (38) W. S. Bakr, J. Gillen, A. Peng, S. Fölling, and M. Greiner, A quantum gas microscope for detecting single atoms in a Hubbard regime optical lattice, Nature 462, 74 (2009).
  * (39) M. Boll, T. A. Hilker, G. Salomon, A. Omran, J. Nespolo, L. Pollet, I. Bloch, and C. Gross, Spin- and density-resolved microscopy of antiferromagnetic correlations in Fermi-Hubbard chains, Science 353, 1257 (2016).
  * (40) D. Mitra, P. Brown, E. Guardado-Sanchez, S. Kondov, T. Devakul, D. Huse, P. Schauss, and W. Bakr, Quantum gas microscopy of an attractive Fermi-Hubbard system, Nat. Phys. 14, 173 (2018).
  * (41) L. Asteria, H. P. Zahn, M. N. Kosch, K. Sengstock, and C. Weitenberg, Quantum gas magnifier for sub-lattice-resolved imaging of three-dimensional quantum systems, Nature 599, 571 (2021).
  * (42) A. Widera, F. Gerbier, S. Fölling, T. Gericke, O. Mandel, and I. Bloch, Precision measurement of spin-dependent interaction strengths for spin-1 and spin-2 87Rb atoms, New J. Phys. 8, 152 (2006).
  * (43) J. J. Kinnunen, J. E. Baarsma, J.-P. Martikainen, and P. Törmä, The Fulde-Ferrell-Larkin-Ovchinnikov state for ultracold fermions in lattice and harmonic potentials: a review, Rep. Prog. Phys. 81, 046401 (2018).
  * (44) K. Tsuruta, M. Mito, Y. Kousaka, J. Akimitsu, J. Kishine, Y. Togawa, H. Ohsumi, and K. Inoue, Discrete change in magnetization by chiral soliton lattice formation in the chiral magnet Cr1/3NbS2, J. Phys. Soc. Jpn. 85, 013707 (2016).


# Approximating monomials using Chebyshev polynomials

Arvind K. Saibaba111Department of Mathematics, North Carolina State
University, Raleigh, NC. Email<EMAIL_ADDRESS>
###### Abstract

This paper considers the approximation of a monomial $x^{n}$ over the interval
$[-1,1]$ by a lower-degree polynomial. This polynomial approximation can be
easily computed analytically and is obtained by truncating the analytical
Chebyshev series expansion of $x^{n}$. The error in the polynomial
approximation in the supremum norm has an exact expression with an interesting
probabilistic interpretation. We use this interpretation along with
concentration inequalities to develop a useful upper bound for the error.

Keywords— Chebyshev polynomials, Polynomial Approximation, Binomial
Coefficients, Concentration inequalities

##  1 Motivation and Introduction

We are interested in approximating the monomial $x^{n}$ by a polynomial of
degree $0\leq k<n$ over the interval $[-1,1]$. The monomials $1,x,x^{2},\dots$
form a basis for $C[-1,1]$, so it seems unlikely that we can represent a
monomial in terms of lower degree polynomials. In Figure 1, we plot a few
functions from the monomial basis over $[0,1]$; the basis function look
increasingly alike as we take higher and higher powers, i.e., they appear to
“lose independence.” Numerical analysts often avoid the monomial basis in
polynomial interpolation since they result in ill-conditioned Vandermonde
matrices, leading to poor numerical performance in finite precision
arithmetic. This loss of independence means that it is reasonable to
approximate the monomial $x^{n}$ as a linear combination of lower order
monomials, i.e., a lower order polynomial approximation. The natural question
to ask, therefore, is: how small can $k$ be so that a well-chosen polynomial
of degree $k$ can accurately approximate $x^{n}$?

Figure 1: Visualization of a few monomials in the interval $[0,1]$.

The surprising answer to this question is that we can approximate the monomial
$x^{n}$ over $[-1,1]$ by a polynomial of small degree, which we will make
precise. Let $\|f\|_{\infty}=\max_{x\in[-1,1]}|f(x)|$ denote the supremum norm
on $C[-1,1]$ and let $\pi_{k}^{*}(\cdot)$ be the best polynomial approximation
to $x^{n}$ in this norm; that is

$E_{n,k}:=\min_{\pi\in\mathcal{P}_{k}}\|x^{n}-\pi(x)\|_{\infty}=\|x^{n}-\pi_{k}^{*}(x)\|_{\infty},$

where $\mathcal{P}_{k}$ is a vector space of polynomials with real
coefficients of degree at most $k$. The minimizer $\pi_{k}^{*}(\cdot)$ exists
and is unique [9, Chapter 10], but does not have a closed form expression.
Newman and Rivlin [7, Theorem 2] showed that222We briefly mention that the
notation our manuscript differs from [7] in that we reverse the roles of $n$
and $k$.

$\frac{p_{n,k}}{4e}\leq\|x^{n}-\pi_{k}^{*}(x)\|_{\infty}\leq p_{n,k},$ (1)

where the term $p_{n,k}$ is given by the formula

$p_{n,k}=\frac{1}{2^{n-1}}\sum_{j=\lfloor(n+k)/2\rfloor+1}^{n}\binom{n}{j}.$

Since $p_{n,k}$ involves the sum of binomial coefficients, it has a
probabilistic interpretation which we explore in Section 4.

To see why a small $k$ is sufficient, consider the upper bound $p_{n,k}$. In
Section 4 we use the probabilistic interpretation to obtain the following
bound $p_{n,k}\leq 2\exp\left(-{k^{2}}/{2n}\right)$. Suppose we are given a
user-defined tolerance $\epsilon>0$. To ensure

$\|x^{n}-\pi_{k}^{*}(x)\|_{\infty}\leq\epsilon,$

we need to choose $k\geq\sqrt{2n\log(2/\epsilon)}$. The accuracy of the
polynomial approximation is visualized in Figure 2, where in the left panel we
plot the monomial $x^{n}$ for $n=75$ and the best polynomial approximation
$\pi_{k}^{*}$ for $k=5,15,25$. The polynomial $\pi_{k}^{*}$ is computed using
the Remez algorithm, implemented in chebfun [3]. We see that for $k=25$, the
polynomial approximation looks very accurate. In the right panel, we display
$p_{n,k}$, which is the upper bound of the best polynomial approximation, as
well as the upper bound for $p_{n,k}$. We see that $p_{n,k}$ and its upper
bound both have sharp decay with increasing $k$. Numerical evidence in [6]
further confirms this analysis; the authors show that the error $E_{n,k}$
behaves approximately like $\frac{1}{2}\text{erfc}(k/\sqrt{n})$, where erfc is
the complementary error function. .

Figure 2: (left) Approximation of the monomial $x^{n}$ for $n=75$ by
$\pi_{k}^{*}$, (right) $p_{n,k}$ and its upper bound $2\exp(-k^{2}/2n)$. The
visualization is restricted to the interval $[0,1]$.

Polynomial and rational approximations to the monomial has received
considerable attention in the approximation theory community, and surveys of
various results can be found in [8, 6]. Polynomial approximations to high
order monomials have many applications in numerical analysis. This key insight
was exploited by Cornelius Lanczos [5] in his “$\tau$-method” for the
numerical solution of differential equations. For a simulating discussion on
this topic, please see [6]. In numerical linear algebra, this has been
exploited to efficiently compute matrix powers and the Schatten p-norm of a
matrix [1, 4].

In this short note, we show to construct a polynomial approximation
$x^{n}\approx\phi_{k}(x)$ using a truncated Chebyshev polynomial expansion.
The error in the truncated representation equals the sum of the discarded
coefficients and is precisely $p_{n,k}$. The polynomial $\phi_{k}$ and the
resulting error can both be computed analytically and, therefore, is of great
practical use. We briefly review Chebyshev polynomials in Section 2 and state
and prove the main result in Section 3. In Section 4, we explore probabilistic
interpretations of $p_{n,k}$ and obtain bounds for partial sums of binomial
coefficients.

##  2 Chebyshev polynomials

The Chebyshev polynomials of the first kind $T_{n}(x)$ for $n=0,1,2,\dots$ can
be represented as

$T_{n}(x)=\cos(n\arccos{x})\qquad x\in[-1,1].$

Starting with $T_{0}(x)=1$, $T_{1}(x)=x$, the Chebyshev polynomials satisfy a
recurrence relationship of the form $T_{n+1}(x)=2xT_{n}(x)-T_{n-1}(x)$ for
$n\geq 1$. The Chebyshev polynomials are orthogonal with respect to the
weighted inner product

$\langle u,v\rangle=\int_{-1}^{1}w(x)u(x)v(x)dx$

where the weight function takes the form $w(x)=(1-x^{2})^{-1/2}$. Any function
$f\in C[-1,1]$ that is Lipschitz continuous can be represented in terms of a
Chebyshev polynomial expansion of the form

$f(x)=\frac{1}{2}c_{0}+\sum_{j=1}^{\infty}c_{j}T_{j}(x)\qquad x\in[-1,1],$

where the coefficients $c_{j}$ are obtained as $c_{j}=\frac{2}{{\pi}}\langle
f(x),T_{j}(x)\rangle$ and the series is uniformly convergent. The monomial
$x^{n}$ is rather special since it has the following exact representation in
terms of the Chebyshev polynomials [2, Section 4]

$x^{n}=\sum_{j=0}^{n}{}^{{}^{\prime}}c_{j}T_{j}(x),$ (2)

where ${}^{{}^{\prime}}$ means the summand corresponding to $j=0$ is halved
(if it appears) and the coefficients $c_{j}$ for $j=0,\dots,n$ are

$c_{j}=\left\\{\begin{array}[]{ll}2^{1-n}\binom{n}{(n-j)/2}&n-j\text{ even}\\\
0&\text{otherwise}.\end{array}\right.$ (3)

Equation (2) takes a more familiar form, when we consider the trigonometric
perspective of Chebyshev polynomials. For example, the well-known
trigonometric identity $\cos(3\theta)=4\cos^{3}\theta-3\cos\theta$, can be
arranged as

$\cos^{3}\theta=\frac{3}{4}\cos\theta+\frac{1}{4}\cos(3\theta)=\frac{1}{2^{2}}\left(\binom{3}{1}\cos\theta+\binom{3}{0}\cos(3\theta)\right).$

With $x=\cos\theta$, we get
$x^{3}=2^{-2}\binom{3}{1}T_{1}(x)+2^{-2}\binom{3}{0}T_{3}(x)$. For
completeness, we provide a derivation of (2) in Appendix A. It is important to
note here that the series in (2) is finite, but can be truncated to obtain an
accurate approximation; see Section 3.

Chebyshev polynomials have many applications in approximation theory and
numerical analysis [9] but we limit ourselves to two such examples here.
First, if the function is differentiable $r$ times or analytic, the Chebyshev
coefficients exhibit decay (algebraic or geometric respectively). Therefore,
the Chebyshev series can be truncated to obtain an polynomial approximation of
the function and the accuracy of the approximation depends on the rate of
decay of the coefficients. Another application of Chebyshev polynomials is in
the theory and practice of polynomial interpolation. The polynomial
${q}_{n-1}^{*}(x)=x^{n}-2^{1-n}T_{n}(x)$ solves the minimax problem

$\min_{q\in\mathcal{P}_{n-1}}\|x^{n}-q(x)\|_{\infty}=2^{1-n}.$ (4)

Based on the minimax characterization, to interpolate a function over $[-1,1]$
by a polynomial of degree $n-1$, the function to be interpolated should be
evaluated at the roots of the Chebyshev polynomial $T_{n}$ given by the points
$x_{j}=\cos\left(\frac{2j+1}{2n}\pi\right)$ for $j=0,\dots,n-1$.

##  3 Main result

We construct the polynomial approximation $x^{n}\approx\phi_{k}(x)$ by
truncating the Chebyshev polynomial expansion in (2) beyond the term $j=k$.
That is

$\phi_{k}(x):=\sum_{j=0}^{k}{}^{{}^{\prime}}c_{j}T_{j}(x).$

Our main result is the following theorem, which quantifies the error in the
polynomial approximation. The proof of this theorem is based on the expression
in (2). We believe this result is new.

######  Theorem 1.

The error in the polynomial approximation $\phi_{k}(x)$ satisfies

$\|x^{n}-\phi_{k}(x)\|_{\infty}=p_{n,k}.$

###### Proof.

From (2), $x^{n}-\phi_{k}(x)=\sum_{j=k+1}^{n}c_{j}T_{j}(x)$. Using triangle
inequality, we find that
$\|x^{n}-\phi_{k}(x)\|_{\infty}\leq\sum_{j=k+1}^{n}c_{j}$ since the
coefficients are nonnegative and the Chebyshev polynomials are bounded as
$|T_{j}(x)|\leq 1$. Substituting the coefficients $c_{j}$ from (3), to get

$\|x^{n}-\phi_{k}(x)\|_{\infty}\leq\frac{1}{2^{n-1}}\sum_{\begin{subarray}{c}j=k+1\\\
n-j\text{ even}\end{subarray}}^{n}\binom{n}{(n-j)/2}.$ (5)

Using the properties of the binomial coefficients, the summation simplifies as

$\sum_{\begin{subarray}{c}j=k+1\\\ n-j\text{
even}\end{subarray}}^{n}\binom{n}{(n-j)/2}=\sum_{\begin{subarray}{c}j=k+1\\\
n+j\text{
even}\end{subarray}}^{n}\binom{n}{(n+j)/2}=\sum_{j=\lfloor(n+k)/2\rfloor+1}^{n}\binom{n}{j}.$

Plug this identity into (5) to get $\|x^{n}-\phi_{k}(x)\|_{\infty}\leq
p_{n,k}$. The bound is clearly achieved at $x=1$, where all the Chebyshev
polynomials take the value $1$. ∎

This theorem shows that the polynomial approximation $\phi_{k}$ is nearly
optimal, and the error due to this approximation is $p_{n,k}$. However, it is
the optimal polynomial for the special case $k=n-1$. It is easy to see that
$x^{n}-\phi_{n-1}(x)=2^{1-n}T_{n}(x)$ and so $\phi_{n-1}$ is the same as the
best polynomial approximation $q^{*}_{n-1}$ in (4). For $k<n-1$, from (1) and
Theorem 1

$\frac{\|x^{n}-\phi_{k}(x)\|_{\infty}}{4e}\leq\|x^{n}-\pi_{k}^{*}(x)\|_{\infty}\leq\|x^{n}-\phi_{k}(x)\|_{\infty},$

so that the error in the Chebyshev polynomial approximation is suboptimal by
at most the factor $4e\approx 10.87$. Therefore, by using $\phi_{k}$ we lose
only one significant digit of accuracy compared to $\pi^{*}_{k}$.

##  4 A probabilistic digression

In Section 1, we saw that the error in the monomial approximation depends on
$p_{n,k}$. Since $p_{n,k}$ depends on the sum of binomial coefficients, it has
a probabilistic interpretation. Newman and Rivlin [7] observed that if a fair
coin is tossed $n$ times, $p_{n,k}$ is the probability that the magnitude of
the difference between the number of heads and the number of tails exceeds
$k$. They used this insight along with the de Moivre-Laplace theorem [10,
Section 1.3] (which is a special case of the Central Limit Theorem) to obtain
the approximation $p_{n,k}\approx 2\text{erfc}(k/\sqrt{n})$.

To convert this into a rigorous inequality for $p_{n,k}$ we use a different
tool from probability, namely, concentration inequalities. The inequalities
are useful in quantifying how much a random variable deviates from its mean.
We start with the following alternative interpretation for $p_{n,k}$: it is
twice the probability that greater than $\lfloor(n+k)/2\rfloor$ coin tosses
result in heads (or equivalently tails). We associate each coin toss with an
independent Bernoulli random variable $X_{i}$ with parameter $p=1/2$ since the
coin is fair. The random variable $X=\sum_{i=1}^{n}X_{i}$ has the Binomial
distribution with parameters $n$ and $p$. Then,

$p_{n,k}=2\mathbb{P}\left(\lfloor(n+k)/2\rfloor+1\leq X\leq n\right)\leq
2\mathbb{P}\left(X\geq(n+k)/2\right).$

Since $X$ has the Binomial distribution, we can once again use the de Moivre-
Laplace theorem, to say that as $n\rightarrow\infty$,

$\frac{X-np}{\sqrt{np(1-p)}}\longrightarrow\mathcal{N}(0,1),\qquad\text{in
distribution}.$

Roughly speaking, this theorem says that $X$ behaves as a normal random
variable with mean $n/2$ and variance $n/4$. Since the tails of normal
distributions decay exponentially, we expect that $X$ lies in the range
$\frac{n}{2}\pm 1.96\sqrt{\frac{n}{4}}$ with nearly $95\%$ probability;
alternatively, the probability that it is outside this range is very small. To
make this more precise, we apply Hoeffding’s concentration inequality [10,
Theorem 2.2.6], to obtain

$\mathbb{P}\left(X\geq(n+k)/2\right)=\mathbb{P}\left(X-\mathbb{E}[X]\geq
k/2\right)\leq\exp\left(-\frac{k^{2}}{2n}\right).$

This gives our desired bound $p_{n,k}\leq 2\exp(-{k^{2}}/{2n})$.

We can use a similar technique to prove the following result which may be of
independent interest. If $0\leq k\leq n/2$, then

$\sum_{j=0}^{k}\binom{n}{j}\leq 2^{n}\exp\left(-\frac{(n-2k)^{2}}{2n}\right).$

Other concentration inequalities such as Chernoff and Bernstein (see [10,
Chapter 2]) also give equally interesting bounds. We invite the reader to
explore such results.

##  5 Acknowledgements

The author would like to thank Alen Alexanderian, Ethan Dudley, Ivy Huang,
Ilse Ipsen, and Nathan Reading for comments and feedback. The work was
supported by the National Science Foundation through the grants DMS-1745654
and DMS-1845406.

##  6 Declaration of Interest

The author has no relevant financial or non-financial competing interests to
report.

##  Appendix A Derivation of the monomial expansion

In this appendix, we provide a short derivation of (2). We take $x=\cos\theta$
and write

$\cos^{n}\theta=\left(\frac{e^{i\theta}+e^{-i\theta}}{2}\right)=\frac{1}{2^{n}}\sum_{j=0}^{n}\binom{n}{j}e^{ij\theta}e^{-i(n-j)\theta}=\frac{1}{2^{n}}\sum_{j=0}^{n}\binom{n}{j}e^{i(2j-n)\theta}.$

We have used Euler’s formula and the binomial theorem. At this point, the
derivation splits into two different paths:

1\. $n$ is odd

$\displaystyle\cos^{n}\theta=$
$\displaystyle\>\frac{1}{2^{n}}\left(\sum_{j=0}^{(n-1)/2}\binom{n}{j}e^{i(2j-n)\theta}+\sum_{j=(n+1)/2}^{n}\binom{n}{j}e^{i(2j-n)\theta}\right)$
$\displaystyle=$
$\displaystyle\>\frac{1}{2^{n-1}}\sum_{j=0}^{(n-1)/2}\binom{n}{j}\cos((n-2j)\theta)$
$\displaystyle=$
$\displaystyle\>\frac{1}{2^{n-1}}\sum_{\begin{subarray}{c}j=0\\\ n-j\text{
even}\end{subarray}}^{n}\binom{n}{(n-j)/2}\cos(j\theta).$

2\. $n$ is even

$\displaystyle\cos^{n}\theta=$
$\displaystyle\>\frac{1}{2^{n}}\left(\sum_{j=0}^{n/2-1}\binom{n}{j}e^{i(2j-n)\theta}+\binom{n}{n/2}+\sum_{j=n/2+1}^{n}\binom{n}{j}e^{i(2j-n)\theta}\right)$
$\displaystyle=$
$\displaystyle\>\frac{1}{2^{n}}\binom{n}{n/2}+\frac{1}{2^{n-1}}\sum_{j=0}^{n/2-1}\binom{n}{j}\cos((n-2j)\theta)$
$\displaystyle=$
$\displaystyle\>\frac{1}{2^{n-1}}\sum_{\begin{subarray}{c}j=0\\\ n-j\text{
even}\end{subarray}}^{n}{}^{{}^{\prime}}\binom{n}{(n-j)/2}\cos(j\theta).$

In either case, substitute $x=\cos\theta$ and $T_{j}(x)=\cos(j\theta)$ to
complete the derivation.

## References

  * [1] H. Avron and S. Toledo.  Randomized algorithms for estimating the trace of an implicit symmetric positive semi-definite matrix.  Journal of the ACM (JACM), 58(2):8, 2011.
  * [2] W. J. Cody.  A survey of practical rational and polynomial approximation of functions.  SIAM Review, 12(3):400–423, 1970.
  * [3] T. A. Driscoll, N. Hale, and L. N. Trefethen.  Chebfun Guide.  Pafnuty Publications, 2014.
  * [4] E. Dudley, A. K. Saibaba, and A. Alexanderian.  Monte Carlo estimators for the Schatten p-norm of symmetric positive semidefinite matrices.  arXiv preprint arXiv:2005.10174, 2020.
  * [5] C. Lanczos.  Applied analysis.  Dover Books on Advanced Mathematics. Dover Publications, Inc., New York, 1988.  Reprint of the 1956 original.
  * [6] Y. Nakatsukasa and L. Trefethen.  Rational approximation of $x^{n}$.  Proceedings of the American Mathematical Society, 146(12):5219–5224, 2018.
  * [7] D. J. Newman and T. J. Rivlin.  Approximation of monomials by lower degree polynomials.  Aequationes Mathematicae, 14(3):451–455, 1976.
  * [8] A. R. Reddy.  Approximations to $x^{n}$ and $|x|$—a survey.  J. Approx. Theory, 51(2):127–137, 1987.
  * [9] L. N. Trefethen.  Approximation theory and approximation practice.  Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2013.
  * [10] R. Vershynin.  High-dimensional probability, volume 47 of Cambridge Series in Statistical and Probabilistic Mathematics.  Cambridge University Press, Cambridge, 2018.