# Re-entrant magic-angle phenomena in twisted bilayer graphene
in integer magnetic fluxes

Yifei Guan  Institute of Physics, Ecole Polytechnique Fédérale de Lausanne
(EPFL), Lausanne, CH-1015, Switzerland  Oleg V. Yazyev  Institute of Physics,
Ecole Polytechnique Fédérale de Lausanne (EPFL), Lausanne, CH-1015,
Switzerland  Alexander Kruchkov  Institute of Physics, Ecole Polytechnique
Fédérale de Lausanne (EPFL), Lausanne, CH-1015, Switzerland  Department of
Physics, Harvard University, Cambridge, MA 02138, USA  Branco Weiss Society in
Science, ETH Zurich, Zurich, Switzerland

###### Abstract

In this work we address the re-entrance of magic-angle phenomena (band
flatness and quantum-geometric transport) in twisted bilayer graphene (TBG)
subjected to strong magnetic fluxes $\pm\Phi_{0}$, $\pm 2\Phi_{0}$, $\pm
3\Phi_{0}$… ($\Phi_{0}=h/e$ is the flux quantum per moiré cell). The moiré
translation invariance is restored at the integer fluxes, for which we
calculate the TBG band structure using accurate atomistic models with lattice
relaxations. Similarly to the zero-flux physics outside the magic angle
condition, the reported effect breaks down rapidly with the twist. We conclude
that the magic-angle physics re-emerges in high magnetic fields, witnessed by
the appearance of flat electronic bands distinct from Landau levels, and
manifesting non-trivial quantum geometry. We further discuss the possible
flat-band quantum geometric contribution to the superfluid weight in strong
magnetic fields (28 T at 1.08∘ twist), according to Peotta-Törmä mechanism.

In 2D systems, the electronic spectrum in magnetic field develops a fractal
structure (”Hofstadter butterfly”),Hofstadter (1976) which lowers the
effective dimensionality, and contributes to suppressing superconductivity
(long-range order is generically destroyed in dimensions lower than $D$=$2$).
The observation of Hofstadter physics requires strong magnetic fluxes
($\sim$$h/e$), which became experimentally accessible only with the advent of
moiré superlattices.Dean  _et al._ (2013); Hunt  _et al._ (2013); Ponomarenko
_et al._ (2013) In those experiments,Dean  _et al._ (2013); Hunt  _et al._
(2013); Ponomarenko  _et al._ (2013) the magnetic fields of nearly 30 T were
employed in the system of graphene monolayer twisted on hexagonal boron
nitride (hBN), resulting into effective fluxes of $\Phi=B/A\sim\Phi_{0}$ ($A$
is the moiré cell area, $\Phi_{0}$$=$$h/e=$$4$ Wb is magnetic flux quantum).
Furthermore, the twisted graphene multilayers provide a natural platform to
test the interplay between the Hofstadter physics and strong correlations.Cao
_et al._ (2018); Hao  _et al._ (2021); Park  _et al._ (2021a, b); Zhang  _et
al._ (2021) In twisted bilayer graphene, the smaller is the twist $\theta$,
the larger is the effective magnetic flux ($\Phi\propto 1/\theta^{2}$) at the
fixed field $\mathbf{B}$: for TBG at the magic angle $1.08^{\circ}$ the
magnetic flux quantum corresponds to $B_{0}\approx 28$ T, which is reachable
in the modern laboratories.Hahn  _et al._ (2019)

The magic-angle graphene heterostructures—two or more graphene sheets twisted
to the angle $\sim$$1^{\circ}$, at which a very narrow band emerges in the
electronic spectrum—have re-attracted significant attention due to re-entrant
superconductivity in strong magnetic fields and reported Pauli limit
violation.Cao  _et al._ (2021); Chaudhary  _et al._ (2021); Shaffer  _et al._
(2021) The re-entrant correlated (Chern) insulator phases were reported at
strong magnetic fluxes,Das  _et al._ (2021); Herzog-Arbeitman  _et al._ (2021)
close to the unit magnetic flux quantum $\Phi_{0}=h/e$ per moiré unit cell
(see also Ref.Sheffer and Stern (2021)). In this paper we show that in all the
integer magnetic flux $\Phi=\pm N\Phi_{0}$, the magic angle physics of TBG is
re-entrant and similar to the physics of TBG in zero magnetic field. Namely,
we find non-trivial flat Chern bands at integer flux $h/e$ (see Fig. 1). These
bands are distinct from Landau levels for three reasons: (i) Outside the magic
angle range the band flatness breaks down; (ii) The magic-angle flat bands are
characterized by the total Chern number $|C|=2$, while Landau levels are
characterized by $|C|=1$; and (iii) The quantum geometry (Fubini-Study metrics
$\mathfrak{G}_{ij}$) of the magic-angle flat bands is highly nontrivial and
different from the Landau level case, and provides the enhanced basis for
quantum-geometric transportPeotta and Törmä (2015) in flat bands without
dispersive contributions due to
$\sum_{\mathbf{k}}\text{Tr}\,\mathfrak{G}_{ij}(\mathbf{k})\approx 3$.

Figure 1: Electronic band structure of magic-angle twisted bilayer graphene in
integer magnetic flux: (a) zero flux ($\Phi=0$); (b) flux one ($\Phi=h/e$);
(c) flux two ($\Phi=2h/e$). All band structures are calculated with tight
binding model including lattice relaxations effects. In integer magnetic flux,
the lower (two out four) flat bands acquire $C=-2$.

In this paper we report that the nontrivial quantum geometry, electronic band
flatness, and conditions for unconventional quantum transport are re-
established at integer magnetic flux (in units of $h/e$) through the moiré
unit cell. Similarly to the zero-flux case, the Fermi velocity and the
bandwidth drop dramatically at the magic angle, to reappear as a dispersive
band both below and above the magic angle. Importantly, the flat bands in
finite flux has non-trivial Chern numbers $|C|$=$2$ (defined at half filling,
see Fig. 1) and non-trivial quantum geometry, which follows the quantum-
geometric flatness criterionKruchkov (2021)

$\displaystyle\text{Tr}\mathcal{G}_{ij}(\mathbf{k})\simeq\mathcal{F}_{xy}(\mathbf{k}).$
(1)

We show numerically with atomistic calculations that in the realistic TBG, the
”ideal band condition” (1) is satisfied in the magnetic moiré Brillouin zone
(mmBZ) regions where the band flatness is pronounced in terms of vanishing
Fermi velocity (e.g. around the K points of mmBZ). Here $\mathcal{G}_{ij}$ and
$F_{xy}(\mathbf{k})$ are the real and imaginary parts of the quantum-geometric
tensor $\mathfrak{G}_{ij}$(see further), determining the quantum distance
between electronic states in the (projected) Hilbert space.Provost and Vallee
(1980)

In this paper, we investigate the effect of strong magnetic fields with
integer flux in twisted bilayer graphene with the help of the accurate
atomistic model including lattice relaxation effects at the magic angle, and
compare the observed results with the established knowledge of the zero-flux
TBG case. The details of the tight-binding Hamiltonian are provided in
Supplementary Materials (SM).SI  The key observation is that the magnetic
translation operators commute at every integer flux $\Phi=N\Phi_{0}$, namely

$\displaystyle\hat{T}_{\mathbf{a}_{1}}\hat{T}_{\mathbf{a}_{2}}=e^{i2\pi\Phi/\Phi_{0}}\hat{T}_{\mathbf{a}_{2}}\hat{T}_{\mathbf{a}_{1}},\
\ \to\ \
[\hat{T}_{\mathbf{a}_{1}},\hat{T}_{\mathbf{a}_{2}}]_{\Phi=N\Phi_{0}}=0.$

Thus the moiré unit cell is restored, and the system flows towards the
electronic band structure defined on the mmBZ, which in the integer flux has
the same periodicity as the moiré Brillouin zone (mBZ) in zero magnetic flux.
But instead of dispersionless Landau levels, we recover the set of dispersive
bands, with its band structure depending crucially on the twist angle (Fig.3).

Figure 2: Hofstader spectrum of magic-angle bilayer graphene in magnetic flux
(in units $h/e$). The digits in gaps indicate the in-gap Chern numbers,
calculated through the edge modes counting. Additionally, the insets below
shows Wannier charge center (WCC) winding at fluxes $h/e$ and $2h/e$,
calculated at the half filling of the flat bands (2 of 4 subbands are
occupied). The WCC winding is nontrivial, revealing $|C|$$=$$2$ at half
filling in the integer flux.

Re-entrant magic angle spectra at integer magnetic flux. We start from
considering the electronic band structure at an integer magnetic flux
$\Phi=Nh/e$. Such a flux provides the $2\pi N$ circulation of magnetic vector
potential, and hence reconstructs the moire Brillouin zone (mBZ); magnetic
translation operators commute at integer flux
$[\hat{T}_{\mathbf{a}_{1}},\hat{T}_{\mathbf{a}_{2}}]=0$. This allows us to re-
introduce momentum as a good quantum number and compute the electronic band
structure starting from the tight binding model for magic-angle twisted
bilayer graphene with atomic relaxations, modified with Peierls substitution,

$\displaystyle t_{ij}\to
t_{ij}e^{-i\frac{e}{h}\int_{\mathbf{r}_{j}}^{\mathbf{r}_{i}}\mathbf{A}(\mathbf{r}^{\prime})d\mathbf{r}^{\prime}}.$
(2)

Comparing to the widely-used continuum models,Dos Santos  _et al._ (2007);
Bistritzer and MacDonald (2011a); Tarnopolsky  _et al._ (2019) the accurate
tight binding model has an important advantage of addressing the magic-angle
physics under realistic conditions of atomic lattice relaxations, proved to be
indispensable in the experiments due to domain formation.Carr  _et al._ (2018)
Worth noting, the finite magnetic flux shifts the effective Brillouin zone
(see Fig. 1), thus re-defining the positions of high-symmetry points.
Otherwise, the original moire BZ and reconstructed magnetic moire Brillouin
zone (mmBZ) have the same orientation and periodicity, which is the main
technical condition to observe the magic angle phenomena.

We report that the characteristic flat band, the hallmark of the magic-angle
graphene, re-appears in every integer magnetic flux, $\pm\Phi_{0}$, $\pm
2\Phi_{0}$, $\pm 3\Phi_{0}$… Figure 1 provides the electronic band spectrum at
$\Phi=0,\Phi_{0},2\Phi_{0}$, while the band structures (and other properties)
in stronger magnetic fluxes are given in the SM.SI  The first observation is
that the flat band reappears exactly at the magic angle, while the higher
bands are dispersive (see Fig.3). We argue below that the magic-angle flat
band (MAFB) in the integer flux are not a consequence of Landau level (LL)
flattening since: (i) It has $|C|$=$2$ at half-filling; (ii) It demonstrates
quantum geometry incompatible with LL physics (Fig.5); and (iii) It becomes
dispersive outside the magic angle. Worth noting, the strong magnetic fields
restore the asymptotic particle-hole symmetry of the low energy states, which
was moderately broken in the zero flux.

