# Floquet superradiance lattices in thermal atoms

Xingqi Xu1,∗, Jiefei Wang1,2∗, Jianhao Dai1, Ruosong Mao1, Han Cai2,1†, Shi-
Yao Zhu1,3, and Da-Wei Wang1,3,4‡ 1Interdisciplinary Center for 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, Zhejiang Province, China;
2College of Optical Science and Engineering, Zhejiang University, Hangzhou,
310027, China;
3Hefei National Laboratory, Hefei 230088, China;
4CAS Center for Excellence in Topological Quantum Computation, University of
Chinese Academy of Sciences, Beijing 100190, China

###### Abstract

Floquet modulation has been widely used in optical lattices for coherent
control of quantum gases, in particular for synthesizing artificial gauge
fields and simulating topological matters. However, such modulation induces
heating which can overwhelm the signal of quantum dynamics in ultracold atoms.
Here we report that the thermal motion, instead of being a noise source,
provides a new control knob in Floquet-modulated superradiance lattices, which
are momentum-space tight-binding lattices of collectively excited states of
atoms. The Doppler shifts combined with Floquet modulation provide effective
forces along arbitrary directions in a lattice in frequency and momentum
dimensions. Dynamic localization, dynamic delocalization and chiral edge
currents can be simultaneously observed from a single transport spectrum of
superradiance lattices in thermal atoms. Our work paves a way for simulating
Floquet topological matters in room-temperature atoms and facilitates their
applications in photonic devices.

Collective interaction between atoms and light plays an important role in
quantum optics Gross1982 , Inouye1999 . A most famous example is the Dicke
superradiance featured by enhanced decay rates of collectively excited atoms
Dicke1954 , in particular for atoms confined in a small volume compared with
the light wavelength. In the opposite limit, i.e., when the size of the atomic
cloud is much larger than the light wavelength, Scully et al. Scully2006 ,
Scully2009  found that the light momentum can be stored in atoms prepared in
the so-called timed Dicke states (TDSs) He2020 , Araujo2016 , Roof2016 . The
stored momentum can be released in the directional superradiant emission
Scully2006 . More interestingly, such TDSs with different momenta can be
coupled by multiple lasers to form momentum-space superradiance lattices (SLs)
Wang2015 , Wang2015optica , which provide a new platform for quantum
simulation.

In contrast to real-space optical lattices that require ultracold temperature,
SLs can be implemented in both cold atoms Chen2018 , Wang2020npj , Mi2021  and
room-temperature atoms Cai2019 , He2021 , Mao2022 , since TDSs are robust to
the center-of-mass motion of atoms Celardo2014 , Damanet2016 , Bromley2016 ,
Weiss2019 . Such a robustness can be understood from a distinct feature of the
momentum-space SL, i.e., its Brillouin zone (BZ) is in real space. Here the
positions of atoms play the same role of the crystal momenta of electrons in
solids. Atoms moving through the real-space BZs follow the same dynamics of
electrons in a DC electric field. Therefore, the atomic motion, instead of
being a source of noise, can simulate an effective electric field in momentum
space, which has been used to measure the geometric Zak phase Mao2022 ,
Zak1989 , Xiao2010 .

Here we show that SLs can also be used to study the rich physics in Floquet
(temporally periodic) systems Eckardt2017 , Oka2019 , Rudner2020 ,
Weitenberg2021 , Torre2021 , where the thermal motion induced DC field
significantly enriches the tunability. Being a widely-used technique in
quantum engineering Oka2009 , Lindner2011 , Rechtsman2013 , Jotzu2014 ,
Rudner2013 , Mukherjee2017 , Cheng2019 , Jiang2011 , Titum2016 ,
Kolodrubetz2018 , Nathan2021 , Li2022 , Floquet driving Floquet , Shirley1965
, Leon2013 , Goldman2014  can extend the 1D SLs into (1+1)D (momentum-
frequency dimensional) lattices, allowing us to simulate higher-dimensional
lattice dynamics Martin2017 , Boyers2020 , Long2021 , Malz2021 . In the (1+1)D
SLs, there is a constant effective electric field in the frequency dimension,
determined by the Floquet modulation frequency Dutt2021 . On the other hand,
the atomic velocity provides a DC field along the momentum dimension with a
strength proportional to the Doppler shift. By resolving the direction of the
total field with a spectroscopic method, we observe dynamic localization,
dynamic delocalization and chiral edge currents Dunlap1986 , Holthaus1992 ,
Madison1998 , Longhi2006 , Lignier2007 , Sias2008 , Szameit2009 , Szameit2010
, Haller2010 , Ma2011 , Chen2011 , Yuan2015 , Mancini2015 , Stuhl2015 ,
Libi2016 , Atala2014 , An2017 , Wang2020  in SLs of thermal atoms.

Figure 1: (a) Schematic configuration of the experimental setup. The coupling
fields are phased-modulated by two EOMs at a RF frequency $\delta_{2}=80$ MHz.
Insets: the atomic levels in the reference frame of the atoms with velocity
$v_{x}$. The coupling fields have opposite Doppler shifts $\pm\delta_{1}$ and
periodically modulated phases $\theta_{\pm}(t)$. (b) Mapping between the
dynamically modulated 1D SL and the (1+1)D Floquet SL, with white and gray
circles labeling the $b$\- and $a$-sublattice sites, respectively. Left: the
driven 1D momentum-space lattices with phase-modulated hopping strengths. The
yellow arrow $\bm{F}_{1}$ indicates the direction of atomic motion induced DC
electric force. Right: the equivalent (1+1)D frequency-momentum lattices. The
yellow arrow $\bm{F}$ indicates the total DC electric force with components
$\bm{F}_{1}$ and $\bm{F}_{2}$ along the momentum and frequency dimensions. The
probed and detected sites are marked schematically.

Experimental setup and theoretical model. The experiment is implemented in a
vapor cell of natural abundance rubidium atoms in a standing-wave-coupled
electromagnetically induced transparency (EIT) configuration (see Fig. 1(a)).
A weak probe field couples the ground state
$|g\rangle\equiv|5^{2}S_{1/2},F=1\rangle$ to an excited state
$|b\rangle\equiv|5^{2}P_{1/2},F=2\rangle$ in 87Rb D1 line. The excited state
is resonantly coupled to a metastable state
$|a\rangle\equiv|5^{2}S_{1/2},F=2\rangle$ by two strong laser fields. The
probe field pumps the atoms to the TDS,

$|b_{\bm{k}_{p}}\rangle=\frac{1}{\sqrt{N}}\sum_{j=1}^{N}e^{i\bm{k}_{p}\cdot\bm{r}_{j}}|g_{1}...b_{j}...g_{N}\rangle,$
(1)

where $\bm{k}_{p}$ is the wavevector of the probe field, $\bm{r}_{j}$ is the
position of the $j$th atom. The standing-wave coupling fields couple
$|b_{\bm{k}_{p}}\rangle$ to other TDSs with discrete momenta to form SLs
Wang2015 , in which the directional superradiant emission from the TDS
$|b_{-\bm{k}_{p}}\rangle$ is collected along $-\bm{k}_{p}$ direction
Scully2006 , Scully2009 , Roof2016 , Araujo2016 , He2020 , such that the
lattice transport from $|b_{\bm{k}_{p}}\rangle$ to $|b_{-\bm{k}_{p}}\rangle$
can be investigated. We can tune the probe frequency to measure the transport
at different energies, similar to tuning the Fermi energy of electrons in
solids.

In our experiment on 1D SLs, the probe field propagates in $+x$ direction with
a wave vector $k_{p}$, parallel to the two coupling laser fields with wave
vectors $\pm k_{c}$ and $k_{p}\approx k_{c}$. The phases of the coupling
fields are modulated by two electric-optic modulators (EOMs) independently,
$\theta_{\pm}(t)=f\sin(\delta_{2}t\pm\phi/2)$, where $f$ and $\delta_{2}$ are
the modulation depth and frequency, and $\phi$ controls the phase delay
between the two radio-frequency (RF) driving fields. The interaction between
the atoms and the coupling fields is described by a Floquet SL (see Fig.
1(b)). The EOMs introduce periodic modulation of the coupling strengths in the
SL Hamiltonian (we set $\hbar=1$ and see details in SM ),

$H_{d}=\sum_{j}\Omega a^{\dagger}_{2j}[e^{-i\theta_{+}(t)}b^{\
}_{2j+1}+e^{-i\theta_{-}(t)}b^{\ }_{2j-1}]+H.c.,$ (2)

where $\Omega$ is the Rabi frequency of the coupling fields,
$a_{2j}^{\dagger}$ and $b_{2j+1}^{\dagger}$ are creation operators in the
$j$th unit cell of the SL, $a_{2j}^{\dagger}\equiv
1/\sqrt{N}\sum_{l}e^{2ijk_{c}x_{l}}|a_{l}\rangle\langle g_{l}|$ and
$b_{2j+1}^{\dagger}\equiv
1/\sqrt{N}\sum_{l}e^{i(2j+1)k_{c}x_{l}}|b_{l}\rangle\langle g_{l}|$ with
$x_{l}$ being the position of the $l$th atom and $N$ the total number of
atoms. They create TDSs from the ground state, e.g.,
$a_{2j}^{\dagger}|g_{1},g_{2},...,g_{N}\rangle=1/\sqrt{N}\sum_{l}\exp(2ijk_{c}x_{l})|g_{1}...a_{l}...g_{N}\rangle$,
which carries a momentum $2jk_{c}$.

Using the Floquet theorem for $H_{d}$, we obtain an effective (1+1)D
Hamiltonian $V_{2}+H_{s}$ SM , including a potential energy in the frequency
dimension Yu2021 ,

$V_{2}=-\delta_{2}\sum_{j,m}m(a^{\dagger}_{2j,m}a^{\
}_{2j,m}+b^{\dagger}_{2j+1,m}b^{\ }_{2j+1,m}),$ (3)

and the static lattice Hamiltonian $H_{s}=\sum_{n}H_{s}^{[n]}$ with

$\begin{split}H_{s}^{[n]}=\sum_{j,m}&\Omega
e^{in\phi/2}a^{\dagger}_{2j,m}[J_{-n}(f)b_{2j+1,m+n}\\\
&+J_{n}(f)b_{2j-1,m-n}]+H.c.,\end{split}$ (4)

which couples lattice sites along the crystal direction
$\hat{e}_{1}+n\hat{e}_{2}$ (denoted by Miller index $[1,n]$) with
$\hat{e}_{1}$ and $\hat{e}_{2}$ being the unit vectors along the momentum and
frequency dimensions. Here $J_{n}$ is the $n$th order Bessel function and
$d_{j,m}^{\dagger}$ ($d=a,b$) is the creation operator for
$|d_{j,m}\rangle\equiv e^{-im\delta_{2}t}|d_{j}\rangle$ Leon2013 , SM . We can
understand $|d_{j,m}\rangle$ as a replica of $|d_{j}\rangle$ that carries
$-m\delta_{2}$ energy from the periodic driving Gao2021 . The $n$th order
hopping coefficient carries a phase $\pm n\phi/2$, which can synthesize a
gauge field in the lattice (see Fig. 1 (b)). The linear potential along the
frequency dimension $V_{2}$ can be characterized by a force
$\bm{F}_{2}=\delta_{2}\hat{e}_{2}$. Throughout the experiment we set
$\delta_{2}=80$ MHz and $\Omega=25$ MHz to guarantee that the driving
frequency is much larger than hopping strengths.

For moving atoms, the photons in the two coupling fields along $\pm x$
directions have extra energies due to the Doppler shifts $\pm\delta_{1}$ with
$\delta_{1}=-k_{c}v_{x}$ and $v_{x}$ being the atomic velocity in $x$
direction. Since the TDSs in the SL are created by absorbing and emitting
photons in the coupling fields, the SL acquires a linear potential along the
momentum dimension (in the lab reference frame),

$V_{1}=-\delta_{1}\sum_{j,m}[2ja_{2j,m}^{\dagger}a_{2j,m}+(2j+1)b_{2j+1,m}^{\dagger}b_{2j+1,m}],$
(5)

