  We investigate galaxy sizes at redshift $z\gtrsim 6$ with the cosmological
radiation-magneto-hydrodynamic simulation suite THESAN(-HR). These simulations
simultaneously capture the reionization of the large-scale intergalactic medium
and resolved galaxy properties. The intrinsic size ($r^{\ast}_{1/2}$) of
simulated galaxies increases moderately with stellar mass at $M_{\ast} \lesssim
10^{8}\,{\rm M}_{\odot}$ and decreases fast at larger masses, resulting in a
hump feature at $M_{\ast}\sim 10^{8}\,{\rm M}_{\odot}$ that is insensitive to
redshift. Low-mass galaxies are in the initial phase of size growth and are
better described by a spherical shell model with feedback-driven gas outflows
competing with the cold inflows. In contrast, massive galaxies fit better with
the disk formation model. They generally experience a phase of rapid compaction
and gas depletion, likely driven by internal disk instability rather than
external processes. We identify four compact quenched galaxies in the
$(95.5\,{\rm cMpc})^{3}$ volume of THESAN-1 at $z\simeq 6$, and their quenching
follows reaching a characteristic stellar surface density akin to the massive
compact galaxies at cosmic noon. Compared to observations, we find that the
median UV effective radius ($R^{\rm UV}_{\rm eff}$) of simulated galaxies is at
least three times larger than the observed ones at $M_{\ast}\lesssim
10^{9}\,{\rm M}_{\odot}$ or $M_{\rm UV}\gtrsim -20$ at $6 \lesssim z \lesssim
10$. This inconsistency, related to the hump feature of the intrinsic
size--mass relation, persists across many other cosmological simulations with
different galaxy formation models and demonstrates the potential of using
galaxy morphology to constrain the physics of galaxy formation at high
redshifts.