Figure 3:  Magic angle signature in integer magnetic flux $\Phi=h/e$ at
different twists. (a) Above the magic angle, $\theta=1.47^{\circ}$. (b) At the
magic angle, $\theta=1.08^{\circ}$. (c) Below the magic angle,
$\theta=0.81^{\circ}$. We observe the similar behavior as the zero-flux TBG
tuned outside of the magic angle.Bistritzer and MacDonald (2011b); Tarnopolsky
_et al._ (2019)

Magnetic spectrum distinct from Landau levels. We further address the
properties of the electronic band spectrum versus magnetic field. For this,
the characteristic quantity to calculate is the Hofstadter diagram, which
traces the energy of electronic states allowed in quantized magnetic field, as
a function of magnetic flux through the unit cell (Fig. 2). We observe that
the flat bands at the integer flux are not stemming from Landau level
Hofstadter physics, but rather from the magic angle physics of the TBG. To
show that they are fundamentally different from the LLs, we calculate the
Chern number through Wilson loop computation in the form of Wannier charge
center winding, see Fig. 2. We find that the flat bands at integer $\Phi$ have
Chern numbers $|C|=2$ incompatible with Chern numbers of Landau levels on the
lattice ($|C|=1$). Furthermore, we fix magnetic flux to $h/e$, and investigate
the change in electronic band spectrum versus the change in twist angle (see
Fig. 3; Cases of higher integer flux are provided in SM.SI ). While the
spectrum is strongly dispersive outside of the magic angle (the bandwidth is
approximately 100 meV at $1.47^{\circ}$ twist, at the magic angle
$1.08^{\circ}$ the bandwidth is just $15$ meV, comparable to the magic angle
bandwidth in zero flux (see also Fig. 1). We conclude that the magic angle
physics is in its essence restored.

Figure 4:  Asymptotic particle-hole symmetry (PHS) in strong magnetic flux is
emergent for $\Phi$$\geq$$2h/e$. The low-energy PH asymmetry of the flat bands
is determined as a maximum energy difference between the conductance and
valence bands at neutrality, in meV. The PH asymmetry is maximal in zero flux
($\Phi=0$) resulting in nearly 4 meV. The PHS violation is however strongly
suppressed with applied magnetic field, resulting to 1 meV at $\Phi=h/e$, and
just 0.07 meV at $\Phi=2h/e$. The suppression of PH asymmetry roughly follows
$\sim\exp[-\Phi/\Phi_{*}]$, with $\Phi_{*}\approx 0.7h/e$.

Emergent particle-hole symmetry. We now focus on the low-energy numerical
analysis of the flat bands. The main observations under strong magnetic fields
is that the low-energy electronic spectrum, featuring flat bands, acquires
asymptotic particle-hole symmetry (Fig. 4) To quantify the particle-hole
asymmetry (PHA), we track the energy difference maximum difference in the flat
band energies between the occupied and free bands at neutrality, and
investigate it versus magnetic field (Fig. 4). We observe that while PHS is
generically violated in zero flux (PHA = 4 meV), in strong magnetic fields
this symmetry reemerges in its asymptotic form above $\Phi\geq 2h/e$ (Fig. 4).
Approximately, the suppression of PH asymmetry is observed as
$\sim\exp[-\Phi/\Phi_{*}]$, with $\Phi_{*}\approx 0.7h/e$ (corresponding to
$\approx 20$ T). This could be qualitatively understood in terms local tight
binding hoppings, which due to strong magnetic fields oscillate rapidly in
space; the system thus performs self-averaging which re-defines effective
hopping parameters. In realistic TBG systems at $B=0$, with sublattice
hoppings and lattice relaxations, the Chern number in zero flux vanishes since
the particle-hole symmetry (PHS) is violated on the atomistic level, chiral
symmetry (CS) is broken explicitly, and time reversal (TRS) is present. At
integer fluxes, we have asymptotic PHS, broken TRS, broken CS, and while
strictly speaking the system belongs to class $A$ topological insulator, its
dynamics flows towards class $C$, characterized by $2Z$ (even-valued)
topological invariants in 2D systems.Altland and Zirnbauer (1997) This gives a
plausible explanation of promotion of the $|C|$=$2$ Chern numbers in the flat
bands at integer flux $|\Phi|\geq\Phi_{0}$, once the two of four subbands are
slightly gapped out (Fig. 1). We observe that the magic-angle TBG at zero flux
$\Phi=0$ and at integer fluxes $\Phi=\pm N\Phi_{0}$ belong to different
topological classes, which should be taken into account for understanding
recent TBG experiments.

Non-trivial quantum geometric properties. We report that the magic-angle flat
band in integer magnetic flux has non-trivial quantum geometry (Fig. 4),
distinct from Landau levels. The quantum geometry is defined for the Bloch
states in the projective Hilbert space; it can be separated into real
(diagonal) part which is Fubini-Study metrics, and imaginary (off-diagonal)
part, with its components being Berry curvature.Provost and Vallee (1980)
Numerically, we compute quantum-geometric tensor for flat bands in TBG by
using its spectral representationKruchkov (2022)

$\displaystyle\mathfrak{G}_{ij}(\mathbf{k})=\sum_{n,m}\frac{\langle
u_{n\mathbf{k}}|\frac{\partial\mathcal{H}_{\mathbf{k}}}{\partial
k_{i}}|u_{m\mathbf{k}}\rangle_{0}\langle
u_{m\mathbf{k}}|\frac{\partial\mathcal{H}_{\mathbf{k}}}{\partial
k_{i}}|u_{n\mathbf{k}}\rangle_{0}}{(\varepsilon_{n\mathbf{k}}-\varepsilon_{m\mathbf{k}})^{2}}.$
(3)

We further introduce $\mathcal{G}_{ij}$=$\text{Re}\mathfrak{G}_{ij}$,
$\mathcal{F}_{ij}$=$-2\,\text{Im}\mathfrak{G}_{ij}$. One can find a basis in
which $\mathcal{G}_{ij}$ is diagonal and $\mathcal{F}_{ij}$ is off-diagonal.
The plots for $\mathcal{G}_{xx},\mathcal{G}_{yy}$ (Fubini-Study metrics)
calculated at half-filling are presented in Fig. 5. We observe that the re-
entrant flat band has nontrivial quantum geometry within the mmBZ, as
manifested in Figs. 5(a,b,c), which is not compatible with a Landau level
quantum geometry (the LL quantum geometry is constant in the whole Brillouin
zone).

For comparison, we plot the Berry curvature $\mathcal{F}_{xy}$ together with
trace of Fubini-Study’s $\mathcal{G}_{ij}$ in Fig. 3d. We observe that the
flat band closely follows the quantum-geometric condition for ideal flat
bandsKruchkov (2021) $\text{Tr}\mathcal{G}_{ij}$=$\mathcal{F}_{xy}$. It is
certainly interesting that this condition is satisfied almost exactly in the
regions of mmBZ, where the band flatness is pronounced in terms of vanishing
Fermi velocity (around the K points of mmBZ). The deviation to this quantum-
geometric bound $\text{Tr}\mathcal{G}_{ij}$=$\mathcal{F}_{xy}$ are observed in
the regions of mmBZ with finite dispersion and significant
$v_{k}=\partial{\varepsilon_{k}}/\partial{k}$ caused by broken CS of the tight
binding calculations. The criterion $\text{Tr}\mathcal{G}_{ij}=F_{xy},$ tests
the closeness of a realistic flat band in TBG to flat band idealization
through holomorphic/meromorphic representation of the flat band wave
functions.Kruchkov (2021) However, since total $|C|=2$, the relevant toy model
for TBG in integer flux cannot be represented by solely a holomorphic
representation of the quasi-LLL TBG (found in Ref.Tarnopolsky  _et al._
(2019)); one needs to take meromorphic flat band contributions into
account.Popov and Milekhin (2021)

We calculate the distribution of the Berry curvature $F_{xy}$ within the mBZ
(the cut is shown in Fig. 3d, the density plots are in SM), which reveals non-
homogeneous structure, which is not consistent with the homogeneous Berry
curvature of the generic Landau levels (LLs are ”Berry flat”). This yet again
provides arguments in the support of the magic-angle physics re-entrance in
integer magnetic flux. Finally, we check numerically that the Berry flux
$F_{xy}$ encompassed by mmBZ sums up to $|C|$=$2.0\pm 0.07$ at the half
filling (see Figs. 1,2).

Quantum-geometric transport. As was first introduced by Peotta and TörmäPeotta
and Törmä (2015), the nontrivial quantum-geometric tensor $\mathfrak{G}_{ij}$
leads to the the finite superfluid current $J_{i}$$=$$-D_{ij}A_{j}$ even in
the limit of perfectly flat band (here $D_{ij}$ is the superfluid weight).
This argument is now understood to apply directly to TBG in zero flux, where
the $\text{Tr}\mathcal{G}_{ij}$ is nonzero due to hidden nontrivial topology
of the flat bands. It was reported with different methodsHu  _et al._ (2019);
Xie  _et al._ (2020); Julku  _et al._ (2020) that the quantum geometric tensor
(QGT) contribution to the superfluid weight $D_{S}$ in TBG is if not dominant,
than at least commensurate with the conventional contributions, thus leading
to the BTK transition temperature estimate
$T_{\text{BKT}}\sim\sum_{\mathbf{k}}\text{Tr}\mathcal{G}_{ij}$. The QGT
contribution holds for different symmetries of the order parameter, and the
argument is valid beyond the mean field.Wang  _et al._ (2020) The essential
physics is captured by Bogoliubov-de-Gennes Hamiltonian

$\displaystyle\mathcal{H}_{\text{BdG}}=\left(\begin{array}[]{cc}H_{\mathbf{k}}&\Delta_{\mathbf{k}}\\\
\Delta^{*}_{-\mathbf{k}}&H^{*}_{-\mathbf{k}}\end{array}\right)$ (6)

Without loss of generality, we consider superconducting order
$\Delta_{\mathbf{k}}=\Delta$. The superfluid weight is then calculated within
Kubo formalism through the current-current correlators. We explicitly
calculate $\mathcal{G}_{ij}$ numerically (Fig. 4) with consequent mmBZ
integration at half filling (for two of four flat bands occupied, we have
$\text{Tr}\mathcal{G}_{ij}\approx{2.8}$), to obtain at the superfluid weight
maximum positioned at the middle of the composite flat band with $|C|=2$,

$\displaystyle
D_{xx}=\frac{2e^{2}}{\hbar^{2}}\Delta\sum_{\mathbf{k}}\text{Tr}\mathcal{G}_{ij}(\mathbf{k})\approx(5.6\pm
0.1)\frac{e^{2}}{\hbar^{2}}\Delta.$ (7)

Symmetry $D_{xx}$$=$$D_{yy}$ is assumed. Here we took into account factor
$\sqrt{\nu(1-\nu)}$ (in notations of Ref.Peotta and Törmä (2015)), where $\nu$
is the filling factor of the composite flat band (indexed by $C=-2$ in Fig.
1).

We can further make estimates for the BKT transition temperature,Berezinskii
(1971); Kosterlitz and Thouless (1973); Nelson and Kosterlitz (1977)
indicating the disappearing of the phase coherence of superconducting order
from expression $\pi\hbar^{2}D(T_{*})/8e^{2}T_{*}=1$. The order-of-magnitude
estimate gives $T_{*}\sim{\hbar^{2}D(0)}/{e^{2}}\sim\Delta$. The remaining
question is of course, what is the value of $\Delta$, which should be found
self-consistently by solving Gorkov equations in magnetic field, or through
indirect experimental data—and it is beyond the scope of this paper. For a
rough estimate, even $\Delta\sim 0.1$ meV will give a physically relevant
$T_{*}\sim 1$ K. Whether or not the superfluid order is re-entrant in strong
magnetic flux remains open,Chaudhary  _et al._ (2021); Shaffer  _et al._
(2021) however it would not be surprising in the view of the reported Pauli-
limit violation and re-entrant superconductivity in strong magnetic parallel
fields in twisted graphene multilayers.Cao  _et al._ (2021); Park  _et al._
(2021b); Zhang  _et al._ (2021) Note the magnetic field may or may not change
the symmetry of the superconducting order (resulting into different
$\Delta_{\mathbf{k}}$), however the quantum-geometric superfluid weight will
still remain finite due to the nontrivial geometry of the underlying flat
bands.

Figure 5: Nontrivial quantum metrics in integer magnetic flux $\Phi=h/e$,
calculated at half-filling (two of four flat bands are occupied). (a-c)
components of Fubini-Study metrics; (c) shows that $\text{Tr}\mathcal{G}_{ij}$
is gauge-invariant (d) Comparison between trace of Fubini-Study tensor
$\text{Tr}\mathcal{G}_{ij}$ and Berry curvature $\mathcal{F}_{xy}$, which
probes how close the band flatness is to the perfect flatness through Eq.(1).

Conclusions. We conclude that there is a re-entrant magic angle physics in
twisted bilayer graphene at every integer magnetic flux quanta $\pm\Phi_{0}$,
$\pm 2\Phi_{0}$, $\pm 3\Phi_{0}$, etc., through the moiré cell. To date, the
practical importance represents the first magnetic quantum $\pm\Phi_{0}$,
which at twist angle $1.08^{\circ}$ corresponds to experimentally-achievable
fields of 28 Tesla. We confirm with accurate atomistic calculations,
incorporating lattice relaxation effects, that at such fields the magic-angle
phenomena re-emerge. This, in particular, could be seen through the re-
emergence of very flat bands at the magic angle distinct from Landau levels,
while beyond the magic angle this physics breaks down, similar to the zero-
flux case.Bistritzer and MacDonald (2011b); Tarnopolsky  _et al._ (2019) These
flat bands at half filling carry nontrivial Chern numbers ($|C|$=$2$) and
nontrivial quantum geometry (Fubini-Study metrics), and are fundamentally
different from conventional Landau levels. We conjecture that, similar to TBG
in the zero flux, there is a nonvanishing contribution to the superfluid
weight coming from the re-entrant quantum geometric properties, and in the
flat topological bands this contribution is significant. Due to the strong
quantum geometry of the TBG flat bands in integer flux
($\sum_{\text{BZ}}\text{Tr}\mathcal{G}_{ij}\approx 3$), and the estimated BKT
temperature is in order of the gap $T_{*}\approx\Delta$, which gives values
significantly elevated with regard to the conventional superconductivity in
geometrically-trivial dispersive bands. The behavior of superfluid order
parameter beyond the conventional Pauli limit (towards integer magnetic flux)
is a subject for further research.

Acknowledgments. The authors thank Ady Stern, Yarden Sheffer, Luiz Santos,
Gaurav Chaudhary, and Marta Brzezińska for useful discussion. The authors
thank B. Andrei Bernevig, Jonah Herzog-Arbeitman, Aaron Chew for further
discussion. The project was supported by the Branco Weiss Society in Science,
ETH Zurich, through the research grant on flat bands, strong interactions and
SYK physics, and by the Swiss National Science Foundation, grants No. 172543.
Computations have been performed at the Swiss National Supercomputing Centre
(CSCS) under project s1008 and the facilities of Scientific IT and Application
Support Center of EPFL.

## References

  * Hofstadter (1976) D. R. Hofstadter, Physical Review B 14, 2239 (1976).
  * Dean  _et al._ (2013) C. R. Dean, L. Wang, P. Maher, C. Forsythe, F. Ghahari, Y. Gao, J. Katoch, M. Ishigami, P. Moon, M. Koshino, _et al._ , Nature 497, 598 (2013).
  * Hunt  _et al._ (2013) B. Hunt, J. D. Sanchez-Yamagishi, A. F. Young, M. Yankowitz, B. J. LeRoy, K. Watanabe, T. Taniguchi, P. Moon, M. Koshino, P. Jarillo-Herrero, _et al._ , Science 340, 1427 (2013).
  * Ponomarenko  _et al._ (2013) L. Ponomarenko, R. Gorbachev, G. Yu, D. Elias, R. Jalil, A. Patel, A. Mishchenko, A. Mayorov, C. Woods, J. Wallbank, _et al._ , Nature 497, 594 (2013).
  * Cao  _et al._ (2018) Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43 (2018).
  * Hao  _et al._ (2021) Z. Hao, A. Zimmerman, P. Ledwith, E. Khalaf, D. H. Najafabadi, K. Watanabe, T. Taniguchi, A. Vishwanath, and P. Kim, Science 371, 1133 (2021).
  * Park  _et al._ (2021a) J. M. Park, Y. Cao, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Nature 590, 249 (2021a).
  * Park  _et al._ (2021b) J. M. Park, Y. Cao, L. Xia, S. Sun, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, arXiv preprint arXiv:2112.10760 (2021b).
  * Zhang  _et al._ (2021) Y. Zhang, R. Polski, C. Lewandowski, A. Thomson, Y. Peng, Y. Choi, H. Kim, K. Watanabe, T. Taniguchi, J. Alicea, _et al._ , arXiv preprint arXiv:2112.09270 (2021).
  * Hahn  _et al._ (2019) S. Hahn, K. Kim, K. Kim, X. Hu, T. Painter, I. Dixon, S. Kim, K. R. Bhattarai, S. Noguchi, J. Jaroszynski, _et al._ , Nature 570, 496 (2019).
  * Cao  _et al._ (2021) Y. Cao, J. M. Park, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Nature 595, 526 (2021).
  * Chaudhary  _et al._ (2021) G. Chaudhary, A. MacDonald, and M. Norman, Physical Review Research 3, 033260 (2021).
  * Shaffer  _et al._ (2021) D. Shaffer, J. Wang, and L. H. Santos, Physical Review B 104, 184501 (2021).
  * Das  _et al._ (2021) I. Das, C. Shen, A. Jaoui, J. Herzog-Arbeitman, A. Chew, C.-W. Cho, K. Watanabe, T. Taniguchi, B. A. Piot, B. A. Bernevig, _et al._ , arXiv preprint arXiv:2111.11341 (2021).
  * Herzog-Arbeitman  _et al._ (2021) J. Herzog-Arbeitman, A. Chew, D. K. Efetov, and B. A. Bernevig, arXiv preprint arXiv:2111.11434 (2021).
  * Sheffer and Stern (2021) Y. Sheffer and A. Stern, Physical Review B 104, L121405 (2021).
  * Peotta and Törmä (2015) S. Peotta and P. Törmä, Nature Communications 6, 1 (2015).
  * Kruchkov (2021) A. Kruchkov, arXiv preprint arXiv:2105.14672 (2021).
  * Provost and Vallee (1980) J. P. Provost and G. Vallee, Communications in Mathematical Physics 76, 289 (1980).
  * (20) Supplementary Materials (SM) are available in the published version of this manuscript .
  * Dos Santos  _et al._ (2007) J. L. Dos Santos, N. Peres, and A. C. Neto, Physical Review Letters 99, 256802 (2007).
  * Bistritzer and MacDonald (2011a) R. Bistritzer and A. H. MacDonald, Proceedings of the National Academy of Sciences 108, 12233 (2011a).
  * Tarnopolsky  _et al._ (2019) G. Tarnopolsky, A. J. Kruchkov, and A. Vishwanath, Physical Review Letters 122, 106405 (2019).
  * Carr  _et al._ (2018) S. Carr, D. Massatt, S. B. Torrisi, P. Cazeaux, M. Luskin, and E. Kaxiras, Physical Review B 98, 224102 (2018).
  * Bistritzer and MacDonald (2011b) R. Bistritzer and A. H. MacDonald, Proceedings of the National Academy of Sciences 108, 12233 (2011b).
  * Altland and Zirnbauer (1997) A. Altland and M. R. Zirnbauer, Physical Review B 55, 1142 (1997).
  * Kruchkov (2022) A. Kruchkov, arXiv:2XXX.XXXXX (2022).
  * Popov and Milekhin (2021) F. K. Popov and A. Milekhin, Physical Review B 103, 155150 (2021).
  * Hu  _et al._ (2019) X. Hu, T. Hyart, D. I. Pikulin, and E. Rossi, Physical Review Letters 123, 237002 (2019).
  * Xie  _et al._ (2020) F. Xie, Z. Song, B. Lian, and B. A. Bernevig, Physical Review Letters 124, 167002 (2020).
  * Julku  _et al._ (2020) A. Julku, T. J. Peltonen, L. Liang, T. T. Heikkilä, and P. Törmä, Physical Review B 101, 060505 (2020).
  * Wang  _et al._ (2020) Z. Wang, G. Chaudhary, Q. Chen, and K. Levin, Physical Review B 102, 184504 (2020).
  * Berezinskii (1971) V. Berezinskii, Sov. Phys. JETP 34, 610 (1971).
  * Kosterlitz and Thouless (1973) J. M. Kosterlitz and D. J. Thouless, Journal of Physics C: Solid State Physics 6, 1181 (1973).
  * Nelson and Kosterlitz (1977) D. R. Nelson and J. Kosterlitz, Physical Review Letters 39, 1201 (1977).


# Confirmation of the centrality of the Huanan market among early
COVID-19 cases
Reply to Stoyan and Chiu (2024)

F. Débarre1 & M. Worobey2

1 Institut d’Écologie et des Sciences de l’Environnement (IEES-Paris, UMR
7618), CNRS, Sorbonne Université, UPEC, IRD, INRAE, Paris, France. ORCID:
0000-0003-2497-833X.
Contact<EMAIL_ADDRESS>
2 Department of Ecology and Evolutionary Biology, University of Arizona,
Tucson, AZ, USA.
Contact<EMAIL_ADDRESS>
## Abstract

The centrality of Wuhan’s Huanan market in maps of December 2019 COVID-19 case
residential locations, established by Worobey et al. (2022a), has recently
been challenged by Stoyan and Chiu (2024, SC2024). SC2024 proposed a
statistical test based on the premise that the measure of central tendency
(hereafter, “centre”) of a sample of case locations must coincide with the
exact point from which local transmission began. Here we show that this
premise is erroneous. SC2024 put forward two alternative centres (centroid and
mode) to the centre-point which was used by Worobey et al. for some analyses,
and proposed a bootstrapping method, based on their premise, to test whether a
particular location is consistent with it being the point source of
transmission. We show that SC2024’s concerns about the use of centre-points
are inconsequential, and that use of centroids for these data is inadvisable.
The mode is an appropriate, even optimal, choice as centre; however, contrary
to SC2024’s results, we demonstrate that with proper implementation of their
methods, the mode falls at the entrance of a parking lot at the market itself,
and the $95\%$ confidence region around the mode includes the market. Thus,
the market cannot be rejected as central even by SC2024’s overly stringent
statistical test. Our results directly contradict SC2024’s and – together with
myriad additional lines of evidence overlooked by SC2024, including crucial
epidemiological information – point to the Huanan market as the early
epicentre of the COVID-19 pandemic.

## 1\. Introduction

While the origin of the COVID-19 pandemic remains widely debated in the media
and the wider public sphere, scientific contributions have established the
central role played by Wuhan’s Huanan Seafood Wholesale Market (hereafter
“Huanan market”) during the early days of the COVID-19 pandemic, given data
currently available (Worobey, 2021; Worobey et al., 2022a; Jiang and Wang,
2022; Crits-Christoph et al., 2023b). Although this market almost immediately
became the main focus of attention because so many of the very first
identified COVID-19 cases worked or spent time there, the identification of
other early cases without a known exposure risk at the market (Huang et al.,
2020) introduced uncertainty about the role of the Huanan market as the
potential source of the outbreak: to some, these cases appeared to be strong
evidence that the epidemic in Wuhan began elsewhere and was perhaps already
geographically widespread in the city in December 2019 (Cohen, 2020).

An important turning point in the scientific understanding of the pandemic’s
origin occurred in 2021, with the release of a report from a joint World
Health Organization-China COVID-19 origins study, the “WHO-convened global
study of origins of SARS-CoV-2: China part” (hereafter “WHO-China report”)
(World Health Organization, 2021). This report contained maps of residential
locations of early COVID-19 cases in Wuhan, those whose symptoms had begun in
December 2019 (hereafter, “December case locations”; there are no earlier
known onset dates). The study report noted that the residences of these cases
were concentrated in the central districts of Wuhan, where the Huanan market
was also located.111“ _There was a concentration of cases, both laboratory-
confirmed and clinically diagnosed, in the central districts (which include
the Huanan market). The earliest cases were mostly resident in the central
districts of Wuhan, but cases began to appear in all districts of Wuhan in
mid-to late December 2019_ ”; World Health Organization (2021, p. 44).

Starting a few months later, Holmes et al. (2021) and Worobey et al. (2022b,
a) manually extracted coordinates from the maps of the WHO-China report and
conducted spatial analyses. Worobey et al. (2022a) established that the early
cases were geographically centred near the Huanan market to a degree that
would be extremely unlikely unless the outbreak had begun there or at another
location near it. Based not just on these spatial findings but also many other
lines of evidence, Worobey et al. (2022a) concluded that the Huanan market was
the “ _early epicenter_ ” of the COVID-19 pandemic.

Of the geographic findings in Worobey et al. (2022a), the most striking and
consequential involved the early cases epidemiologically unlinked to the
market (i.e., those who had not been to the market and who reported no known
contact with anyone who had). (Hereafter, as in Worobey et al. (2022a), we use
“linked cases” and “unlinked cases” only to refer to ascertained
epidemiological linkage, regardless of whether cases were geographically
associated with the market.) Analysed separately from linked cases, who worked
at or were otherwise knowingly connected to the market, not only was the
geographical centre of the unlinked case locations revealed to be closer to
the market than expected from a null population distribution, but it was also
closer to the Huanan market than the centre of linked case locations (Worobey
et al., 2022a). It had previously been argued that the initial requirement of
a market link for case definition could have biased the early case data
towards the Huanan market. While it is important to take into account
potential ascertainment bias, multiple arguments indicate that it was more
limited than some claimed. First, case definitions were rapidly updated and
did not include a market link (World Health Organization, 2021, Annex E3).
Second, a large fraction of the cases with symptom onset in December were
identified retrospectively, after case definitions had been updated. Third,
the very existence of large numbers of unlinked December cases – about twice
as many as linked cases – directly disproves the claim that the cases with
onset in December 2019 were identified via their link to the market. Fourth,
the spatial structure of these early, unlinked cases provided compelling
evidence for a connection between the Huanan market and the onset of the
outbreak in Wuhan (Worobey et al., 2022a).

Stoyan and Chiu (2024, hereafter, SC2024) contested some of the results
presented by Worobey et al. (2022a, hereafter, W2022). SC2024’s criticisms,
however, stem from  (i) profound misunderstandings of W2022’s study and
objectives,  (ii) a crude technical implementation of a key function, and
(iii) a very narrow scientific perspective.

## 2\. The source of an outbreak is not expected to be at the exact centre of
case locations

### 2.1 A profound misunderstanding of Worobey et al. (2022a)’s claims

While W2022 inferred, based on both geographical and other evidence, that the
market was the “ _epicenter_ ” of the pandemic, they never claimed that the
market was positioned exactly at the centre of the cloud of the December case
locations.

Even in cases where a contagious disease is initially transmitted from a point
source, there are many reasons why that source will not fall exactly at the
centre of case residential locations. Geographic constraints, like the
arrangement of streets or unequal population densities in different directions
from the source, are obvious reasons. Long distance dispersal events as people
move around an urban landscape is another. There is also not a single possible
case map, and therefore not a single centre of case locations: workplaces
could for instance have been mapped instead of residences – there would have
been many more points right at the Huanan market. The case map would also have
been different if the locations of all individuals at a precise time on a
precise day had been mapped. In addition, as time goes by, the distribution of
cases will change, and not necessarily isotropically. Indeed, as expected as
neighbourhood transmission gives way to widespread transmission across a city,
in Wuhan the distribution shifted dramatically over time from the
neighbourhoods around Huanan market towards the centre of the city’s high
density areas, southeast of the market (see Figs. 1D–E in Worobey et al.,
2022a, and our Figure 1). This illustrates that even if an outbreak is
initially clustered around its point source, time can be the biggest enemy to
detecting that signal.

$(a)$ December 2019 cases ($n=155$) |  $(b)$ Unlinked December 2019 cases only ($n=120$)
---|---
|
$(c)$ January 2020 Weibo cases ($n=440$) |  $(d)$ February 2020 Weibo cases ($n=233$)
|
Figure 1: Positions of the modes of case locations and the associated $95\%$
confidence surfaces, inferred by computing the centres of bootstrap
pseudoreplicates, using Worobey et al. (2022a)’s KDE function. $(a)$ all
December case locations from W2022; $(b)$ unlinked cases only; $(c)$ Weibo
cases with confirmation time in January 2020; $(d)$ Weibo cases with
confirmation time in February 2020. The December case locations ($(a)$–$(b)$)
were compiled by W2022 from World Health Organization (2021); the Weibo data
($(c)$–$(d)$) were compiled by Xu et al. (2020).

For these reasons, the rejection of SC2024’s null hypothesis of exact
centrality would not actually tell us much about the role of the Huanan market
as the possible source. It is therefore preposterous to exclude the market on
the basis of a too-stringent definition of “centre”.

### 2.2 Impact on the choice of statistical tests

An important divergence between W2022 and SC2024 relates to the question of
how best to test the hypothesis that a site – identified in advance of spatial
analyses by compelling non-geographic evidence – may have been the place from
which the early outbreak spread. For COVID-19, the Huanan market is the
preeminent such site.

W2022 asked whether the centre of the December case locations was closer to
the Huanan market than other possible “epicenter” locations in the city were.
SC2024 asked whether December case locations were centred on the Huanan
market. As detailed already, there is no reason to expect that the centre of
early case locations will necessarily be exactly centred on the initial
source. SC2024’s null hypothesis is therefore not appropriate, and rejecting
it is not informative. W2022’s tests, on the other hand, compared the
distribution of December case locations to chosen null distributions, and they
characterized the distributions through their proximity to the Huanan market.

W2022 considered two types of null distributions. The first was Wuhan’s
population, using age-adjusted data from worldpop.org. This null distribution
assumes that the disease has had so much time to spread that any signal of its
initial source has vanished, and the distribution of cases matches the
underlying distribution of the population of susceptible hosts. However,
population density data for Wuhan from worldpop.org may not be as accurate as
other available sources,222Daniel A. Walker, 2023-08-23, “ _The worldpop data
are highly discordant with other sources regarding Wuhan population density_
”, https://archive.is/WVSDv. and we therefore use another source in our
analysis here (Peng et al., 2020b, shown in Figure S1).

The second type of null distribution considered by W2022 consisted of self-
reported COVID-19 cases on a social media service (Weibo), with onsets in Jan-
Feb 2020 (from Peng et al., 2020a). This distribution provided a snapshot of
case distributions later in the spread of the disease, but from a point when
clustering of cases may still have existed. In our analysis, we use another
dataset from the same original Weibo source (Xu et al., 2020), for which onset
dates were available, unlike for the Weibo dataset used in W2022 (see Figure
1c; we use the January data as null distribution).

In their first set of tests, W2022 first compared the median distance to
Huanan of sets of locations drawn from a null distribution, to the median
distance to Huanan of the December case locations. Tests presented in W2022’s
main text were based on a null distribution defined using worldpop.org
population density data, and as such, did not include a notion of clustering:
again, these were tests of the scenario outlined above, that the outbreak had
had enough time to spread such that the spatial distribution of cases matched
that of the population of hosts, and therefore lacked clustering. Although
this is not a hypothesis that W2022 thought was likely to be true, in part
because it is refuted by genomic evidence indicating that the Wuhan outbreak
likely began only in late November (Pekar et al., 2022), the idea of a hidden
outbreak that long predated cases at the Huanan market has been proposed by
others (e.g., Cohen, 2020). Despite the limitations of this first test
conducted by W2022, the rejection of this null hypothesis thus provided an
important confirmation that this was not the case. Finally, unlike the
population density data, and contrary to SC2024’s assertions, the second null
distribution considered by W2022, Weibo case data, did encompass case
clustering and human-to-human transmission (see Worobey et al., 2022a, Table
S4, “median distance”).

Case clustering was also encompassed in the second, and more important, set of
tests conducted by W2022. These tests compared the distance from the centre of
the December case locations to the market against the distances to the Huanan
market of single points drawn from a null distribution of plausible starting
points of the outbreak. On Figure 2, this amounts to comparing distance
$d_{DH}$ to the distribution of distances $d_{FH}$. Implicit in this test was
the idea that the single point drawn in each replicate represented the first
human infection in Wuhan, which then acted as the ultimate source of all
subsequent infections, and that the centre of the resulting cluster of cases
would have been near the position of the first infection. (It was assumed, but
not explicitly modelled, that simulating a clustered outbreak around each
first infection, then inferring its centre, would have produced nearly
identical null distributions, just with some jitter around the location of
each first infection.) For both null distributions considered by W2022
(worldpop.org population density, and Weibo cases), the centre of the December
case locations was significantly closer to the Huanan market than were the
putative centres of case locations drawn from each null distribution (see
Worobey et al., 2022a, Table S4, “center-point distance”). As part of the
current study, we reproduced this test and its results with different data
sources for the two null distributions considered by W2022, and using the mode
as centre (see below for a description of the mode). The mode of the December
case locations was significantly closer to the Huanan market, whether we
compared it to first infection locations drawn from Wuhan’s population density
($p=8.9\times 10^{-5}$; $10^{6}$ draws), or from the January 2020 Weibo cases
($p<0.0023$).

We also implemented a version of the test in which we compared the distance
between the market and the centre of the December case locations ($d_{DH}$ in
Figure 2) against the distances between the market and the centres of sets of
the same number of cases ($n=155$) drawn from a null distribution (i.e., the
distribution of distances $d_{CH}$ in Figure 2). We used the mode as centre.
This test is meant to provide a more apples-to-apples comparison than the
previous ones described above, by comparing the positions of putative centres
inferred from the null distribution of cases with the centre of the observed
cases. It is additionally meant to mitigate the clustering issue – the
distribution of the centres of sets of points is much more condensed than the
distribution from which these points are drawn (see Figure S2). We found that
the centre of the December case locations was significantly closer to the
Huanan market than were the centres of randomly sampled cases, whether we used
used sets of ($n=155$) locations from Wuhan’s population density as the null
distribution ($p<5\times 10^{-4}$; $1999$ replicates; Figure S2a) or from the
January 2020 Weibo cases ($p<5\times 10^{-4}$; $1999$ replicates; Figure S2b).

Figure 2: Schematic of the different distances compared in the various tests
of market closeness and centrality.

While W2022 asked whether the centre of the December case locations was closer
to the market than (appropriate) random locations, SC2024 asked a different
question: whether the December case locations were centred on the market.
SC2024’s question did not involve the identification of a null distribution;
instead, SC2024 bootstrapped December case locations. They compared the
distance between the Huanan market and the centre of the December case
locations ($d_{DH}$ in Figure 2), against the distances between centres of
bootstrapped sets of December case locations and the centre of the December
case locations ($d_{DB}$ in Figure 2). Their tests (and in particular, as we
will see below, their method for computing the mode) led to the rejection of
the hypothesis that the market was the centre of the December case locations,
because the market was just outside of the region covered by the bootstrapped
centres. Again, though, the source of the outbreak is not expected to be
exactly at the centre of the December case locations; SC2024’s comparison is
therefore based on a false premise. We nevertheless explore SC2024’s approach
to further highlight its shortcomings.

## 3\. Defining an appropriate centre of case locations and computing its
location

SC2024 discussed different types of centres of a cloud of points. They
asserted that a centre-point, defined by W2022 as the coordinate-wise median,
is “ _a questionable choice_ ”, because of its lack of rotational invariance.
SC2024 proposed two alternative centres, which are rotationally invariant: the
centroid (defined as the coordinate-wise mean) and the mode (defined as the
location of the maximum of a corresponding spatial density function). We will
discuss the relevance of the three types of proposed centres, and will show
that SC2024’s rejection of the Huanan market as centre was due, specifically,
to a crude implementation of a kernel density estimation (KDE).

### 3.1 The impact of the lack of rotation-invariance of the centre-point is
limited

A rotation-invariant centre does not depend on the orientation of the axes,
which, all else being equal, seems to be a desirable property for the centre
of a cloud of points – while natural, the latitude-longitude orientation is
not the only orientation of axes one could use to reference locations. We
tested how the position of the centre-point changes with the orientation of
the axes. We find that, while the computed position of the centre-point does
indeed change slightly, it remains in a circumscribed area, in the vicinity of
the Huanan market (Figure S3). The centre-point may not be strictly speaking
rotation-invariant, but the overall conclusions of W2022 are not affected.

### 3.2 A centroid is affected by extreme values

SC2024 proposed to use a centroid (coordinate-wise means) as the first of two
rotationally invariant alternatives to a centre-point. W2022 had however
chosen a centre-point because a median is more robust to extreme values than a
mean. The December case locations data used by W2022 contain extreme values,
and are not symmetric (Figure S4), so a centre based on medians was preferred
by W2022.333“ _We used medians rather than means for our analyses so as to not
give undue influence to outliers like those that can be seen in fig. S8._ ”,
Worobey et al. (2022a, Supplementary text).

In addition, because their focus was on Wuhan, the clear epicentre city of the
pandemic, W2022’s data did not include ten additional cases mentioned in the
China-WHO report (World Health Organization, 2021, Fig. 23, p.45; annex Fig.
3, p.147). These cases were in seven cities in Hubei province, but outside of
Wuhan, and their features (sequenced or not, type of case confirmation) were
not described. To investigate the impact of minor changes to data on centroids
and centre-points, we digitised the positions of these cities, and added them
to the December case locations data from W2022 (Figure S5). We then repeated
SC2024’s methodology and computed the new positions of the centroid and the
centre-point, and their $95\%$ confidence regions using SC2024’s bootstrap
methodology (Figure S6; $1999$ bootstrap resamples). As expected since it is a
median-based centre, the centre-point and its $95\%$ confidence surface are
barely affected by the inclusion of the seven Hubei cases (Figure S6a–b).
However, because the new cases are distant from the centre of Wuhan (about 150
km for the furthest away), their inclusion changes the position of the
centroid, but also considerably widens the $95\%$ confidence surface (Figure
S6c–d). While the market was outside of the confidence region when only cases
within Wuhan were considered, the market is now clearly inside when the seven
Hubei locations are added ($p=0.2$), and its centrality is not rejected
anymore by SC2024’s (albeit dubious) test.

Worobey et al. (2022a) chose centre-points to avoid just such issues with
extreme values. Indeed, it is not advisable to rely upon a measure of central
tendency that is, like the centroid, so sensitive to minor changes to the
data.

### 3.3 Stoyan and Chiu’s exclusion of the Huanan market is due to the use of
a circular kernel with too high a bandwidth

The mode, i.e. the peak of the underlying spatial density distribution, was
the second alternative centre proposed by SC2024. This suggestion, described
as rotation-invariant (unlike the centre-point), is (like the centre-point but
unlike the centroid) relatively immune to extreme values. As such, in
important ways it combines the main advantages of both the centre-point and
the centroid. The determination of a mode, however, requires the choice of a
kernel function and, specifically, of a bandwidth for kernel density
estimation.

SC2024 used a circular Gaussian kernel and a large value for their bandwidth,
$3000$ m, which had been considered in the preprint version of Worobey et al.
(2022b), in which kernel density estimates (KDE) were computed using ArcGIS
Online. The published peer-reviewed article (Worobey et al., 2022a), however,
employed different software and methods than the preprint (R, function kde
from the ks package (Duong, 2024)), and, importantly, automatic bandwidth
selection, yielding a matrix (via the Hpi function) instead of a scalar
(Figure S7 illustrates the impacts of these choices on the shape of the
spatial density function). Rather than attempting to replicate the methods in
the published Science paper, SC2024 emulated (partially) the Worobey et al.
(2022b) preprint methodology, using a circular kernel (with the bkde2D
function from the KernSmooth package (Wand, 2023), and setting the bandwidth
to $h=3000$ m).

Implementing SC2024’s bootstrap test but with W2022’s more sophisticated KDE
function, we find that the mode lies at the entrance of a parking area at the
Huanan market, and that the market is within the mode’s $95\%$ confidence
region (Figure 1a; $p=0.89$; and see Figure S8 for a detailed view). The
market is thus in fact not rejected as the centre of the distribution of the
December case locations, even using SC2024’s preferred but inappropriately
stringent hypothesis test. A similar non-rejection is obtained with Stoyan and
Chiu (2024)’s KDE function when using a bandwidth automatically determined
using a principled rather than arbitrary approach ($h=866$; Figure S9b;
$p=0.42$).

## 4\. Stoyan and Chiu overlooked key epidemiological data

### 4.1 A limited perspective

A key feature of SC2024’s analysis is its limited perspective. Indeed, the
authors overlooked almost the entirety of the scientific evidence available
about the early period of the pandemic, considering only the December case
locations. Other pieces of evidence are however the reason that the Huanan
market has, since almost the first inklings that a new disease may have
emerged, been the prime suspect as the source of Wuhan’s COVID-19 outbreak.
The link to the Huanan market of early patients with pneumonias of unknown
aetiology cases that would turn out to be COVID-19 helped lead to the
discovery of the new disease, and the market was rapidly closed down (The-nCoV
Outbreak Joint Field Epidemiology Investigation Team and Li, 2020; Yang,
2024); this single market accounted for a substantial proportion of the
earliest-onset COVID-19 patients in what is a very large city (Worobey, 2021);
it had been identified years before the pandemic as a site where viruses with
pandemic potential might jump from animals into humans (Zhang and Holmes,
2020); and it was one of only four markets in Wuhan with sustained sales of
live mammals from intermediate host species known to harbour SARS-related
coronavirus (Xiao et al., 2021). Additionally, although the data were not
public until after Worobey et al.’s study was conceived and published, genetic
traces of these animals were found in environmental samples from the market,
in stalls where SARS-CoV-2 was also detected (Liu et al., 2023; Crits-
Christoph et al., 2023a, b; Débarre, 2024).

### 4.2 Linked vs. unlinked cases

Notably, SC2024 also overlooked the fundamental difference between the linked
and the unlinked December 2019 cases as it relates to spatial analyses. Cases
epidemiologically linked to the market worked there or had shopped there, and
for some were clearly involved in human-to-human transmission chains within
the market (Li et al., 2020). These linked cases are thus expected to be found
around the Huanan market only because people often shop close to where they
live, and often live not too far from where they work. Unlinked cases, on the
other hand, are only expected to be found to reside near the market if local
transmission originated from the market, and only at an early point in the
outbreak, before the virus had spread widely. The observed spatial pattern for
unlinked cases was therefore particularly noteworthy. And yet, even though
information on market links was available in the case dataset used by SC2024,
they did not use it in their analyses. With Stoyan and Chiu (2024)’s
methodology and W2022’s KDE function, we find that analysing only unlinked
cases gives the same results as considering all cases regardless of linkage:
the Huanan market cannot be rejected as being at the very centre of the
unlinked cases (Figure 1b; $p=0.57$).

### 4.3 Missing the tree for the forest

In their paper, SC2024 proffer multiple other landmarks, in the vicinity of
the Huanan market, that could have played a role in the origin of the
pandemic, based on the argument that they are close to one of the centres they
computed from the December case locations. It is of course true that, when
considering the December case locations data in an epistemic vacuum, it is not
possible to rule out alternative sources also in the vicinity of their centre.
Other data, however, help assess the plausibility of the other landmarks
proposed by SC2024. This is for instance the case for a complex called “Wanda
plaza”.444We note that this complex, as well as other landmarks, were
improperly located by SC2024. SC2024 seem to have used Google Maps to extract
coordinates, but did not take into account the offset on maps of China shown
on Google Maps. We present corrected locations (see Figure S10). We also note
that the complex does not just include a shopping centre, but also residences.
The social media check-in data used by SC2024 do not seem to differentiate
between the different uses of the place. SC2024’s post hoc reasoning, in
inferring a centre of the December case locations and then looking for a
landmark near it, completely misses the fact that human-to-human transmission
did happen in relation to the Huanan market (Li et al., 2020). The market did
play an epidemiological role, and is of interest not just because of its
position in Wuhan: it was of interest because of its role in the early days of
COVID-19, before it was known that the centre of the December case locations
would be revealed to be so close to it – as already described above. Landmarks
are to be chosen by taking into account actual potential sources. While the
market is one such potential source, Wanda Plaza has, to our knowledge, never
been mentioned as such until SC2024. Highlighting the existence of other data
linking cases to the Huanan market, and the absence of such evidence for other
landmarks in Wuhan, was actually the point of the “No other location except
the Huanan market clearly epidemiologically linked to early COVID-19 cases”
paragraph in W2022’s supplementary text, that SC2024 seem to have
misunderstood.

SC2024 used social media check-in data to argue that Wanda Plaza was a much
more visited place than the Huanan market. Putting aside the fact that Wanda
Plaza is more than just a shopping centre (see footnote 4), SC2024’s argument
reveals a profound misunderstanding of the use of social media check-in data
by W2022. The Huanan market being much less visited than Wanda Plaza actually
emphasises the extraordinary role played by the Huanan market in the early
days of COVID-19. The market was not a hotspot in Wuhan in terms of visitation
check-ins – it was one of hundreds of possible early sites for case
clustering, and even amongst markets across the city, there were 70 with (up
to hundred times) more check-ins (Worobey et al., 2022a). And yet, the Huanan
market was the hotspot for December cases and in the distribution of the
December case locations.

A new location of the Wuhan Centre for Disease control (WCDC; see Figure S11)
was also proposed as the potential source of the Wuhan outbreak. WCDC was the
focus of one of the earliest speculations about a potential non-natural origin
of SARS-CoV-2 (Xiao and Xiao, 2020), and was identified because of its
proximity to the market (i.e., post-hoc). However, WCDC had moved to its new
location close to the market on December 2nd, 2019 (World Health Organization,
2021, p.122 and Annex D5), leaving little time or no time for its BSL2
laboratory to be operational – and this laboratory did not conduct elaborate
virology experiments like the Wuhan Institute of Virology. While WCDC had been
involved in wildlife sampling, no mammalian samples were moved between the old
and new locations (Holmes, 2024). Finally, as noted in Worobey et al. (2022b,
a), only one WCDC staff member had a positive serology result (World Health
Organization, 2021, Annex D5), to be contrasted with the at least $30$ Huanan
market vendors who were identified as cases. Moreover, this individual was
reported to have been infected by a family member and not at WCDC (World
Health Organization, 2021, Annex D5).

Finally, it is useful to zoom out and to consider the distribution of the
December case locations in Wuhan relative to the locations of the campuses of
the Wuhan Institute of Virology (WIV). This institute has been at the centre
of speculations around a potential lab leak since early 2020. However, Figure
3 shows the clear disconnect between the December case locations and the WIV
campuses.

Figure 3: Zooming out. The two campuses of the Wuhan Institute of Virology
(WIV) are shown.

## 5\. Discussion

Using the same methodology as Stoyan and Chiu (2024, (SC2024)) on the data
extracted by Worobey et al. (2022a, (W2022)), but with more appropriate
parameters and implementation, we arrive at opposite conclusions, and show
that Wuhan’s Huanan market cannot, after all, be rejected as the source of the
COVID-19 pandemic – even with an overly-stringent test.

This result is not surprising given the amount of other evidence pointing to
the Huanan market. In their article, SC2024 dismissed these other pieces of
information, characterising them as “ _some non-statistical argument_ ”. The
question of the origin of SARS-CoV-2, however, requires the consideration of
multiple lines of evidence – and these include other statistical arguments.
Epidemiological links of numerous early cases to the market were identified
early on (Worobey, 2021; Yang, 2024), and our study re-emphasizes the spatial
proximity between the market and December case locations, including cases for
which no link to the market had been identified (Worobey et al., 2022a).

A second main line of evidence involves wildlife sales: animals susceptible to
SARS-CoV-2 were demonstrated to have been in the Huanan market (Xiao et al.,
2021; Crits-Christoph et al., 2023a, b; Liu et al., 2023; Débarre, 2024),
which was the most active of only four markets in Wuhan with consistent sales
of live mammals that are the most plausible intermediate hosts of SARS-CoV-2
(Xiao et al., 2021). In addition, there was a high concentration of SARS-
CoV-2-positive samples in the corner of the market where most live wildlife
was sold (Wu, 2020; Worobey et al., 2022a; Crits-Christoph et al., 2023b), and
genetic traces of SARS-CoV-2 and animals were detected in the same stalls and
in drains directly below and downstream of these stalls (Liu et al., 2023).

A third line of evidence involves SARS-CoV-2’s genetic diversity. There were
two early SARS-CoV-2 lineages, A and B. Sequenced market-linked cases were of
lineage B, but the earliest-known lineage A cases were close to the market (Lu
et al., 2020; Worobey et al., 2022a), and both lineages A and B were detected
in environmental samples at the market (Liu et al., 2023). In addition,
statistical analyses of the SARS-CoV-2 molecular clock strongly reject a
single introduction into humans of a bat virus-like lineage A-haplotype
ancestor from which lineage B evolved within humans (Pekar et al., 2022). The
inconsistency between ancestor inference with molecular clock methods
(rejecting lineage A) vs. with the consideration of related sarbecoviruses
(favouring lineage A), led to the suggestion that there may not have been only
a single introduction of the progenitor virus into humans, but multiple ones
(Pekar et al., 2022) – a scenario that is fully expected at a wildlife market,
but not with a lab origin. Importantly, a market origin is not contingent on
multiple spillovers, and any origin scenario needs to account for the presence
of early lineage A and lineage B both in and near the market.

Finally, temporal arguments also align with a market origin: phylogenetic and
epidemiological analyses suggested that the COVID-19 outbreak in Wuhan began
only shortly before the earliest known cases at the market became symptomatic
(Jijón et al., 2024; Pekar et al., 2022). Consistent with this, serological
testing of stored blood samples, collected from September to December 2019
from more than 32,000 individuals in Wuhan, revealed that not a single one had
neutralizing antibodies against SARS-CoV-2 (Chang et al., 2023), i.e., that
SARS-CoV-2 was not widely circulating in Wuhan in the last quarter of 2019.
And it was not until December 27th–29th 2019 that the first trickle of reports
of COVID-19 cases was recognized by local health authorities as a concern
(Yang, 2024). All but one of the hospitals first reporting pneumonia of
aetiology cases were very near the Huanan market (Yang, 2024).

We agree with SC2024 that there are of course many ways that the spatial
pattern of cases in an emerging epidemic of an infectious respiratory disease
could have deviated from the Huanan-centred pattern observed here. Indeed,
data of individuals self-reporting their COVID-19 infection on social media
(Weibo) to get help, mostly in January and February 2020, indicate that by
that time cases had spread widely across the city, with locations of these
cases largely recapitulating patterns of population density, particularly for
older people (Peng et al., 2020a, and our Figures 1 and S1). By then, the
pattern placing the Huanan market at the centre of the earlier, December
2019-onset cases, had been obscured by subsequent spread (Worobey et al.,
2022a).

Even in cases of non-transmissible diseases, a source may not be located at
the centre of a cloud of case locations, when environmental factors come into
play. For instance, the anthrax cases near the Sverdlovsk military facility in
Spring 1979 were to the south-east of the building, reflecting the direction
of the wind (Meselson et al., 1994).

Conversely, in the classic case of ‘John Snow and the Broad Street pump’,
which established the field of spatial epidemiology, the Broad street pump was
a well-identified source during the 1854 London cholera outbreak, with quite a
clearly central position compared to case fatality locations (Snow, 1855;
Falcone et al., 2020). As with live wildlife markets and COVID-19, water pumps
were pre-identified by John Snow as potential point sources of water-borne
disease. Snow did not go for a pump, of all possible features in the landscape
of the city, at random; he focused on pumps because he had other reasons to
suspect a role of water in spreading cholera. But even then, inappropriate
criteria of the sort employed by SC2024 would exclude the pump as the source.
Indeed, we applied SC2024’s methodology and their KDE function on the Snow
data; the pump was outside of the $95\%$ confidence surfaces (see Figure 4a)
and thus rejected as the centre of that cholera outbreak. Yet we all
understand that the source of the cholera outbreak was not a cobblestone near
the pump, nor a building across the street from the pump, but the contaminated
water flowing from the pump itself. When we used W2022’s KDE on the Snow data
on the other hand, the pump was not rejected as central (Figure 4b).

$(a)$ Mode, SC2024’s KDE, $h=500$ |  $(b)$ Mode, W2022’s KDE
---|---
|
Figure 4: Identification (or not) of the Broad street pump using John Snow’s
data for the 1854 London cholera outbreak, using Stoyan and Chiu (2024)’s
bootstrap method (the background map shows modern London). The panels show
results for different centres and KDE functions. $p$ values: $(a)$ Mode with
SC2024’s KDE and $h=500$, $p=0.001$; $(b)$ Mode with W2022’s KDE $p=0.12$.

Real-life data are always finite. However, while data on the early days of
COVID-19 could have been shared more quickly, in more accessible formats, and
while some data are known to exist but have not been shared yet (Holmes,
2024), the amount of information available is unprecedented for the emergence
of a new disease. Data from various sources and of various types, all indicate
that wildlife sales in the Huanan market played a key role in the onset of the
COVID-19 pandemic.

In this way, the Huanan market hypothesis is fundamentally distinct from lab
leak hypotheses, which were important to consider (Bloom et al., 2021), but
which encompass many unsupported, highly unlikely, and mutually exclusive
scenarios. (For example, if the virus escaped from the Wuhan CDC, then it did
not escape from the Wuhan Institute of Virology; if it came from a mine in
southern China then it did not come from Laos; if this was a fieldwork
accident then it was not an engineered virus, etc.) The Huanan market
hypothesis, in contrast, has given rise to multiple predictions that were
later confirmed as additional evidence arose. These include  (i) that the
SARS-CoV-2 lineage A, found to be geographically associated with the Huanan
market (Lu et al., 2020; Worobey, 2021; Worobey et al., 2022a), would
eventually be found to have been at the market (Liu et al., 2023), and  (ii)
that cases unlinked epidemiologically to the market would be found to reside
around the market (and see Débarre, 2024, for other examples of such
predictions).  Contrary to what SC2024 contended to have rebutted,555Their
article’s title reads “ _Statistics did not prove that the Huanan Seafood
Wholesale Market was the early epicentre of the COVID-19 pandemic_ ”. Worobey
et al. (2022a)’s study did not claim to “prove” that the pandemic started at
the market. Statistics do not offer “proof”, but rather a means to think
clearly about how probable different hypotheses are. Our results show that it
is highly probable that community transmission in Wuhan began in the
neighbourhoods directly surrounding the Huanan market. Other evidence shows
that the Huanan market was deeply unlikely to be at the centre of the Wuhan
outbreak by pure chance (i.e., if the Wuhan outbreak had not begun there;
Worobey et al., 2022a). We should not lose sight of how remarkable it is that
we have such unique insights into the origin of this pandemic, and we would do
well to be mindful of the lesson all this evidence provides: that to prevent
SARS-CoV-3 from emerging, we must reduce opportunities for pathogens with
pandemic potential to be brought into the heart of large human populations via
the live animals that harbour them.

## Acknowledgements

We thank Alex Crits-Christoph, Zach Hensel, Peter Jacobs, Josh Levy, Lorena
Malpica, and Marc Suchard for discussions and/or comments.
This project has been funded in part with federal funds from the National
Institute of Allergy and Infectious Diseases, National Institutes of Health
(NIH), Department of Health and Human Services (contract no. 75N93021C00015 to
M.W.)

## Competing interests

M.W. has received consulting fees from GLG and compensation for expert
testimony from Endurance Specialty Insurance.

## Data and methods

We based our analyses on Stoyan and Chiu (2024)’s R script (provided as
supplementary information) and on Worobey et al. (2022a)’s data and code,
shared on Zenodo (https://zenodo.org/records/6908012;
scripts/construct_KDE_contours.R). The data we used and our R (R Core Team,
2023) scripts are available on Zenodo at
https://doi.org/10.5281/zenodo.10779463. The maps were drawn thanks to the R
leaflet package (Cheng et al., 2023).

### Data sources

Like Stoyan and Chiu (2024), we used case location data shared by Worobey et
al. (2022a), which were manually extracted from figures in the China-WHO Joint
mission report (World Health Organization, 2021) (see Worobey et al. (2022b,
a) for the extraction method). Following a similar methodology, we manually
extracted the seven additional locations shown in annex Fig. 3 of the China-
WHO Joint mission report (World Health Organization, 2021) (there are ten
additional cases, but seven locations, and we do not know which location(s)
had multiple cases; we therefore only placed seven additional points). Because
some maps in the WHO-China report were too blurry, the linked or unlinked
status of some cases is uncertain; the number of identified linked December
case locations in W2022 – $35$ out of the $155$ cases recovered from a map
that contained $164$ cases within Wuhan (World Health Organization, 2021,
Annex E4) – is lower than the reported number of linked cases in the report –
$55$ out of $168$, including one or more Huanan-linked cases outside of Wuhan;
the $168$ number excludes $6$ cases, out of the total of $174$ December 2019
cases, whose exposure history was unknown (World Health Organization, 2021,
Annex E4). W2022 tested the robustness of their results to this (and other)
sources of uncertainty (Worobey et al., 2022a, Supplementary information).

The Snow data were from the dataset provided by Falcone et al. (2020). The
Weibo case data that we used were compiled by Xu et al. (2020).

Population density data were manually digitised from Peng et al. (2020b, Fig.
10) (This imperfect solution was chosen for lack of access to data). We first
increased contrast between consecutive colors using Gimp, and then extracted
the positions of squares with WebPlotDigitizer (Rohatgi, 2024). We
georeferenced the figure using the positions of two landmarks (tips of
islands). We programmatically re-aligned all positions.

### Landmarks

The position of Stoyan and Chiu (2024)’s additional landmarks was clearly
erroneous (as revealed by a comparison of their Figures 1 and 2), so we
repositioned these landmarks using Baidu Maps and Open Street Maps.

### Kernel Density Estimation (KDE)

Stoyan and Chiu (2024) used the bkde2D function from the KernSmooth package
for 2-dimensional kernel smoothing (Wand, 2023), with bandwidth $h=3000$ m.
Stoyan and Chiu (2024) used a circular Gaussian kernel, i.e. the same
bandwidth value in the horizontal and vertical directions. Given that the
choice of a bandwidth can be a bit arbitrary, we used the bw.diggle function
for automatic bandwidth selection (spatstat package (Baddeley and Turner,
2005)), and obtained a value of $866$ m (i.e., yielding less smoothing).

Worobey et al. (2022a), in the peer-reviewed version of their work, used a
different KDE function: the kde function from the ks package (Duong, 2024).
Like W2022 did, we automatically selected a bandwidth matrix using the
package’s Hpi function. The matrix $H$ affects the amount and the orientation
of smoothing. This matrix is not necessarily (and is not, in our case)
diagonal, i.e. it is not necessarily aligned along the coordinate axes.

### Bootstrap analyses

When implementing Stoyan and Chiu (2024)’s test, we estimate $p$ values using
resampled datasets of December case locations. We draw the same number of data
points as in the original datasets, with replacement, $1999$ times. The $p$
value is given by the relative rank of the distance between the original
centre and the landmark of interest, compared to the distances between the
original centre and the centres of the resampled datasets.

## References

  * Baddeley and Turner (2005) Baddeley A, Turner R. 2005. spatstat: An R package for analyzing spatial point patterns. Journal of Statistical Software. 12:1–42.
  * Bloom et al. (2021) Bloom JD, Chan YA, Baric RS, Bjorkman PJ, Cobey S, Deverman BE, Fisman DN, Gupta R, Iwasaki A, Lipsitch M et al. 2021. Investigate the origins of COVID-19. Science. 372:694–694.
  * Chang et al. (2023) Chang L, Zhao L, Xiao Y, Xu T, Chen L, Cai Y, Dong X, Wang C, Xiao X, Ren L et al. 2023. Serosurvey for sars-cov-2 among blood donors in wuhan, china from september to december 2019. Protein & Cell. 14:28–36.
  * Cheng et al. (2023) Cheng J, Schloerke B, Karambelkar B, Xie Y. 2023. leaflet: Create Interactive Web Maps with the JavaScript ’Leaflet’ Library. R package version 2.2.1.
  * Cohen (2020) Cohen J. 2020. Wuhan seafood market may not be source of novel virus spreading globally. Science. https://www.science.org/content/article/wuhan-seafood-market-may-not-be-source-novel-virus-spreading-globally.
  * Crits-Christoph et al. (2023a) Crits-Christoph A, Gangavarapu K, Pekar JE, Moshiri N, Singh R, Levy JI, Goldstein SA, Suchard MA, Popescu S, Robertson DL et al. 2023a. Genetic evidence of susceptible wildlife in SARS-CoV-2 positive samples at the Huanan Wholesale Seafood Market, Wuhan: Analysis and interpretation of data released by the Chinese Center for Disease Control. Technical report. Zenodo. https://zenodo.org/record/7754299.
  * Crits-Christoph et al. (2023b) Crits-Christoph A, Levy JI, Pekar JE, Goldstein SA, Singh R, Hensel Z, Gangavarapu K, Rogers MB, Moshiri N, Garry RF et al. 2023b. Genetic tracing of market wildlife and viruses at the epicenter of the COVID-19 pandemic. Preprint. Genomics. https://www.biorxiv.org/content/10.1101/2023.09.13.557637v1.full.
  * Débarre (2024) Débarre F. 2024. What we can and cannot learn from SARS-CoV-2 and animals in metagenomic samples from the Huanan market. Virus Evolution. 10:vead077. https://academic.oup.com/ve/article/10/1/vead077/7503693.
  * Duong (2024) Duong T. 2024. ks: Kernel Smoothing. R package version 1.14.2.
  * Falcone et al. (2020) Falcone M, Koschinsky J, Vinten-Johansen P, Coleman T, Anselin L. 2020. John Snow & the Cholera Epidemic in Mid-19th Century London: 8 Datasets With Documentation for Use in GeoDa{}. https://geodacenter.github.io/data-and-lab/data/snow_documentation.pdf.
  * Holmes (2024) Holmes EC. 2024. The emergence and evolution of SARS-CoV-2. Annu. Rev.Virol.. p. In Press.
  * Holmes et al. (2021) Holmes EC, Goldstein SA, Rasmussen AL, Robertson DL, Crits-Christoph A, Wertheim JO, Anthony SJ, Barclay WS, Boni MF, Doherty PC et al. 2021\. The origins of SARS-CoV-2: A critical review. Cell. 184:4848–4856. https://www.cell.com/cell/fulltext/S0092-8674(21)00991-0.
  * Huang et al. (2020) Huang C, Wang Y, Li X, Ren L, Zhao J, Hu Y, Zhang L, Fan G, Xu J, Gu X et al. 2020. Clinical features of patients infected with 2019 novel coronavirus in Wuhan, China. The Lancet. 395:497–506. https://www.thelancet.com/journals/lancet/article/PIIS0140-6736(20)30183-5/fulltext.
  * Jiang and Wang (2022) Jiang X, Wang R. 2022. Wildlife trade is likely the source of sars-cov-2. Science. 377:925–926. https://www.science.org/doi/abs/10.1126/science.add8384.
  * Jijón et al. (2024) Jijón S, Czuppon P, Blanquart F, Débarre F. 2024. Using early detection data to estimate the date of emergence of an epidemic outbreak. PLOS Computational Biology. 20:1–19. https://doi.org/10.1371/journal.pcbi.1011934.
  * Li et al. (2020) Li Q, Guan X, Wu P, Wang X, Zhou L, Tong Y, Ren R, Leung KS, Lau EH, Wong JY et al. 2020. Early Transmission Dynamics in Wuhan, China, of Novel Coronavirus–Infected Pneumonia. New England Journal of Medicine. 382:1199–1207. https://www.nejm.org/doi/10.1056/nejmoa2001316.
  * Liu et al. (2023) Liu WJ, Liu P, Lei W, Jia Z, He X, Shi W, Tan Y, Zou S, Wong G, Wang J et al. 2023. Surveillance of SARS-CoV-2 at the Huanan Seafood Market. Nature. https://www.nature.com/articles/s41586-023-06043-2.
  * Lu et al. (2020) Lu R, Zhao X, Li J, Niu P, Yang B, Wu H, Wang W, Song H, Huang B, Zhu N et al. 2020. Genomic characterisation and epidemiology of 2019 novel coronavirus: implications for virus origins and receptor binding. The Lancet. 395:565–574.
  * Meselson et al. (1994) Meselson M, Guillemin J, Hugh-Jones M, Langmuir A, Popova I, Shelokov A, Yampolskaya O. 1994. The Sverdlovsk Anthrax Outbreak of 1979. Science. 266:1202–1208. https://www.science.org/doi/10.1126/science.7973702.
  * Pekar et al. (2022) Pekar JE, Magee A, Parker E, Moshiri N, Izhikevich K, Havens JL, Gangavarapu K, Malpica Serrano LM, Crits-Christoph A, Matteson NL et al. 2022. The molecular epidemiology of multiple zoonotic origins of SARS-CoV-2. Science. 377:960–966. https://www.science.org/doi/10.1126/science.abp8337.
  * Peng et al. (2020a) Peng Z, Wang R, Liu L, Wu H. 2020a. Exploring Urban Spatial Features of COVID-19 Transmission in Wuhan Based on Social Media Data. ISPRS International Journal of Geo-Information. 9:402. https://www.mdpi.com/2220-9964/9/6/402.
  * Peng et al. (2020b) Peng Z, Wang R, Liu L, Wu H. 2020b. Fine-Scale Dasymetric Population Mapping with Mobile Phone and Building Use Data Based on Grid Voronoi Method. ISPRS International Journal of Geo-Information. 9:344. https://www.mdpi.com/2220-9964/9/6/344.
  * R Core Team (2023) R Core Team. 2023. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing. Vienna, Austria.
  * Rohatgi (2024) Rohatgi A. 2024. Webplotdigitizer: Version 4.7.
  * Snow (1855) Snow J. 1855. On the Mode of Communication of Cholera. John Churchill. London. second edition. https://archive.org/stream/b28985266#page/n3/mode/2up.
  * Stoyan and Chiu (2024) Stoyan D, Chiu SN. 2024. Statistics did not prove that the Huanan Seafood Wholesale Market was the early epicentre of the COVID-19 pandemic. Journal of the Royal Statistical Society Series A: Statistics in Society. p. qnad139. https://doi.org/10.1093/jrsssa/qnad139.
  * The-nCoV Outbreak Joint Field Epidemiology Investigation Team and Li (2020) The-nCoV Outbreak Joint Field Epidemiology Investigation Team, Li Q. 2020. An outbreak of NCIP (2019-nCoV) infection in China–Wuhan, Hubei province, 2019–2020. China CDC Weekly. 2:79.
  * Wand (2023) Wand M. 2023. KernSmooth: Functions for Kernel Smoothing Supporting Wand & Jones (1995). R package version 2.23-22.
  * World Health Organization (2021) World Health Organization. 2021. WHO-convened Global Study of Origins of SARS-CoV-2: China Part: Joint WHO-China Study, 14 January-10 February 2021 : Joint Report. WHO. https://www.who.int/publications/i/item/who-convened-global-study-of-origins-of-sars-cov-2-china-part.
  * Worobey (2021) Worobey M. 2021. Dissecting the early COVID-19 cases in Wuhan. Science. 374:1202–1204. https://www.science.org/doi/10.1126/science.abm4454.
  * Worobey et al. (2022a) Worobey M, Levy JI, Serrano LM, Crits-Christoph A, Pekar JE, Goldstein SA, Rasmussen AL, Kraemer MUG, Newman C, Koopmans MPG et al. 2022a. The Huanan Seafood Wholesale Market in Wuhan was the early epicenter of the COVID-19 pandemic. Science. 377:951–959. https://www.science.org/doi/10.1126/science.abp8715.
  * Worobey et al. (2022b) Worobey M, Levy JI, Serrano LMM, Crits-Christoph A, Pekar JE, Goldstein SA, Rasmussen AL, Kraemer MUG, Newman C, Koopmans MPG et al. 2022b. The Huanan market was the epicenter of SARS-CoV-2 emergence. https://zenodo.org/record/6299600.
  * Wu (2020) Wu G. 2020. Report from the Chinese Center for Disease Control and Prevention’s Office of Virus Control and Prevention, Regarding the Results of Environmental Sample Testing for the Novel Coronavirus Epidemic in Wuhan. Technical Report 53. Chinese Center for Disease Control and Prevention. https://archive.ph/xPBFD.
  * Xiao and Xiao (2020) Xiao B, Xiao L. 2020. The possible origins of 2019-nCoV coronavirus. ResearchGate. 6.
  * Xiao et al. (2021) Xiao X, Newman C, Buesching CD, Macdonald DW, Zhou ZM. 2021. Animal sales from Wuhan wet markets immediately prior to the COVID-19 pandemic. Scientific Reports. 11:1–7. https://www.nature.com/articles/s41598-021-91470-2.
  * Xu et al. (2020) Xu Q, Shen Z, Shah N, Cuomo R, Cai M, Brown M, Li J, Mackey T et al. 2020\. Characterizing weibo social media posts from wuhan, china during the early stages of the covid-19 pandemic: qualitative content analysis. JMIR public health and surveillance. 6:e24125.
  * Yang (2024) Yang DL. 2024. Wuhan: How the COVID-19 Outbreak in Wuhan, China Spiralled out of Control. Oxford University Press. New York, NY.
  * Zhang and Holmes (2020) Zhang YZ, Holmes EC. 2020. A Genomic Perspective on the Origin and Emergence of SARS-CoV-2. Cell. 181:223–227. https://www.cell.com/cell/fulltext/S0092-8674(20)30328-7.

## Appendix

Figure S1: Heatmap of population density in Wuhan, manually digitised from Peng et al. (2020b)’s Figure 10.  $(a)$ With population density data |  $(b)$ With Jan Weibo cases
---|---
|
Figure S2: Positions of the modes of sets of $155$ locations, drawn from $(a)$
Wuhan’s population density or $(b)$ from January 2020 Weibo cases; $1999$
resampled sets. Figure S3: Positions of the centre-point of December case
locations when the coordinates axes are rotated. The centre-point with regular
latitude-longitude axes is the centre of rotation, and we tested rotations of
$101$ evenly spaced angles between $0$ and $2\pi$; the resulting centre-points
are shown in lighter colour.  $(a)$ Distribution of longitudes of the December
case locations
---
$(b)$ Distribution of latitudes of the December case locations
Figure S4: Distribution of longitudes and latitudes of the December case locations in the dataset used by Worobey et al. (2022a). The vertical lines correspond to the means (dashed lines) and medians (full lines) of the distributions. Figure S5: All case residential locations, including seven Hubei case locations outside of Wuhan (the city contour is delimited in black). $(a)$ Centre-point, cases inside of Wuhan |  $(b)$ Centre-point, adding the seven Hubei locations
---|---
|
$(c)$ Centroid, cases inside of Wuhan |  $(d)$ Centroid, adding the seven Hubei locations
|
Figure S6: Effect of the addition of new cases on the position of the centre-point and centroid and on their $95\%$ confidence surfaces. The full distribution of cases is shown in Figure S5. The contour delimitates the $95\%$ closest bootstrapped centre locations. $p$ values: $(a)$ centre-point, cases inside Wuhan, $p=0.0085$; $(b)$ centre-point, with Hubei cases, $p=0.012$. $(c)$ centroid, cases inside Wuhan, $p=0.022$; $(d)$ centroid, with Hubei cases, $p=0.2$; $1999$ resamples of the data.  $(a)$ SC2024’s KDE, $h=3000$ m |  $(b)$ W2022’s KDE, bandwidth matrix $H$
---|---
|
$(c)$ SC2024’s KDE, $h=866$ m |
|
Figure S7: The impacts of KDE functions and bandwidth choice ($h$) on kernel density estimation , Wuhan data. SC2024 used $h=3000$ but a coarser grid, $501\times 501$, than ours ($1001\times 1001$), which is why the mode is at the exact same position (but is still very similar; compare the positions of the modes in their Figure 3a to our Figure S7a). In $(b)$, the bandwidth matrix $H$ is automatically determined by the Hpi function; in $(c)$ the bandwidth value $h$ is automatically determined by the bw.diggle function.  $(a)$ Satellite view |  $(b)$ Close-up
---|---
|
Figure S8: Zoom on the position of the mode, computed with W2022’s KDE function, relative to the Huanan market. In panel $(a)$, the red marker between the two sides of the market is the position used throughout our study. Panel $(b)$ further zooms in on the position of the mode (green star in $(a)$); the snapshot is from a picture taken in December 2021, uploaded by user “logan logan” on Google Maps, and available as dynamical image at https://maps.app.goo.gl/q5CmZjrTS3bkSGqW8.  $(a)$ SC2024’s KDE, $h=3000$ |  $(b)$ SC2024’s KDE, $h=866$
---|---
|
Figure S9: Position of the mode and the associated $95\%$ confidence surface,
inferred by computing the centres of bootstrap pseudoreplicates, as a function
of the bandwidth $h$, using SC2024’s KDE function implementation in R. Panel
$(a)$: $h=3000$, which is the bandwidth value used by SC2024 ($p=0.031$).
Panel $(b)$: $h=866$; the bandwidth value was chosen automatically via the
bw.diggle function from the spatstat package (Baddeley and Turner, 2005)
($p=0.42$). The $95\%$ confidence surfaces were determined after $1999$
bootstrap resamples in each panel. We used a finer discretisation grid than
SC2024 ($1001\times 1001$ instead of $501\times 501$).  Figure S10: Positions
of landmarks in Wuhan, using latitude and longitude coordinates shared by
SC2024 (“SC” on the figure) as supplementary information, and showing the
corrected positions. As an example, the Hankou railway station is in reality
to the North-West of the Huanan market (as shown on SC2024’s Figure 2 map). It
appears in the South-West of the market on their Figure 1 and here, because
their coordinates did not take into account the offset on China maps. The
position of Wuhan CDC in the figure is different from SC2024’s Figure 1b,
because we used their latitude and longitude coordinates, which are not
consistent with the UTM coordinates in SC2024’s data (UTM coordinates for
Wuhan CDC seem to have been manually changed in SC2024’s data.) Figure S11:
Equivalent of Figure 1a, showing additional landmarks (with locations
corrected compared to SC2024, see Figure S10).