which has the effect of a force $\bm{F}_{1}=\delta_{1}\hat{e}_{1}$.
Consequently, the force associated to the total linear potential $V_{1}+V_{2}$
is $\bm{F}=\bm{F}_{1}+\bm{F}_{2}=\delta_{1}\hat{e}_{1}+\delta_{2}\hat{e}_{2}$.
In our experiment, the hopping of excitations along the force direction is
inhibited because the potential difference is much larger than the coupling
strengths.

Figure 2: Correspondence between peaks at $\Delta_{p}=n\delta_{2}$ in the
reflection spectrum and the dynamics along the crystal direction $[1,n]$
(dashed lines). The $n$th peak is contributed by the atoms in the velocity
group with a Doppler shift $\delta_{1}=-n\delta_{2}$. The effective force
$\bm{F}=-\Delta_{p}\hat{e}_{1}+\delta_{2}\hat{e}_{2}$ makes the dynamics of
(1+1)D lattices dominated by transitions along the crystal direction $[1,n]$.
The reflection signal shows the transport from the probed state
$|b_{1,0}\rangle$ to the detected site $|b_{-1,-2n}\rangle$. The phase delay
$\phi=\pi$ and modulation depth $f=3$.

Lattice transport measurement. We can tune the probe frequency to measure the
transport in SLs with different values of $\bm{F}$. The interaction
Hamiltonian between the probe field and the atoms is
$H_{p}=\sqrt{N}\Omega_{p}e^{-i\Delta_{p}t}b^{\dagger}_{1,0}+H.c.$, where
$\Delta_{p}=\nu_{p}-\omega_{bg}$ with $\nu_{p}$ being the probe field
frequency and $\omega_{bg}$ being the atomic transition frequency in the lab
reference frame. The probe field couples the ground state to a TDS
$|b_{1,0}\rangle$, which has a potential energy $-\delta_{1}$ proportional to
the atomic velocity, such that we need to detune the probe field
$\Delta_{p}\approx-\delta_{1}$ to excite the TDSs of atoms with the
corresponding velocities. Since $\bm{F}_{1}$ is also proportional to the
atomic velocity, we can measure lattice dynamics with different $\bm{F}_{1}$
at the corresponding probe detunings. In another word, with the probe detuning
$\Delta_{p}$, we measure the transport from $|b_{1,0}\rangle$ to
$|b_{-1,m}\rangle$ in an effective force
$\bm{F}=-\Delta_{p}\hat{e}_{1}+\delta_{2}\hat{e}_{2}$. Such transport is
experimentally measured from the superradiant directional emission of
$|b_{-1,m}\rangle$ along the $-x$ direction, which satisfies the phase-
matching condition $-2k_{c}+k_{p}\approx-k_{p}$ SM . Therefore, we can
simultaneously measure the dynamics with different effective forces in a
single transport spectrum by detuning the probe field frequency. The coherent
probe field may also create more than one excitations which, however, have the
same dynamic response as the TDSs due to the bosonic nature of atomic
excitations in a weak probe field Wang2015 , Cai2019 .

Figure 3: Dynamic localization and delocalization. (a) Experimental reflection
spectra as a function of the probe detuning $\Delta_{p}$ and the modulation
depth $f$. We denote the points for Stark localization, dynamic localization,
and dynamic delocalization with dots in blue, red and green colors. (b) The
relfections at $\Delta_{p}=n\delta_{2}$ characterizing transport along crystal
directions $[1,n]$ ($n=0,1,2,3,4$) are measured as functions of $f$, where the
transport is completely inhibited at the zero points (marked by black arrow)
of $J_{n}(f)$. We set the phase delay $\phi=\pi$.

We show a typical spectrum for the transport from $|b_{1,0}\rangle$ to
$|b_{-1,m}\rangle$ in Fig. 2. Due to the Stark localization, the transport is
only significant when the two states have the same potential energy, i.e.,
when the force points in a direction perpendicular to a crystal direction of
the (1+1)D SL. As a result, the reflection is featured by discrete peaks at
$\Delta_{p}=n\delta_{2}$. The reflection peak at $\Delta_{p}=n\delta_{2}$ is
contributed by the transport from $|b_{1},0\rangle$ to $|b_{-1,-2n}\rangle$
along the crystal direction $[1,n]$ (dotted lines in the insets),
perpendicular to $\bm{F}$. From the reflection spectra we are able to
investigate the dynamic localization and delocalization by simultaneously
tuning the modulation depths of the two EOMs while keep their phase delay
$\phi=\pi$ fixed. We can also tune $\phi$ to introduce effective magnetic
fluxes in the (1+1)D SLs to investigate the Floquet chiral edge currents.

Figure 4: Floquet chiral currents. (a) The schematics of the effective
Hamiltonian along a composite dimension of momentum and frequency. The NNN
hoppings (green lines with arrows indicating the transition directions
attached with a phase $\pi/2$) along crystal direction $[1,n]$ ($n=0,1,2$) are
induced by all possible second-order transition paths (insets). (b) The zigzag
chiral ladder along $[1,n]$ (replotted from the reduced 1D lattices in (a))
with an effective magnetic flux $\Phi_{n}$ in a unit loop. The two sublattices
are the two edges which host chiral edge currents with directions (along thick
arrows) determined by $\Phi_{n}$. The green arrows denote the transition
directions attached with a phase $\pi/2$. (c) Phase diagram of $\Phi_{n}$ as
functions of $\phi$ and $f$. (d) The phase diagram of chiral currents measured
by $\delta R_{n}\equiv R(n\delta_{1})-R^{\prime}(-n\delta_{1})$, which are
consistent with the corresponding phase diagram of $\Phi_{n}$ with phase
boundaries indicated by solid lines.

Dynamic localization and delocalization. We measure the reflection spectra as
a function of the probe detuning $\Delta_{p}$ and the modulation depth $f$
(see Fig. 3(a)). The transport along the crystal direction $[1,n]$ is
demonstrated near the detuning $\Delta_{p}=n\delta_{2}$ and governed by the
interaction Hamiltonian $H_{s}^{[n]}$ in Eq. (4). Remarkably, the transport is
totally inhibited when the modulation depth $f$ reaches to the zeros of the
corresponding Bessel functions $J_{n}(f)$ (marked by black arrows in Fig.
3(b)), where the effective hopping strength is zero in $H_{s}^{[n]}$. In
particular, at the detuning $\Delta_{p}=0$, the reflectivity vanishes at
$f=2.4$, $5.5$, and $8.7$, which are the first three zeros of $J_{0}(f)$,
characterizing the dynamic localization Dunlap1986 , Holthaus1992 . For
parameters in our experiment, localization at zeros of $J_{n}(f)$ can be
clearly identified up to $|n|=4$, a significant improvement compared with
other platforms Sias2008 , Szameit2009 . The order $n$ is limited by the width
of the Doppler broadening, which has a FWHM of $500$ MHz at the experiment
temperature $60^{\circ}$C.

The measured transport can be understood as the interplay between the lattice
dynamics in AC and DC field Dunlap1986 , Holthaus1992 , Madison1998 ,
Longhi2006 , Lignier2007 , Sias2008 , Szameit2010 , Haller2010 , Ma2011 ,
Chen2011 , Yuan2015  in a 1D SL. The oscillating phases $\theta_{\pm}(t)$
introduce a uniform AC field $(\partial\theta_{-}/\partial
t-\partial\theta_{+}/\partial t)/2=f\delta_{2}\sin(\delta_{2}t)$ Leon2013 . We
observe three important effects in the presence of both AC and DC field as
being marked in Fig. 3(a). Without AC and DC fields, i.e., for
$\delta_{1}=f=0$ near $\Delta_{p}=0$, the 1D SL transport signal is strong
from $|b_{1}\rangle$ to $|b_{-1}\rangle$. When we keep $f=0$ and detuning
$\Delta_{p}$ away from zero, the signal disappears due to Stark localization
in a DC field (see the blue point at $\Delta_{p}=\delta_{2}/2$ and $f=0$).
When we keep $\Delta_{p}=0$ and increase $f$, an effective AC field leads to
the dynamic localization at $f=2.4$ Dunlap1986 , Holthaus1992 , Madison1998 ,
Longhi2006 , Lignier2007 , Szameit2010 . With both AC and DC fields, the
transport is recovered when the modulation frequency of the AC field is on
resonance with the Bloch frequency of the DC field (see the green dot for
$\Delta_{p}=\delta_{2}$), which is called dynamic delocalization or photon-
assisted tunnelling Sias2008 , Haller2010 , Ma2011 , Chen2011 , Yuan2015 .

Phase diagram of chiral currents. Chiral edge currents have been observed in
zigzag SLs Cai2019 . Floquet engineering can bring new control knobs to induce
such unidirectional transport Mancini2015 , Stuhl2015 , Libi2016 , Atala2014 ,
An2017 , Wang2020 . Near $\Delta_{p}=n\delta_{2}$, the Hamiltonian
$H_{s}^{[n]}$ introduces equipotential transitions while $H_{s}^{[m]}$’s with
$m\neq n$ provide second-order transitions via intermediate states (see green
lines in Fig. 4(a) and derivation in SM ),

$\displaystyle K^{[n]}$
$\displaystyle=\sum_{j,m}i\kappa_{n}(b^{\dagger}_{2j-1,m-n}b^{\ }_{2j+1,m+n}$
(6) $\displaystyle-a^{\dagger}_{2j,m}a^{\ }_{2j+2,m-2n})+H.c.,$

with
$\kappa_{n}=-\sum^{\infty}_{p=1}{2\sin(p\phi)}\Omega^{2}J_{p+n}(f)J_{p-n}(f)/{p\delta_{2}}$
SM . When $\phi=\pi$, we obtain $\kappa_{n}=0$ for dynamic localization and
delocalization. In general, for $\phi\neq 0$ and $\pi$, $H_{s}^{[n]}$ and
$K^{[n]}$ induce the nearest-neighbor and next-nearest-neighbor (NNN) hoppings
(see red lines and green arrows in Fig. 4(a) and (b)). The phase factors in
$K^{[n]}$ bring a magnetic flux in a unit closed loop (circular arrows in Fig.
4(b)) Cai2019 , Hugel2014 ,

$\Phi_{n}=\frac{\pi}{2}(-1)^{n}\text{sgn}[\kappa_{n}],$ (7)

which depends on both $f$ and $\phi$ and exhibits rich structures along
different crystal directions (see Fig. 4(c)).

The direction of the synthesized magnetic field, i.e., the sign of $\Phi_{n}$,
determines the unidirectional transport along $b$-sublattice, regardless of
the strength $\kappa_{n}$ Cai2019 . The chiral current flows to the right and
left for $\Phi_{n}=\pi/2$ and $-\pi/2$, which results in asymmetric
distribution of steady-state probabilities in the SL with respect to the
probed site. Such an asymmetry can be measured by comparing the reflection
signals of two probe fields Cai2019 , He2022  in opposite directions and with
opposite detunings. We define $R$ and $R^{\prime}$ as the reflectivities of
the probe field incident along $+x$ and $-x$ directions. $R(n\delta_{2})$
characterizes the transport from the site $|b_{1,0}\rangle$ to
$|b_{-1,-2n}\rangle$ while $R^{\prime}(-n\delta_{2})$ describes the transport
from the site $|b_{-1,0}\rangle$ to $|b_{1,2n}\rangle$ SM . The chiral current
along the crystal direction $[1,n]$ are measured by $\delta R_{n}\equiv
R(n\delta_{2})-R^{\prime}(-n\delta_{2})$ (see Fig. 4(d)), which agrees with
the phase diagrams in Fig. 4(c).

In conclusion, we demonstrate that the SL provides a highly tunable platform
to study the Floquet dynamics in thermal atoms. We simultaneously observe
dynamic localization and delocalization from a single transport spectrum. We
also experimentally measure the phase diagrams of chiral edge currents in a
(1+1)D lattice with both artificial magnetic fluxes and electric fields. We
can increase hopping strengths between TDSs to overwhelm the linear
potentials, such that the motion of excitations is governed by the underlying
2D Bloch bands with the effective electric force being a perturbation. In such
a regime topological non-equilibrium phenomena such as quantized energy
pumping Martin2017 , Kolodrubetz2018 , Boyers2020 , Long2021 , Malz2021 ,
Nathan2021  and space-time crystals Morimoto2017 , Xu2018 , Peng2019 , Gao2021
, Peng2022  can be investigated. Furthermore, we can introduce interaction in
SLs by using additional lasers to couple atoms to the Rydberg states Li2020
to realize many-body phenomena beyond the short-range interaction SM , such as
the interplay between long-range nonlocal interaction and dynamic
localization-delocalization transition Cao2022 , Toh2022 .

We thank Luqi Yuan for insightful suggestions. This work was supported by the
National Natural Science Foundation of China (Grants No. U21A20437, No.
11874322 and No. 11934011), the National Key Research and Development Program
of China (Grants No. 2019YFA0308100, No. 2018YFA0307200 and No.
2017YFA0304202), Zhejiang Province Key Research and Development Program (Grant
No. 2020C01019), Innovation Program for Quantum Science and Technology (Grant
No. 2021ZD0303200), the Strategic Priority Research Program of Chinese Academy
of Sciences (Grant No. XDB28000000), and the Fundamental Research Funds for
the Central Universities.

∗These authors contributed equally to this work.

<EMAIL_ADDRESS>
<EMAIL_ADDRESS>
## References

  * [1] M. Gross and S. Haroche,  Superradiance: An essay on the theory of collective spontaneous emission,  Phys. Rep. 93, 301 (1982).
  * [2] S. Inouye, A. P. Chikkatur, D. M. Stamper-Kurn, J. Stenger, D. E. Pritchard, W. Ketterle,  Superradiant Rayleigh Scattering from a Bose-Einstein Condensate,  Science 285, 571 (1999).
  * [3] R. H. Dicke,  Coherence in spontaneous radiation processes,  Phys. Rev. 93, 99 (1954).
  * [4] M. O. Scully, E. S. Fry, C. H. Ooi, and K. Wódkiewicz,  Directed spontaneous emission from an extended ensemble of N atoms: Timing is everything,  Phys. Rev. Lett. 96, 010501 (2006).
  * [5] M. O. Scully and A. A. Svdzinsky,  The super of superradiance.  Science 325, 1510 (2009).
  * [6] M. O. Araújo, I. Krešić, R. Kaiser, and W. Guerin,  Superradiance in a Large Cloud of Cold Atoms in the Linear-Optics Regime,  Phys. Rev. Lett. 117, 073002 (2016).
  * [7] S. J. Roof, K. J. Kemp, M. D. Havey, and I. M. Sokolov,  Observation of Single-Photon Superradiance and the Cooperative Lamb Shift in an Extended Sample of Cold Atoms,  Phys. Rev. Lett. 117, 073003 (2016).
  * [8] Y. He, L. Ji, Y. Wang, L. Qiu, J. Zhao, Y. Ma, X. Huang, S. Wu, and D. E. Chang,  Geometric Control of Collective Spontaneous Emission,  Phys. Rev. Lett. 125, 213602 (2020).
  * [9] D. W. Wang, R. B. Liu, S. Y. Zhu, and M. O. Scully,  Superradiance Lattice,  Phys. Rev. Lett. 114, 043602 (2015).
  * [10] D. W. Wang, H. Cai, L. Yuan, S. Y. Zhu, and R. B. Liu,  Topological phase transitions in superradiance lattices,  Optica 2, 712 (2015).
  * [11] L. Chen, P. Wang, Z. Meng, L. Huang, H. Cai, D. W. Wang, S. Y. Zhu, and J. Zhang,  Experimental Observation of One-Dimensional Superradiance Lattices in Ultracold Atoms,  Phys. Rev. Lett. 120, 193601 (2018).
  * [12] P. Wang, L. Chen, C. Mi, Z. Meng, L. Huang, S. N. Khan, H. Cai, D. W. Wang, S. Y. Zhu, and J. Zhang,  Synthesized magnetic field of a sawtooth superradiance lattice in Bose-Einstein condensates,  npj Quantum Inf. 6, 18 (2020).
  * [13] C. Mi, K. S. Nawaz, L. Chen, P. Wang, H. Cai, D. W. Wang, S. Y. Zhu, and J. Zhang,  Time-resolved interplay between superradiant and subradiant states in superradiance lattices of Bose-Einstein condensates,  Phys. Rev. A 104, 043326 (2021).
  * [14] H. Cai, J. Liu, J. Wu, Y. He, S. Y. Zhu, J. X. Zhang, and D. W. Wang,  Experimental Observation of Momentum-Space Chiral Edge Currents in Room-Temperature Atoms,  Phys. Rev. Lett. 122, 023601 (2019).
  * [15] Y. He, R. Mao, H. Cai, J. X. Zhang, Y. Li, L. Yuan, S. Y. Zhu, and D. W. Wang,  Flat-Band Localization in Creutz Superradiance Lattices,  Phys. Rev. Lett. 126, 103601 (2021).
  * [16] R. Mao, X. Xu, J. Wang, C. Xu, G. Qian, H. Cai, S. Y. Zhu, and D. W. Wang,  Measuring Zak phase in room-temperature atoms,  Light Sci. Appl. 11, 291 (2022).
  * [17] G. L. Celardo, G. G. Giusteri, and F. Borgonovi,  Cooperative robustness to static disorder: Superradiance and localization in a nanoscale ring to model light-harvesting systems found in nature,  Phys. Rev. B 90, 075113 (2014).
  * [18] F. Damanet, D. Braun, and J. Martin,  Cooperative spontaneous emission from indistinguishable atoms in arbitrary motional quantum states,  Phys. Rev. A 94, 033838 (2016).
  * [19] S. L. Bromley, B. Zhu, M. Bishof, X. Zhang, T. Bothwell, J. Schachenmayer, T. L. Nicholson, R. Kaiser, S. F. Yelin, M. D. Lukin, A. M. Rey, and J. Ye,  Collective atomic scattering and motional effects in a dense coherent medium,  Nat. Commun. 7, 11039 (2016).
  * [20] P. Weiss, A. Cipris, M. O. Araújo, R. Kaiser, and W. Guerin,  Robustness of Dicke subradiance against thermal decoherence,  Phys. Rev. A 100, 033833 (2019).
  * [21] J. Zak,  Berry’s phase for energy bands in solids,  Phys. Rev. Lett. 62, 2747 (1989).
  * [22] D. Xiao, M. C. Chang, and Q. Niu,  Berry phase effects on electronic properties,  Rev. Mod. Phys. 82, 1959 (2010).
  * [23] A. Eckardt,  Colloquium: Atomic quantum gases in periodically driven optical lattices,  Rev. Mod. Phys. 89, 011004 (2017).
  * [24] T. Oka and S. Kitamura,  Floquet Engineering of Quantum Materials,  Annu. Rev. Condens. Matter Phys. 10,387 (2019).
  * [25] M. S. Rudner, N. H. Lindner,  Band structure engineering and non-equilibrium dynamics in Floquet topological insulators,  Nat. Rev. Phys. 2, 229 (2020).
  * [26] C. Weitenberg and J. Simonet,  Tailoring quantum gases by Floquet engineering,  Nat. Phys. 17, 1342 (2021).
  * [27] A. de la Torre, D. M. Kennes, M. Claassen, S. Gerber, J. W. McIver, and M. A. Sentef  Colloquium: Nonthermal pathways to ultrafast control in quantum materials,  Rev. Mod. Phys. 93, 041002 (2021).
  * [28] T. Oka and H. Aoki,  Photovoltaic Hall effect in graphene,  Phys. Rev. B 79, 081406(R) (2009).
  * [29] N. H. Lindner, G. Refael, and V. Galitski,  Floquet topological insulator in semiconductor quantum wells,  Nat. Phys. 7, 490 (2011).
  * [30] G. Jotzu, et al.,  Experimental realization of the topological Haldane model with ultracold fermions,  Nature 515, 237 (2014).
  * [31] M. C. Rechtsman, et al.,  Photonic Floquet topological insulators,  Nature 496, 196 (2013).
  * [32] M. S. Rudner, N. H. Lindner, E. Berg, and M. Levin  Anomalous Edge States and the Bulk-Edge Correspondence for Periodically Driven Two-Dimensional Systems,  Phys. Rev. X 3, 031005 (2013).
  * [33] S. Mukherjee, et al.,  Experimental observation of anomalous topological edge modes in a slowly driven photonic lattice,  Nat. Commun. 8, 13918 (2017).
  * [34] Q. Cheng, Y. Pan, H. Wang, C. Zhang, D. Yu, A. Gover, H. Zhang. T. Li, L. Zhou, and S. Zhu,  Observation of Anomalous $\pi$ Modes in Photonic Floquet Engineering,  Phys. Rev. Lett. 122, 173901 (2019).
  * [35] L. Jiang, T. Kitagawa, J. Alicea, A. R. Akhmerov, D. Pekker, G. Refael, J. I. Cirac, E. Demler, M. D. Lukin, and P. Zoller,  Majorana Fermions in Equilibrium and in Driven Cold-Atom Quantum Wires,  Phys. Rev. Lett. 106, 220402 (2011).
  * [36] P. Titum, E. Berg, M. S. Rudner, G. Refael, and N. H. Lindner,  Anomalous Floquet-Anderson Insulator as a Nonadiabatic Quantized Charge Pump,  Phys. Rev. X 6, 021013 (2016).
  * [37] M. H. Kolodrubetz, F. Nathan, S. Gazit, T. Morimoto, and J. E. Moore,  Topological Floquet-Thouless Energy Pump,  Phys. Rev. Lett. 120, 150601 (2018).
  * [38] F. Nathan, R. Ge, S. Gazit, M. Rudner, and M. Kolodrubetz,  Quasiperiodic Floquet-Thouless Energy Pump,  Phys. Rev. Lett. 127, 166804 (2021).
  * [39] G. Li, L. Wang, R. Ye, S. Liu, Y. Zheng, L. Yuan, and X. Chen,  Observation of flat-band and band transition in the synthetic space,  Adv. Photon. 4, 036002 (2022).
  * [40] G. Floquet,  Sur les equations differentielles lineaires a coefficients periodiques,  Ann. Ecole Norm. Superieure 12, 47 (1883).
  * [41] J. H. Shirley,  Solution of the Schrödinger Equation with a Hamiltonian Periodic in Time,  Phys. Rev. 138, B979 (1965).
  * [42] N. Goldman and J. Dalibard,  Periodically Driven Quantum Systems: Effective Hamiltonians and Engineered Gauge Fields,  Phys. Rev. X 4, 031027 (2014).
  * [43] A. Gómez-León and G. Platero,  Floquet-Bloch Theory and Topology in Periodically Driven Lattices,  Phys. Rev. Lett. 110, 200403 (2013).
  * [44] I. Martin, G. Refael, and B. Halperin,  Topological Frequency Conversion in Strongly Driven Quantum Systems,  Phys. Rev. X 7, 041008 (2017).
  * [45] E. Boyers, P. J. D. Crowley, A. Chandran, and A. Sushkov,  Exploring 2D Synthetic Quantum Hall Physics with a Quasiperiodically Driven Qubit,  Phys. Rev. Lett. 125, 160505 (2020).
  * [46] D. M. Long, P. J. D. Crowley, and A. Chandran,  Nonadiabatic Topological Energy Pumps with Quasiperiodic Driving,  Phys. Rev. Lett. 126, 106805 (2021).
  * [47] D. Malz and A. Smith,  Topological Two-Dimensional Floquet Lattice on a Single Superconducting Qubit,  Phys. Rev. Lett. 126, 163602 (2021).
  * [48] A. Dutt, Q. Lin, L. Yuan, M. Minkov, M. Xiao, S. Fan  A single photonic cavity with two independent physical synthetic dimensions,  Science 367, 59-64 (2020).
  * [49] D. Dunlap and V. Kenkre,  Dynamic localization of a charged particle moving under the influence of an electric field,  Phys. Rev. B 34, 3625 (1986).
  * [50] M. Holthaus,  Collapse of minibands in far-infrared irradiated superlattices,  Phys. Rev. Lett. 69, 351 (1992).
  * [51] C. Sias, H. Lignier, Y. P. Singh, A. Zenesini, D. Ciampini, O. Morsch, and E. Arimondo,  Observation of Photon-Assisted Tunneling in Optical Lattices,  Phys. Rev. Lett. 100, 040404 (2008).
  * [52] A. Szameit, I. L. Garanovich, M. Heinrich, A. A. Sukhorukkov, F. Dreisow, T. Pertsch, S. Nolte, A Tünnermann, and Y. S. Kivshar,  Polychromatic dynamic localization in curved photonic lattices,  Nat. Phys. 5, 271 (2009).
  * [53] K. Madison, M. Fischer, R. Diener, Q. Niu, and M. G. Raizen,  Dynamical bloch band suppression in an optical lattice,  Phys. Rev. Lett. 81, 5093 (1998).
  * [54] S. Longhi, M. Marangoni, M. Lobino, R. Ramponi, P. Laporta, E. Cianci, and V. Foglietti,  Observation of dynamic localization in periodically curved waveguide arrays,  Phys. Rev. Lett. 96, 243901 (2006).
  * [55] H. Lignier, C. Sias, D. Ciampini, Y. Singh, A. Zenesini, O. Morsch, and E. Arimondo,  Dynamical Control of Matter-Wave Tunneling in Periodic Potentials,  Phys. Rev. Lett. 99, 220403 (2007).
  * [56] A. Szameit, I. L. Garanovich, M. Heinrich, A. A. Sukhorukov, F. Dreisow, T. Pertsch, S. Nolte, A. Tünnermann, S. Longhi, and Y. S. Kivshar,  Observation of Two-Dimensional Dynamic Localization of Light,  Phys. Rev. Lett. 104, 223903 (2010).
  * [57] E. Haller, R. Hart, M. J. Mark, J. G. Danzl, L. R. llner, and H. Nägerl,  Inducing Transport in a Dissipation-Free Lattice with Super Bloch Oscillations,  Phys. Rev. Lett. 104, 200403 (2010).
  * [58] R. Ma, M. E. Tai, P. M. Preiss, W. S. Bkar, J. Simon, and M. Greiner,  Photon-Assisted Tunneling in a Biased Strongly Correlated Bose Gas,  Phys. Rev. Lett. 107, 095301 (2011).
  * [59] Y. A. Chen, S. Nascimbène, M. Aidelsburger, M. Atala, S. Trotzky, and I. Bloch  Controlling Correlated Tunneling and Superexchange Interactions with ac-Driven Optical Lattices,  Phys. Rev. Lett. 107, 210405 (2011).
  * [60] L. Yuan and S. Fan,  Three-Dimensional Dynamic Localization of Light from a Time-Dependent Effective Gauge Field for Photons,  Phys. Rev. Lett. 114, 243901 (2015).
  * [61] M. Mancini, G. Pagano, G. Cappellini, L. Livi, M. Rider, J. Catani, C. Sias, P. Zoller, M. Inguscio, M. Dalmonte, and L. Fallani,  Observation of chiral edge states with neutral fermions in synthetic Hall ribbons,  Science 349, 1510 (2015).
  * [62] B. K. Stuhl, H. I. Lu, L. M. Aycock, D. Genkina, and I. B. Spielman,  Visualizing edge states with an atomic Bose gas in the quantum Hall regime,  Science 349, 1514 (2015).
  * [63] L. F. Livi, G. Cappellini, M. Diem, L. Franchi, C. Clivati, M. Frittelli, F. Levi, D. Calonico, J. Catani, M. Inguscio, and L. Fallani,  Synthetic Dimensions and Spin-Orbit Coupling with an Optical Clock Transition,  Phys. Rev. Lett. 117, 220401 (2016).
  * [64] M. Atala, M. Aidelsburger, M. Lohse, J. T. Barreiro, B. Paredes, and I. Bloch,  Observation of chiral currents with ultracold atoms in bosonic ladders,  Nat. Phys. 10, 588 (2014).
  * [65] F. A. An, E. J. Meier, and B. Gadway,  Direct observation of chiral currents and magnetic reflection in atomic flux lattices,  Sci. Adv. 3, e1602685 (2017).
  * [66] P. Wang, L. Chen, C. Mi, Z. Meng, L. Huang, K. S. Nawaz, H. Cai, D. W. Wang, S. Y. Zhu, and J. Zhang,  Synthesized magnetic field of a sawtooth superradiance lattice in Bose-Einstein condensates,  npj Quantum Information 6, 18 (2020).
  * [67] See Supplemental Material for the experimental setup, the derivation of the Hamiltonian and the proposal on introducing interactions. It includes Refs. Sambe1973 , Cohen1998 , Ripka2018 , Ding2022 , Henkel2010 , Pupillo2010 , Honer2010 .
  * [68] H. Sambe,  Steady States and Quasienergies of a Quantum-Mechanical System in an Oscillating Field,  Phys. Rev. A 7, 2203 (1973).
  * [69] C. Cohen-Tannoudji, J. Dupont-Roc, and G. Grynberg,  Atom-Photon Interactions (Wiley-VCH Verlag GmbH, Weinheim, Germany, 1998).
  * [70] F. Ripka, H. Kübler, R. Löw, T. Pfau,  A room-temperature single-photon source based on strongly interacting Rydberg atoms,  Science 362, 446 (2018).
  * [71] D. S. Ding, Z. K. Liu, B. S. Shi, G. C, Guo, K. Mølmer, and C. S. Adams,  Enhanced metrology at the critical point of a many-body Rydberg atomic system,  Nat. Phys. (2022). https://doi.org/10.1038/s41567-022-01777-8
  * [72] N. Henkel, R. Nath, and T. Pohl,  Three-Dimensional Roton Excitations and Supersolid Formation in Rydberg-Excited Bose-Einstein Condensates,  Phys. Rev. Lett. 104, 195302 (2010).
  * [73] G. Pupillo, A. Micheli, M. Boninsegni, I. Lesanovsky, and P. Zoller,  Strongly Correlated Gases of Rydberg-Dressed Atoms: Quantum and Classical Dynamics,  Phys. Rev. Lett. 104, 223002 (2010).
  * [74] J. Honer, H. Weimer, T. Pfau, and H. P. Büchler,  Collective Many-Body Interaction in Rydberg Dressed Atoms,  Phys. Rev. Lett. 105, 160404 (2010).
  * [75] D. Yu, B. Peng, X. Chen, and X.-J. Liu, and L. Yuan  Topological holographic quench dynamics in a synthetic frequency dimension,  Light Sci. Appl. 10, 209 (2021).
  * [76] Q. Gao and Q. Niu,  Floquet-Bloch Oscillations and Intraband Zener Tunneling in an Oblique Spacetime Crystal,  Phys. Rev. Lett. 127, 036401 (2021).
  * [77] D. Hügel and B. Paredes,  Chiral ladders and the edges of quantum Hall insulators,  Phys. Rev. A 89, 023619 (2014).
  * [78] Y. He, J. Wu, Y. Hu, J. X. Zhang, and S. Y. Zhu,  Unidirectional reflectionless anti-parity-symmetric photonic lattices of thermal atoms,  Phys. Rev. A 105, 043712 (2022).
  * [79] T. Morimoto, H. C. Po, and A. Vishwanath,  Floquet topological phases protected by time glide symmetry,  Phys. Rev. B 95, 195155 (2017).
  * [80] S. Xu and C. Wu,  Space-Time Crystal and Space-Time Group,  Phys. Rev. Lett. 120, 096401 (2018).
  * [81] Y. Peng and G. Refael,  Floquet Second-Order Topological Insulators from Nonsymmorphic Space-Time Symmetries,  Phys. Rev. Lett. 123, 016806 (2019).
  * [82] Y. Peng,  Topological Space-Time Crystal,  Phys. Rev. Lett. 128, 186802 (2022).
  * [83] Y. Li, H. Cai, D. W. Wang, L. Li, J. Yuan, and W. Li,  Many-Body Chiral Edge Currents and Sliding Phases of Atomic Spin Waves in Momentum-Space Lattice,  Phys. Rev. Lett. 124, 140401 (2020).
  * [84] A. Cao, R. Sajjad, H. Mas, E. Q. Simmons, J. L. Tanlimco, E. Nolasco-Martinez, T. Shimasaki, H. Esat Kondakci, V. Galitski, and D. M. Weld,  Interaction-driven breakdown of dynamical localization in a kicked quantum gas,  Nat. Phys. 18, 1302 (2022).
  * [85] J. H. S. Toh, K. C. McCormick, X. Tang, Y. Su, X. W. Luo, C. Zhang, and S. Gupta,  Many-body dynamical delocalization in a kicked one-dimensional ultracold gas,  Nat. Phys. 18, 1297 (2022).


# Non-Robustness of the Cluster-Robust Inference:
with a Proposal of a New Robust Method††thanks: We benefited from discussions
with A. Colin Cameron and Ulrich Müller. We would like to thank Georgina
Cisneros, Pratap Khattri, Kalika Likhi, Perry Sharaf, Hong Xu and Siyuan Xu
for excellent research assistance. All remaining errors are ours.

Yuya Sasaki and Yulong Wang  Associate professor of economics, Vanderbilt
University. Email<EMAIL_ADDRESS>professor of economics,
Syracuse University. Email<EMAIL_ADDRESS>
###### Abstract

The conventional cluster-robust (CR) standard errors may not be robust. They
are vulnerable to data that contain a small number of large clusters. When a
researcher uses the 51 states in the U.S. as clusters, the largest cluster
(California) consists of about 10% of the total sample. Such a case in fact
violates the assumptions under which the widely used CR methods are guaranteed
to work. We formally show that the conventional CR method fails if the
distribution of cluster sizes follows a power law with exponent less than two.
Besides the example of 51 state clusters, some examples are drawn from a list
of recent original research articles published in a top journal. In light of
these negative results about the existing CR methods, we propose a weighted CR
(WCR) method as a simple fix. Simulation studies support our arguments that
the WCR method is robust while the conventional CR method is not.

Keywords:  cluster-robust inference, non-robustness, power law
JEL Code:  C12, C18

##  1 Introduction

Cluster-robust (CR) standard errors account for within-cluster correlations.
Such correlations often arise by construction within an industry (Hersch,
1998) or within a state (Bertrand et al., 2004), to list a couple of the
earliest examples. Today, even if a model may not induce cluster dependence by
construction, applying CR methods by observable group identifiers is quite
common in practice.

While CR methods have been used for decades, “it is only recently that
theoretical foundations for the use of these methods in many empirically
relevant situations have been developed” (MacKinnon et al., 2022). The initial
theory fixes the (maximum) cluster size $N_{g}$ to a small number and assumes
some large number $G$ of clusters (White, 1984; Liang and Zeger, 1986;
Arellano, 1987). More recent theory (Djogbenou et al., 2019; Hansen and Lee,
2019) allows for large $N_{g}$, but still imposes the restriction
$\sup_{g}N_{g}/N\rightarrow 0$ of vanishing maximum cluster size relative to
the whole sample size $N=\sum_{g=1}^{G}N_{g}$ as $G\rightarrow\infty$.

These conventional assumptions may not be always satisfied. In fact, they do
not hold for certain important situations, such as the case in which a
researcher uses the 51 states in the U.S. as clusters. As documented later in
this article, the distribution of such cluster sizes $N_{g}$ appears to follow
the power law. Consequently, $\sup_{g}N_{g}/N\rightarrow 0$ will not hold. To
fix ideas, think of $\sup_{g}N_{g}/N\approx 0.1\gg 0$ for the state of
California as an outlier in terms of the cluster size in typical survey data.

We show that the assumption of bounded second moments for the CR score fails
if the distribution of $N_{g}$ follows the power law with the tail exponent
less than two. For the example of using the 51 states as clusters,
unfortunately, we cannot rule out the possibility of the tail exponent being
less than two. Therefore, we do not recommend using the conventional CR
methods. We also highlight several examples from a list of original research
articles recently published in a top journal.

We propose a simple fix in light of these adverse findings concerning the
conventional CR methods. Specifically, we propose a weighted CR (WCR) method
to suppress heavy tails in the distributions of cluster sums of scores. We
provide supporting theories to guarantee that this proposed method works.
Simulation studies demonstrate that our proposed WCR method achieves
significantly more accurate coverage than the conventional CR method, thus
supporting our arguments that the WCR method is robust while the conventional
CR method is not.

##  2 Relation to the Literature

The initial theory (White, 1984; Liang and Zeger, 1986; Arellano, 1987) for CR
methods assumes small cluster sizes $N_{g}$ as $G\rightarrow\infty$. It is
implemented by the ‘cluster()’ option or the ‘vce(cluster)’ option by Stata,
and is used by almost all, if not all, empirical papers that report CR
standard errors. It has been known that a large cluster size $N_{g}$ in data
could lead to a large CR standard error (e.g., Cameron and Miller, 2015, p.
324).

More recently, Hansen (2007), Ibragimov and Müller (2010, 2016), Bester,
Conley, and Hansen (2011), and Canay, Santos, and Shaikh (2021) propose
alternative theory to accommodate large $N_{g}$ while $G$ is fixed. Their
frameworks differ from ours in a few aspects. First, they exploit within-
cluster central limit theorem (CLT) which is suitable for panel data, while we
focus on the cross-sectional settings which accommodate arbitrary within-
cluster dependence as opposed to weak dependence. Second, they consider
sequences of increasing cluster sizes $N_{g}$ for all $g$ which is suitable
for long panels, while we treat $N_{g}$ as random variables drawn from a
widely supported distribution to accommodate both small clusters (like the
state of Wyoming) and large clusters (like the state of California) in cross-
sectional settings. Third, they consider small $G$ while we consider large
$G$.

As such, more closely related to this article are the recent developments by
Djogbenou, MacKinnon, and Nielsen (2019) and Hansen and Lee (2019). Both their
frameworks and our framework consider large $G$ with arbitrary within-cluster
dependence as opposed to weak dependence. As mentioned in the introductory
section, however, the existing theory requires $\sup_{g}N_{g}/N\rightarrow 0$,
while $\sup_{g}N_{g}/N\approx 0.1\gg 0$ for the state clusters in typical
survey data in the U.S. Consequently, we provide a fix to the conventional CR
method allowing for cases with $\sup_{g}N_{g}/N\not\rightarrow 0$ unlike the
existing literature. Moreover, treating $N_{g}$ as a random variable, we allow
for substantial heterogeneity in $N_{g}$ across $g$, which is empirically
relevant. Again in typical survey data, we find that
$\sup_{g}N_{g}/\inf_{g}N_{g}\approx 70$ for the state clusters in the U.S.

Besides those discussed above, there are many important papers in the
extensive literature of CR methods. We refer readers to the comprehensive
surveys by Cameron and Miller (2015, 2022) and MacKinnon, Nielsen, and Webb
(2022).

##  3 Review of the Cluster-Robust Inference

Often available for empirical research is a sample of $G$ clusters, where the
$g$-th cluster consists of $N_{g}\in\mathbb{N}$ observations for each
$g\in\\{1,\cdots,G\\}$. A common assumption in this setting is that
observations within a cluster may be arbitrarily correlated, but they are
assumed to be independent across clusters. For instance, the fifty-one states
are often treated as $G=51$ clusters based on a random sample drawn from the
population in the U.S., allowing state-specific factors to induce dependence
among observations within each state possibly.

To fix ideas, suppose that a researcher uses a clustered sample
$\\{\\{(Y_{gi},X_{gi}^{\prime})^{\prime}\\}_{i=1}^{N_{g}}\\}_{g=1}^{G}$ to
study the linear model

$\displaystyle
Y_{gi}=X_{gi}^{\prime}\theta+U_{gi},\qquad\mathbb{E}[U_{g}|X_{g}]=0,$

where $U_{g}=(U_{g1},\cdots,U_{gN_{g}})^{\prime}$ and
$X_{g}=(X_{g1},\cdots,X_{gN_{g}})^{\prime}$ for each $g\in\\{1,\cdots,G\\}$.
Also write $Y_{g}=(Y_{g1},\cdots,Y_{gN_{g}})^{\prime}$ for each
$g\in\\{1,\cdots,G\\}$. Then, the ordinary least squares (OLS) estimator for
the parameter vector $\theta$ under the cluster sampling takes the form of

$\displaystyle\widehat{\theta}=\left(\sum_{g=1}^{G}\sum_{i=1}^{N_{g}}X_{gi}X_{gi}^{\prime}\right)^{-1}\left(\sum_{g=1}^{G}\sum_{i=1}^{N_{g}}X_{gi}Y_{gi}\right)=\left(\sum_{g=1}^{G}X_{g}^{\prime}X_{g}\right)^{-1}\left(\sum_{g=1}^{G}X_{g}^{\prime}Y_{g}\right).$
(1)

Furthermore, commonly employed cluster-robust (CR) variance estimators take
the form of

$\displaystyle\widehat{V}_{\hat{\theta}}^{\text{CR}}=a_{n}\left(\sum_{g=1}^{G}X_{g}^{\prime}X_{g}\right)^{-1}\left(\sum_{g=1}^{G}\widehat{S}_{g}\widehat{S}_{g}^{\prime}\right)\left(\sum_{g=1}^{G}X_{g}^{\prime}X_{g}\right)^{-1},$
(2)

where $a_{n}\rightarrow 1$ is a suitable finite-sample adjustment and
$\widehat{S}_{g}=\sum_{i=1}^{N_{g}}X_{gi}\widehat{U}_{gi}$ with
$\widehat{U}_{gi}=Y_{gi}-X_{gi}^{\prime}\widehat{\theta}.$ In particular, a
majority of empirical economics papers is based on the ‘cluster()’ option or
the ‘vce(cluster)’ option in Stata, which uses (1) and (2) with the finite-
sample adjustment

$\displaystyle
a_{n}=\left(\frac{\sum_{g=1}^{G}N_{g}-1}{\sum_{g=1}^{G}N_{g}-p}\right)\left(\frac{G}{G-1}\right),$

where $p$ denotes the dimension of $X_{gi}$. Note that this adjustment factor
satisfies $a_{n}\rightarrow 1$ as $G\rightarrow\infty$.

##  4 A Pitfall in the Cluster-Robust Inference

While it provides robustness to cluster dependence, the CR method is
vulnerable to common cross-sectional situations in which there are a small
number of large clusters.

The standard econometric theory to guarantee that the OLS estimator (1) and
the CR variance estimator (2) behave well under the cluster sampling are based
on the asymptotic property

$\displaystyle\sqrt{G}\left(\widehat{\theta}-\theta\right)=\Big{(}\underbrace{\frac{1}{G}\sum_{g=1}^{G}\Xi_{g}}_{\stackrel{{\scriptstyle
p}}{{\longrightarrow}}\
Q}\Big{)}^{-1}\Big{(}\underbrace{\frac{1}{\sqrt{G}}\sum_{g=1}^{G}S_{g}}_{\stackrel{{\scriptstyle
d}}{{\longrightarrow}}\mathcal{N}(0,V)}\Big{)}\stackrel{{\scriptstyle
d}}{{\longrightarrow}}\mathcal{N}(0,Q^{-1}VQ^{-1})$ (3)

as $G\rightarrow\infty$, where $\stackrel{{\scriptstyle p}}{{\rightarrow}}$
stands for convergence in probability, $\stackrel{{\scriptstyle
d}}{{\rightarrow}}$ stands for convergence in distribution, and

$\displaystyle Q=\mathbb{E}[\Xi_{g}],$
$\displaystyle\Xi_{g}=\sum_{i=1}^{N_{g}}X_{gi}X_{gi}^{\prime},$ $\displaystyle
V=\text{Var}(S_{g}),$ $\displaystyle S_{g}=\sum_{i=1}^{N_{g}}X_{gi}U_{gi}.$

Indeed, the cluster-robust approach allows for robustness against within-
cluster dependence. However, this robustness is not cost-free. Namely,
suitable moments of $\Xi_{g}$ and $S_{g}$ need to be finite for the weak law
of large numbers and the central limit theorem to be invoked in (3).
Specifically, the asymptotic normality (3) and the validity of the CR variance
estimator (2) under the cluster sampling require $S_{g}$ to have finite second
moments. If the second moment of $S_{g}$ does not exist, then
$\sqrt{G}\left(\widehat{\theta}-\theta\right)$ diverges and hence would not
converge in distribution as in (3). However, we are going to argue that this
bounded second moment condition for the cluster score $S_{g}$ can fail under
the cluster sampling even if the second moment of the individual’s score
$X_{gi}U_{gi}$ were finite.

For ease of illustration, we first introduce a few definitions and notations.
For a given $k\in\\{1,\cdots,p\\}$, let $\Sigma_{g}$ and $Z_{gi}$ be short-
hand notations for the $k$-the coordinate of $S_{g}$ and the $k$-th coordinate
of $X_{gi}U_{gi}$, respectively – we omit the dependence on $k$ in these
notations for brevity. For any distribution function $F$, we say that $F$ is
regularly varying (RV) at infinity if it satisfies the following property that

$\frac{1-F(xt)}{1-F(t)}\rightarrow x^{-\alpha}\text{ as }t\rightarrow\infty,$

for any $x>0$ and some constant $\alpha>0$. The constant $\alpha$ is referred
to as the tail exponent, which measures the tail heaviness of $F$. Let $F$
denote the marginal distribution of $Z_{gi}$. For each $n\in\\{2,3,\cdots\\}$,
let $C^{n}$ denote the copula that characterizes the joint distribution of
$(Z_{g1},\dots,Z_{gn})$ such that

$\mathbb{P}\left(Z_{g1}\leq z_{1},\dots,Z_{gn}\leq
z_{n}\right)=C^{n}(F(z_{1}),\dots,F(z_{n})).$

With these definitions and notations, we make the following set of
assumptions.

######  Assumption 1.

The following conditions are satisfied.

  1. 1.

$\\{Z_{gi}\\}_{gi}$ is identically distributed with CDF $F$, which is RV at
infinity with coefficient $\alpha>1$.

  2. 2.

The copula density
$c(u_{1},\dots,u_{n})={\partial^{n}C(u_{1},\dots,u_{n})}/{\partial
u_{1}\cdots\partial u_{n}}$ exists and is uniformly bounded.

  3. 3.

$N_{g}$ is independent of $(Z_{g1},Z_{g2},....)$ and its distribution $H$ is
RV at infinity with coefficient $\beta>1$.

We provide some discussions about Assumption 1. Assumption 1.1 requires that
the CDF of $Z_{gi}$ implies a finite mean and its tail is regularly varying.
The regularly varying tail is mild and satisfied by many commonly used heavy-
tailed distributions such as Pareto, Student-t, Cauchy, F distributions.
Without loss of generality, we focus on the right tail and assume that
$Z_{gi}$ is non-negative. Otherwise, our conclusion still holds provided

$\lim_{z\rightarrow\infty}\frac{\mathbb{P}(Z_{gi}<-z)}{Z_{gi}>z}=0.$

Given the regularly varying tail, we can elegantly characterize the moment
conditions (e.g., Mikosch, 1999, Proposition 1.3.2) as

$\displaystyle\mathbb{E}[|Z_{gi}|^{r}]$ $\displaystyle=\infty\text{ for any
}r<\alpha$ $\displaystyle\mathbb{E}[|Z_{gi}|^{r}]$ $\displaystyle<\infty\text{
for any }r>\alpha.$

For conciseness, we suppress the distributions with bounded infinite-order
moments (such as Gaussian distribution) and those with a compact support. They
are commonly referred to as thin-tailed distributions, which are still covered
by our conclusion. See the simulation studies in Section 8. Assumption 1.2
allows for dependence among $Z_{g1},\dots,Z_{gN_{g}}$ within each cluster $g$.
This condition allows for very general formats of correlations within each
cluster. Assumption 1.3 requires that the distribution of the cluster size
$N_{g}$ is also regularly varying. This condition is plausibly satisfied as
shown in the following section. Since $N_{g}$ takes integer values, we may
consider it as the integer part of some continuous random variable, say
$N^{*}_{g}$ whose tail is regularly varying. Given that $N^{*}_{g}-1\leq
N_{g}\leq N^{*}_{g}$, we have that
$\mathbb{P}(N_{g}>y)\sim\mathbb{P}(N^{*}_{g}>y)$ as $y\rightarrow\infty$ and
hence continue with $N_{g}$ for conciseness.

Suppose that $\beta<\alpha$ under Assumption 1.1 and 1.3. In other words, the
tail of the distribution of $N_{g}$ dominates that of $Z_{gi}$. In this case,
the tail heaviness of the summation $\Sigma_{g}=\sum_{i=1}^{N_{g}}Z_{gi}$ is
dictated by that of $N_{g}$. Therefore, even if $Z_{gi}$ has a finite $r$-th
moment for $r<\alpha$, the $r$-th moment of $\Sigma_{g}$ might still be
infinite if $r>\beta$. The following theorem formalizes this argument.

######  Theorem 1.

Suppose Assumption 1 holds. If $\beta<\alpha$, then for any $z>0$,

$\frac{\mathbb{P}(\Sigma_{g}>zt)}{\mathbb{P}(\Sigma_{g}>t)}=z^{-\beta}\text{
as }t\rightarrow\infty.$

We provide a proof in Appendix A.1.

Under i.i.d. sampling, the ‘robust’ standard errors (e.g., those based on
Eicker-Huber-White variance estimator or the robust option in Stata) exist and
behave well if $Z_{gi}$ has a bounded second moment. Under cluster sampling,
however, the cluster-‘robust’ standard errors (i.e., those based on (2) or the
cluster()/vce(cluster) options in Stata) may fail to exist even if $Z_{gi}$
has a bounded second moment. Theorem 1 implies that $\Sigma_{g}$ may have an
infinite second moment even if $Z_{gi}$ has a bounded second moment. This
occurs when $\beta<2<\alpha$.

##  5 Example: Problem with the 51 States as Custers

A natural question now is whether the distribution of cluster sizes $N_{g}$
has a heavy tail with tail index $\beta<2$. An answer will of course depend on
data and the context of empirical research. Let us focus on one of the most
common situations of empirical economic analyses with clustering.
Specifically, suppose that a researcher clusters a sample of individuals in
the U.S. Panel Study of Income Dynamics (PSID) by the 51 states. We extract
the sample of all the male individuals from the most recent wave of 2019. This
sample consists of $4808$ observations across $G=51$ states.

Figure 1: The log-log plot for the subsamples of male individuals by state
clusters in the U.S. PSID wave of 2019.

In this data set, the five largest clusters have their sizes of $N_{4}=629$
(California), $N_{42}=430$ (Texas), $N_{32}=366$ (North Carolina),
$N_{21}=310$ (Michigan), and $N_{39}=304$ (South Carolina). Figure 1 plots the
logarithm of the rank of $N_{g}$ (in descending order) against the logarithm
of $N_{g}$ for all the 51 clusters. The gray line indicates the fitted line
using the largest 25 observations, i.e., 50% of the 51 states. Note that these
largest 25 states follow the line of slope $-1.76$, while the remaining 26
states follow another line of a less steep slope. This shape of the log-log
plot implies that the distribution of the cluster sizes $N_{g}$ is
approximately asymmetric double Pareto. The slope of $-1.76$ can be
interpreted as the negative value of the Pareto exponent $\beta$ in the right
tail of the distribution $H$ of $N_{g}$.

Figure 2: The Hill plot for the largest 25 (50%) of the 51 state clusters in
the U.S. PSID wave of 2019: men.

From this graphical analysis, it seems quite possible that $\beta<2$. We now
resort to a more formal and widely practiced test to reinforce this heuristic
observation. Figure 2 displays the so-called Hill plot (e.g., Drees et al.,
2000), consisting of 95% confidence intervals for the Parteto exponent $\beta$
of the distribution of $N_{g}$ using the largest $k$ order statistics for
$k\in\\{2,\cdots,25\\}$. Observe that we cannot rule out the possibility of
$\beta<2$ for any $k$.

Recall from Theorem 1 that the Pareto exponent $\beta$ of the distribution $H$
of $N_{g}$ dictates the heaviness of the distribution of the coordinates
$\Sigma_{g}$ of the score $S_{g}$. If $\beta<2$, then the asymptotic
convergence (3) fails even if $\alpha>2$ holds for each coordinate $Z_{gi}$ of
the score $X_{gi}U_{gi}$. In light of the conclusion from the previous two
paragraphs, we cannot rule out $\beta<2$. Therefore, it follows that the
widely employed CR variance estimator (3) may lead to misleading inference. In
summary, we do not recommend the common practice of using the CR standard
error with the 51 states as clusters for this data set from the U.S. PSID.

##  6 Examples from Econometrica, 2020

We highlight several examples based on recent original research articles
published in 2020 in Econometrica that use cluster-robust methods of inference
for which we cannot rule out the possibility of $\beta<2$. Namely we showcase
in point with state clusters in the U.S. (Burstein et al., 2020), region
clusters in Russia (Enikolopov et al., 2020), and NGO branch clusters (Alfonsi
et al., 2020).

We would like to emphasize that the main point of this section is not to
criticize a specific list of papers – the authors of these papers are in fact
great for allowing us to replicate their analyses more easily than many
others. Our selection from one journal from only one recent year reflects only
the tip of the iceberg – furthermore, there are even more articles published
in Econometrica and other journals for which data are unavailable to us or
replication was difficult for us. We suspect that an enormous number of
empirical research papers from a wide variety of journals and from a long
history of the literature are in fact subject to this problem of non-
robustness.

###  6.1 State Clusters in the U.S. (Burstein et al., 2020)

Burstein et al. (2020) use 56 state clusters in the U.S.111The cluster
variable is the state as in Section 5. We selected this example, however,
because the data are different and the numbers of state categories are
different too. for their analysis summarized in their Tables 1–2. We extracted
their cluster variable and draw its Hill plot for the largest 28 (50%) in
Figure 3.

Figure 3: The Hill plot for the largest 28 (50%) of the 56 state clusters for
(Burstein et al., 2020).

For every number $k$ of the top order statistics except for $k=20$ and 21, the
95% confidence interval does not exclude the possibility of $\beta<2$.
Therefore, we cannot plausibly assume $\beta>2$ for this data set to use the
conventional CR method. Furthermore, using the formal statistical test
presented in Appendix B, we reject the hypothesis of
$\mathbb{E}[S_{g}^{2}]<\infty$ for a number of regressions whose results are
reported in their Tables 1–2.

###  6.2 Region Clusters in Russia (Enikolopov et al., 2020)

Enikolopov et al. (2020) use 78 region clusters in Russia for their analysis
summarized in their Tables 1–3. We extract their cluster variable and draw its
Hill plot for the largest 39 (50%) in Figure 4.

Figure 4: The Hill plot for the largest 39 (50%) of the 78 region clusters for
(Enikolopov et al., 2020).

For every number $k$ of the top order statistics, the 95% confidence interval
does not exclude the possibility of $\beta<2$. Therefore, we cannot plausibly
assume $\beta>2$ for this data set to use the conventional CR method.
Furthermore, using the formal statistical test presented in Appendix B, we
reject the hypothesis of $\mathbb{E}[S_{g}^{2}]<\infty$ for a number of
regressions whose results are reported in their Tables 1–3.

###  6.3 Branch Clusters (Alfonsi et al., 2020)

Alfonsi et al. (2020) use 108 NGO (BRAC) branch clusters for for their
analysis summarized in their Table 7. We extract their cluster variable and
draw its Hill plot for the largest 54 (50%) in Figure 5.

Figure 5: The Hill plot for the largest 54 (50%) of the 108 BRAC branch
clusters for (Alfonsi et al., 2020).

For every number $k$ of the top order statistics, the 95% confidence interval
does not exclude the possibility of $\beta<2$. Therefore, we cannot plausibly
assume $\beta>2$ for this data set to use the conventional CR method.
Furthermore, using the formal statistical test presented in Appendix B, we
reject the hypothesis of $\mathbb{E}[S_{g}^{2}]<\infty$ for a number of
regressions whose results are reported in their Table 7.

###  6.4 How Many Research Papers Cluster Standard Errors?

Once again, we want to stress that the above list of examples is just the tip
of the iceberg. Our selection of the three papers is due only to the relative
ease of replication thanks to the efforts by the authors of these papers.
There were many other papers for which replication was difficult or
infeasible, sometimes due to limited data availability. It is worthy to note
that there are at least nine papers that cluster standard errors among those
published in American Economic Review in year 2020, and at least eight papers
that cluster standard errors among those published in Econometrica in year
2020, among others.222The results of the test of the null hypothesis
$\mathbb{E}[S_{g}^{2}]<\infty$ for all these papers, as well as those in other
top five journals, based on the method presented in Appendix B, will become
available in spreadsheet upon request from the authors under certain
conditions. We thank our graduate research assistants for their great efforts
to replicate the analyses by the large number of papers.

##  7 Our Proposal of a New Robust Approach

This section proposes a new approach to cluster-robust inference. It does not
suffer from the aforementioned problem with the standard CR method
characterized by Theorem 1. Our proposal is simple, but we provide a complete
theoretical rationale for why it works unlike the standard CR approach.

One simple way to fix the problem of the non-robustness of the standard CR
method is to modify (1) and (2) by

$\displaystyle\widehat{\theta}^{\text{WCR}}=\left(\sum_{g=1}^{G}N_{g}^{-1}X_{g}^{\prime}X_{g}\right)^{-1}\left(\sum_{g=1}^{G}N_{g}^{-1}X_{g}^{\prime}Y_{g}\right)\qquad\text{and}$
$\displaystyle\widehat{V}_{\hat{\theta}}^{\text{WCR}}=a_{n}\left(\sum_{g=1}^{G}N_{g}^{-1}X_{g}^{\prime}X_{g}\right)^{-1}\left(\sum_{g=1}^{G}N_{g}^{-2}\widehat{S}_{g}\widehat{S}_{g}^{\prime}\right)\left(\sum_{g=1}^{G}N_{g}^{-1}X_{g}^{\prime}X_{g}\right)^{-1},$

respectively, where $a_{n}\rightarrow 1$.

As in Section 4, for given $k\in\\{1,\cdots,p\\}$, let $\Sigma_{g}$ and
$Z_{gi}$ be short-hand notations for the $k$-the coordinate of $S_{g}$ and the
$k$-th coordinate of $X_{gi}U_{gi}$, respectively. Define
$\widetilde{\Sigma}_{g}:=\Sigma_{g}/N_{g}=\sum_{i=1}^{N_{g}}Z_{gi}$. In
contrast to Theorem 1, the following theorem establishes that
$\widetilde{\Sigma}_{g}$ has the same tail index as the original score
$Z_{gi}$.

######  Theorem 2.

If Assumption 1 holds, then for any $z>0$,

$\frac{\mathbb{P}(\widetilde{\Sigma}_{g}>zt)}{\mathbb{P}(\widetilde{\Sigma}_{g}>t)}=z^{-\alpha}\text{
as }t\rightarrow\infty.$

We present a proof in Appendix A.2.

Recall that the non-robust methods (e.g., those based on Eicker-Huber-White
standard errors) requires only bounded second moments of the score
$X_{gi}U_{gi}$. Also recall from Section 4 that the conventional CR method can
fail when the second moment of $N_{g}$ is infinite even if the second moments
of the score $X_{gi}U_{gi}$ were finite. Now, Theorem 2 states that our
proposed robust approach works as far as the second moments of the score
$X_{gi}U_{gi}$ is bounded regardless of whether the second moment of $N_{g}$
if finite or not.

######  Assumption 2.

1\. $(N_{g},X_{g},Y_{g})$ is i.i.d. across $g$. 2\.
$\mathbb{E}[N_{g}^{-1}X_{g}^{\prime}X_{g}]$ is non-singular.

Part 1 of this assumption requires the i.i.d. sampling across clusters, as is
assumed under cluster sampling environments. It does, however, allow for
arbitrary dependence within each cluster, again, as is assumed under cluster
sampling environments. Part 2 of the above assumption requires rules out
multi-collinearity, as is usually assumed.

Recall that our main claim from Section 4 is that the asymptotic normality (3)
for the conventional CR inference may not work if $\beta<2$ even if
$\alpha>2$. Section 5 demonstrates that we cannot rule out such a pathetic
case with $\beta<2$ even for the most common setting of cluster sampling,
namely using 51 states as clusters. On the other hand, the following theorem
shows that our proposed CR approach based on $\widehat{\theta}^{\text{WCR}}$
and$\widehat{V}_{\hat{\theta}}^{\text{WCR}}$ works as far as $\alpha>2$ is
true, regardless of whether $\beta<2$ is true or not.

######  Theorem 3.

If Assumptions 1 and 2 hold with $\alpha>2$, then

$\sqrt{G}(\widehat{\theta}^{\text{WCR}}-\theta)\xrightarrow{p}\mathcal{N}(0,V^{\text{WCR}})$

holds as $G\rightarrow\infty$, where

$V^{\text{WCR}}=(\mathbb{E}[N_{g}^{-1}X_{g}^{\prime}X_{g}])^{-1}(\mathbb{E}[N_{g}^{-2}X_{g}^{\prime}X_{g}U^{2}_{g}])(\mathbb{E}[N_{g}^{-1}X_{g}^{\prime}X_{g}])^{-1}.$

Furthermore, we have
$G\widehat{V}_{\hat{\theta}}^{\text{WCR}}\xrightarrow{p}V^{\text{WCR}}$ as
$G\rightarrow\infty$

A proof is presented in Appendix A.3. We remark that the RV condition in
Assumption 1.1 is assumed only for concise statement and is not essential. We
can more trivially accommodate thin-tailed distributions, such as the Gaussian
distribution, for the individual score.

As a consequence of this theorem, one can use

$\widehat{\theta}^{\text{WCR}}\pm 1.96\widehat{V}_{\hat{\theta}}^{\text{WCR}}$

as a robust 95% confidence interval even if $\beta<2$ as is likely the case
with 51 states used as clusters.

##  8 Simulations

In this section, we present simulation studies to evaluate the finite sample
performance of our proposed method (WCR) in comparison with the conventional
cluster-robust (CR) method.

The data generating design is defined as follows. We draw cluster sizes

$N_{g}\sim\lceil 10\cdot\text{Pareto}(1,\beta)\rceil$

independently for $g\in\\{1,\cdots,G\\}$, where $\lceil a\rceil$ denotes the
smallest integer greater than or equal to $a$. For each
$g\in\\{1,\cdots,G\\}$, we draw two $N_{g}$-variate vectors

$\displaystyle(X_{g11},\cdots,X_{gN_{g}1})^{\prime}$ $\displaystyle\sim
N(0,\Omega)\qquad\text{and}$
$\displaystyle(U_{g1},\cdots,U_{gN_{g}})^{\prime}$ $\displaystyle\sim
N(0,\Omega),$

where $\Omega$ is an $N_{g}\times N_{g}$ variance-covariance matrix such that
$\Omega_{ii}=1$ for all $i\in\\{1,\cdots,N_{g}\\}$ and $\Omega_{ij}=\rho$
whenever $i\neq j$. Note that the correlation parameter $\rho$ models the
dependence within each cluster. We then construct

$Y_{gi}=\theta_{0}+\theta_{1}X_{gi1}+U_{gi}$

for each individual $i\in\\{1,\cdots,N_{g}\\}$ for each cluster
$g\in\\{1,\cdots,G\\}$.

We vary values of the exponent parameter $\beta\in\\{1.25,1.50,1.75,2.00\\}$
and the correlation parameter $\rho\in\\{0.25,0.50,0.75\\}$ across sets of
simulations. The regression coefficients are fixed to
$(\theta_{0},\theta_{1})^{\prime}=(1,1)^{\prime}$ throughout. The sample sizes
(i.e., the numbers of clusters) are varied as $G\in\\{25,50\\}$ across sets of
simulations. Each set of simulations consists of 1,000 Monte Carlo iterations.

Coverage Frequencies for the Slope Parameter $\theta_{1}$ by the 95%
Confidence Intervals
---
|  | Method | $\beta=2.00$ | $\beta=1.75$ | $\beta=1.50$ | $\beta=1.25$
$G=25$ | $\rho=0.25$ | Conventional CR | 0.90 | 0.89 | 0.89 | 0.85
|  | New WCR | 0.93 | 0.93 | 0.93 | 0.94
| $\rho=0.50$ | Conventional CR | 0.88 | 0.88 | 0.88 | 0.83
|  | New WCR | 0.91 | 0.92 | 0.93 | 0.94
| $\rho=0.75$ | Conventional CR | 0.86 | 0.86 | 0.86 | 0.81
|  | New WCR | 0.89 | 0.92 | 0.92 | 0.93
$G=50$ | $\rho=0.25$ | Conventional CR | 0.92 | 0.91 | 0.89 | 0.84
|  | New WCR | 0.94 | 0.94 | 0.93 | 0.95
| $\rho=0.50$ | Conventional CR | 0.90 | 0.89 | 0.87 | 0.83
|  | New WCR | 0.93 | 0.93 | 0.93 | 0.95
| $\rho=0.75$ | Conventional CR | 0.90 | 0.89 | 0.86 | 0.82
|  | New WCR | 0.93 | 0.93 | 0.92 | 0.94

Table 1: Coverage frequencies of the slope parameter $\theta_{1}$ by the 95%
confidence intervals based on the conventional cluster-robust (CR) standard
error and our proposed method (WCR). The results are based on 1,000 Monte
Carlo iterations.

Table 1 reports coverage frequencies of the slope parameter $\theta_{1}$ by
the 95% confidence intervals based on the conventional cluster-robust (CR)
standard error and our proposed method (WCR). Observe that the conventional CR
method yields significant under-coverage or over-rejection. Its results tend
to exacerbate as the within-cluster correlation $\rho$ increases and as the
tail heaviness $1/\beta$ of the distribution of $N_{g}$ increases. On the
other hand, our proposed WCR method achieves more accurate coverage for the
nominal probability of 95%, and its results are less sensitive to variations
in $\rho$ or $1/\beta$. Overall, these results demonstrate the outperformance
of the WCR method over the conventional CR method across various data
generating designs.

##  9 Summary and Discussions

CR standard errors are used extensively in the practice of empirical economic
analyses. Despite its wide practice, the literature has been silent about a
potential pitfall of the conventional CR methods. The theory for the
conventional CR methods assumes that (i) $N_{g}$ is fixed while
$G\rightarrow\infty$, or at best (ii) $\sup_{g}N_{g}/N\rightarrow 0$ as
$G\rightarrow\infty$. Even the most common empirical settings, such as the
case of using the 51 states in the U.S. as clusters, violate these existing
assumptions.

We establish three theoretical results in this paper. First, if the
distribution of $N_{g}$ follows the power low with exponent less than two,
then the second moment of the CR score fails to exist even if the second
moment of the individual score were finite. Second, the second moment of our
proposed WCR score exists as far as the second moment of the individual score
exists regardless of the exponent of the distribution of $N_{g}$. Third,
consequently, our proposed WCR method enjoys robust inference unlike the
conventional CR method. It is also worthy to note that, in the special case
where the tail exponent $\beta$ is large, the WCR method behaves similarly to
the conventional CR method.

Indeed, we could rely on self-normalized CLT to establish an alternative
asymptotic theory for the conventional CR method even if the distribution of
$N_{g}$ had its tail exponent less than two. In that case, however, the rate
of convergence will be slower than $\sqrt{G}$ whereas our proposed WCR
estimator achieves the $\sqrt{G}$ rate of convergence. As evidenced by our
simulation results, inference based on the conventional CR method is
inaccurate unlike that based on our proposed WCR method. Therefore, we do not
recommend using the conventional CR method anyway.

## References

  * Alfonsi et al. (2020) Alfonsi, L., O. Bandiera, V. Bassi, R. Burgess, I. Rasul, M. Sulaiman, and A. Vitali (2020): “Tackling youth unemployment: Evidence from a labor market experiment in Uganda,” _Econometrica_ , 88, 2369–2414.
  * Arellano (1987) Arellano, M. (1987): “Computing robust standard errors for within-groups estimators,” _Oxford Bulletin of Economics and Statistics_ , 49, 431–434.
  * Bertrand et al. (2004) Bertrand, M., E. Duflo, and S. Mullainathan (2004): “How much should we trust differences-in-differences estimates?” _Quarterly Journal of Economics_ , 119, 249–275.
  * Bester et al. (2011) Bester, C. A., T. G. Conley, and C. B. Hansen (2011): “Inference with dependent data using cluster covariance estimators,” _Journal of Econometrics_ , 165, 137–151.
  * Burstein et al. (2020) Burstein, A., G. Hanson, L. Tian, and J. Vogel (2020): “Tradability and the Labor-Market Impact of Immigration: Theory and Evidence From the United States,” _Econometrica_ , 88, 1071–1112.
  * Cameron and Miller (2015) Cameron, A. C. and D. L. Miller (2015): “A practitioner’s guide to cluster-robust inference,” _Journal of Human Resources_ , 50, 317–372.
  * Cameron and Miller (2022) ——— (2022): “Recent developments in cluster-robust inference,” _Working Paper_.
  * Canay et al. (2021) Canay, I. A., A. Santos, and A. M. Shaikh (2021): “The wild bootstrap with a “small” number of “large” clusters,” _Review of Economics and Statistics_ , 103, 346–363.
  * Djogbenou et al. (2019) Djogbenou, A. A., J. G. MacKinnon, and M. Ø. Nielsen (2019): “Asymptotic theory and wild bootstrap inference with clustered errors,” _Journal of Econometrics_ , 212, 393–412.
  * Drees et al. (2000) Drees, H., S. Resnick, and L. De Haan (2000): “How to make a Hill plot,” _Annuals of Statistics_ , 28, 254–274.
  * Enikolopov et al. (2020) Enikolopov, R., A. Makarin, and M. Petrova (2020): “Social media and protest participation: Evidence from Russia,” _Econometrica_ , 88, 1479–1514.
  * Hansen and Lee (2019) Hansen, B. E. and S. Lee (2019): “Asymptotic theory for clustered samples,” _Journal of econometrics_ , 210, 268–290.
  * Hansen (2007) Hansen, C. B. (2007): “Asymptotic properties of a robust variance matrix estimator for panel data when T is large,” _Journal of Econometrics_ , 141, 597–620.
  * Hersch (1998) Hersch, J. (1998): “Compensating differentials for gender-specific job injury risks,” _American Economic Review_ , 88, 598–607.
  * Ibragimov and Müller (2010) Ibragimov, R. and U. K. Müller (2010): “t-Statistic based correlation and heterogeneity robust inference,” _Journal of Business & Economic Statistics_, 28, 453–468.
  * Ibragimov and Müller (2016) ——— (2016): “Inference with few heterogeneous clusters,” _Review of Economics and Statistics_ , 98, 83–96.
  * Liang and Zeger (1986) Liang, K.-Y. and S. L. Zeger (1986): “Longitudinal data analysis using generalized linear models,” _Biometrika_ , 73, 13–22.
  * MacKinnon et al. (2022) MacKinnon, J. G., M. Ø. Nielsen, and M. D. Webb (2022): “Cluster-robust inference: A guide to empirical practice,” _Journal of Econometrics_.
  * Mikosch (1999) Mikosch, T. (1999): _Regular variation, subexponentiality and their applications in probability theory_ , vol. 99, Eindhoven University of Technology Eindhoven, The Netherlands.
  * Sasaki and Wang (2022) Sasaki, Y. and Y. Wang (2022): “Diagnostic testing of finite moment conditions for the consistency and root-n asymptotic normality of the gmm and m estimators,” _Journal of Business & Economic Statistics (forthcoming)_.
  * White (1984) White, H. (1984): _Asymptotic Theory for Econometricians_ , Academic Press.

## Appendix

##  Appendix A Proofs

###  A.1 Proof of Theorem 1

###### Proof.

We say that sequence $\\{Z_{gi}\\}_{i=1}^{n}$ of random variables is extended
negatively dependent (END) if there exists a finite constant $M>0$ such that
both

$\mathbb{P}\left(Z_{g1}\leq z_{1},...,Z_{gn}\leq z_{n}\right)\leq
M\prod_{i=1}^{n}\mathbb{P}\left(Z_{gi}\leq z_{i}\right)$

and

$\mathbb{P}\left(Z_{g1}>z_{1},...,Z_{gn}>z_{n}\right)\leq
M\prod_{i=1}^{n}\mathbb{P}\left(Z_{gi}>z_{i}\right)$

hold for all $z_{1},...,z_{n}$.

Assumptions 1.1 implies $\mathbb{E}\left[N_{g}\right]<\infty$. Assumption 1.2
implies that $Z_{g1},\dots,Z_{gn}$ are END for every $n\geq 1$ (Liu, 2009,
Example 4.2). Finally, Assumption 1.3 and $\beta<\alpha$ imply
$\overline{F}\left(x\right)=1-F(x)=o\left(\overline{H}\left(x\right)\right)$.
Then, Theorem 3.2 of Feng (2015) yields

$\mathbb{P}(\Sigma_{g}>x)\sim\overline{H}\left(x/\mathbb{E}\left[N_{g}\right]\right)\sim\overline{H}\left(x\right).$
(4)

Note also that Assumption 1.3 implies

$\frac{\overline{H}\left(xt\right)}{\overline{H}\left(t\right)}=x^{-\beta}\text{
as }t\rightarrow\infty.$ (5)

A proof of the statement is now complete by combining (4) and (5). ∎

###  A.2 Proof of Theorem 2

###### Proof.

Throughout the proof, we let $C$ denote a generic finite constant, whose value
could change line-by-line. Define

$\theta_{i}=\frac{1}{N_{g}}\mathbf{1}[i\leq N_{g}]$

for each $i\in\mathbb{N}$, where $\mathbf{1}[\cdot]$ denotes the indicator
function. Then, we can write

$\widetilde{\Sigma}_{g}=\sum_{i=1}^{N}\frac{Z_{gi}}{N_{g}}=\sum_{i=1}^{\infty}\theta_{i}Z_{gi}.$

We first show that

$\sum_{i=1}^{\infty}\left(\mathbb{E}[\theta_{i}^{\alpha-\epsilon}]\right)^{\frac{1}{\alpha+\epsilon}}<\infty$
(6)

holds for some $\epsilon\in[0,\alpha]$. To this end, we use Assumption 1.3 to
obtain

$\displaystyle\mathbb{E}\left[\theta_{i}^{\alpha-\epsilon}\right]$
$\displaystyle=\mathbb{E}\left[\frac{1}{N_{g}^{\alpha-\epsilon}}\mathbf{1}[i\leq
N_{g}]\right]$
$\displaystyle\leq\mathbb{E}\left[\frac{1}{i^{\alpha-\epsilon}}\mathbf{1}[i\leq
N_{g}]\right]$ $\displaystyle=\frac{1}{i^{\alpha-\epsilon}}\mathbb{P}[i\leq
N_{g}]$ $\displaystyle\leq Ci^{-\beta-\alpha+\epsilon}.$

Thus, the inequality (6) follows by $\beta>\alpha>1$.

Second, we show that

$\sum_{i=1}^{\infty}\mathbb{P}\left(\theta_{i}Z_{gi}>z\right)<\infty$ (7)

holds for any $z$. Using Markov’s inequality and the fact that
$\mathbb{E}[Z_{gi}]<\infty$ holds by Assumption 1.1, we obtain

$\displaystyle\sum_{i=1}^{\infty}\mathbb{P}\left(\theta_{i}Z_{gi}>z\right)\leq$
$\displaystyle\sum_{i=1}^{\infty}\frac{\mathbb{E}\left[N_{g}^{-1}\mathbf{1}[i\leq
N_{g}]Z_{gi}\right]}{z}$ $\displaystyle\leq$
$\displaystyle\sum_{i=1}^{\infty}\frac{C\mathbb{E}[Z_{gi}]i^{-1-\beta}}{z}$
$\displaystyle\leq$ $\displaystyle C<\infty.$

Finally, the extended negatively dependence condition is sufficient for the
tail independence as required by equation (2.9) in Zhang, Shen, and Wang
(2009). By our equations (6) and (7), their Conditions H1 and H2 are
satisfied. Then, their Theorem 3.1 and Remark 3.2 yield

$\mathbb{P}\left(\widetilde{\Sigma}_{g}>z\right)\sim(1-F(z))\sum_{i=1}^{\infty}\mathbb{E}[\theta_{i}^{\alpha}]\text{
as }z\rightarrow\infty,$

where $\sum_{i=1}^{\infty}\mathbb{E}[\theta_{i}^{\alpha}]$ is finite by
similar lines of argument to those used to derive (6) above. A proof is now
complete by Assumption 1.3. ∎

###  A.3 Proof of Theorem 3

###### Proof.

We first establish the asymptotic normality of
$\widehat{\theta}^{\text{WCR}}$. Given Assumption 2, it suffices to establish

$\mathbb{E}[N_{g}^{-1}X^{\prime}_{g}X_{g}]<\infty$ (8)
$\mathbb{E}[N_{g}^{-2}X^{\prime}_{g}X_{g}U^{2}_{g}]<\infty.$ (9)

Without loss of generality, we show the case with a scalar $X$. For (8), the
identical distribution of $Z_{gi}$ across $i$ (Assumption 1.1) and the
independence between $N_{g}$ and $Z_{gi}$ (Assumption 1.3) imply

$\mathbb{E}[N_{g}^{-1}X^{\prime}_{g}X_{g}]=\mathbb{E}[N_{g}^{-1}\sum^{N_{g}}_{i=1}X^{2}_{gi}]=\mathbb{E}[X^{2}_{gi}]<\infty.$

For (9), we have

$\mathbb{E}[N_{g}^{-2}X^{\prime}_{g}X_{g}U^{2}_{g}]=\mathbb{E}[N_{g}^{-2}\sum^{N_{g}}_{i=1}X^{2}_{gi}U^{2}_{gi}]\overset{(1)}{=}\mathbb{E}[N_{g}^{-1}]\mathbb{E}[X^{2}_{gi}U^{2}_{gi}]\overset{(2)}{<}\mathbb{E}[X^{2}_{gi}U^{2}_{gi}]\overset{(3)}{<}\infty,$

where (1) follows from the identical distribution of $Z_{gi}$ across $i$
(Assumption 1.1) and the independence between $N_{g}$ and $Z_{gi}$ (Assumption
1.3), (2) follows from the fact that $N_{g}\geq 1$ almost surely, and (3)
follows from Assumption 1.1 with $\alpha>2$.

Now we show the consistency of the variance estimator, that is,

$G\widehat{V}_{\hat{\theta}}^{\text{WCR}}\xrightarrow{p}V^{\text{WCR}}.$

The convergence
$G^{-1}\sum_{g=1}^{G}N_{g}^{-1}X_{g}^{\prime}X_{g}\overset{p}{\rightarrow}\mathbb{E}[N_{g}^{-1}X^{\prime}_{g}X_{g}]$
follows from the law of large numbers given (8) and Assumption 2.1. Given that
$a_{n}\rightarrow 1$ almost surely, it remains to show that

$G^{-1}\sum_{g=1}^{G}N_{g}^{-2}\widehat{S}_{g}\widehat{S}_{g}^{\prime}\overset{p}{\rightarrow}\mathbb{E}[N_{g}^{-2}X^{\prime}_{g}X_{g}U^{2}_{g}].$

To this end, we first write

$\displaystyle\widehat{S}_{g}=$
$\displaystyle\,\sum_{i=1}^{N_{g}}X_{gi}\widehat{U}_{gi}$ $\displaystyle=$
$\displaystyle\,\sum_{i=1}^{N_{g}}X_{gi}U_{gi}-\sum_{i=1}^{N_{g}}X_{gi}^{2}\,(\widehat{\theta}-\theta)$
$\displaystyle=$ $\displaystyle\,S_{g}-M_{g}\,(\widehat{\theta}-\theta),$

where we again illustrate with the scalar $X$ case. Then, we have

$\displaystyle
G^{-1}\sum_{g=1}^{G}N_{g}^{-2}\widehat{S}_{g}\widehat{S}_{g}^{\prime}$
$\displaystyle=$
$\displaystyle\,G^{-1}\sum_{g=1}^{G}N_{g}^{-2}S_{g}^{2}-2G^{-1}\sum_{g=1}^{G}N_{g}^{-2}S_{g}M_{g}(\widehat{\theta}-\theta)$
$\displaystyle\,+G^{-1}\sum_{g=1}^{G}N_{g}^{-2}M_{g}^{2}(\widehat{\theta}-\theta)^{2}$
$\displaystyle\overset{(1)}{=}$
$\displaystyle\,\mathbb{E}[N_{g}^{-2}S_{g}^{2}]-2\mathbb{E}[N_{g}^{-2}S_{g}M_{g}](\widehat{\theta}-\theta)+\mathbb{E}[N_{g}^{-2}M_{g}^{2}](\widehat{\theta}-\theta)^{2}+o_{p}(1)$
$\displaystyle\overset{(2)}{=}$
$\displaystyle\,\mathbb{E}[N_{g}^{-2}S_{g}^{2}]+o_{p}(1),$

where (1) follows from the law of large numbers given (8), (9), Assumption 1.1
and Assumption 2.1, and (2) follows from the consistency of
$\widehat{\theta}$. ∎

##  Appendix B Testing the Robustness

As discussed in Section 4, essential for credible CR variance estimation under
cluster sampling based on existing theory is a finite second moment of
$S_{g}$. Thus, we are interested in testing the hypotheses

$\displaystyle H_{0}:\mathbb{E}[||S_{g}||^{2}]<\infty$ (10)

Note that a naïve application of common tail index or Pareto exponent
estimators will not work to test this hypothesis, since the score $S_{g}$ in
general involves unknown parameters that need to be estimated from data.
However, this hypothesis (10) can be still formally tested by applying some
existing theories and methods. This appendix section briefly describes the
testing procedure.

In general, we are interested in testing the hypothesis

$\displaystyle H_{0}:\mathbb{E}[||S_{g}||^{r}]<\infty\qquad$ (11)

It can be tested by applying a modified version of the method of Sasaki and
Wang (2022) suitably adapted to the cluster sampling framework. A proposed
testing procedure is outlined below.

Define
$\widehat{A}^{r}_{g}=({\widehat{S}_{g}}^{\prime}{\widehat{S}_{g}})^{r/2}$ for
$r\in\\{1,2\\}$. Sort $\\{\widehat{A}^{r}_{g}\\}_{g=1}^{G}$ in the descending
order into
$\widehat{A}^{r}_{(1)}\geq\widehat{A}^{r}_{(2)}\geq\cdots\geq\widehat{A}^{r}_{(G)}$,
i.e., $\widehat{A}^{r}_{(g)}$ denotes the $g$-th order statistic of
$\\{\widehat{A}^{r}_{g}\\}_{g=1}^{G}$ for each $g\in\\{1,\cdots,G\\}$. For a
predetermined integer $k\geq 3$, collect the largest $k$ order statistics
$\widehat{\mathbf{A}}^{r}=(\widehat{A}^{r}_{(1)},\cdots,A^{r}_{(k)})^{\prime}$.
Implement the location- and scale-normalization of $\widehat{\mathbf{A}}^{r}$
to obtain

$\widehat{\mathbf{A}}^{r}_{\ast}=\frac{\widehat{\mathbf{A}}^{r}-\widehat{A}^{r}_{(k)}}{\widehat{A}^{r}_{(1)}-\widehat{A}^{r}_{(k)}}.$

We reject the null hypothesis (11) of a finite $r$-th moment of $S_{g}$ if

$\frac{\int_{1}^{2}f_{{\mathbf{V}}_{\ast}}(\widehat{\mathbf{A}}^{r}_{\ast};\xi)dW(\xi)}{\int_{0}^{1}f_{{\mathbf{V}}_{\ast}}(\widehat{\mathbf{A}}^{r}_{\ast};\xi)d\Lambda(\xi)}>1,$

where $W$ denotes the uniform probability measure, $\Lambda$ is a pre-
determined measure that transforms the composite null into a simple one while
subsuming the critical value for size control, and $f_{{\mathbf{V}}_{\ast}}$
is defined by

$f_{{\mathbf{V}}_{\ast}}({\mathbf{v}}_{\ast};\xi)=\Gamma(k)\int_{0}^{\infty}s^{k-2}\exp\left((-1-\xi^{-1})\left(\log(1+\xi
s)+\sum_{i=2}^{k-1}\log(1+\xi v_{\ast i}s)\right)\right)ds$

with $\Gamma$ denoting the gamma function. See Sasaki and Wang (2022) for
details. The Stata command ‘testout’ associated with this article implements
this test.333This Stata command is available through SSC by the command line:
ssc install testout.