  An important characteristic of cosmic hydrogen reionization is the growth of
ionized gas bubbles surrounding early luminous objects. Ionized bubble sizes
are beginning to be probed using Lyman-$\alpha$ emission from high-redshift
galaxies, and will also be probed by upcoming 21-cm maps. We present results
from a study of bubble sizes using the state-of-the-art THESAN
radiation-hydrodynamics simulation suite, which self-consistently models
radiation transport and realistic galaxy formation. We employ the mean-free
path method, and track the evolution of the effective ionized bubble size at
each point ($R_{\rm eff}$) throughout the Epoch of Reionization. We show there
is a slow growth period for regions ionized early, but a rapid "flash
ionization" process for regions ionized later as they immediately enter a
large, pre-existing bubble. We also find that bright sources are preferentially
in larger bubbles, and find consistency with recent observational constraints
at $z \gtrsim 9$, but tension with idealized Lyman-$\alpha$ damping-wing models
at $z \approx 7$. We find that high overdensity regions have larger
characteristic bubble sizes, but the correlation decreases as reionization
progresses, likely due to runaway formation of large percolated bubbles.
Finally, we compare the redshift at which a region transitions from neutral to
ionized ($z_{\rm reion}$) with the time it takes to reach a given bubble size
and conclude that $z_{\rm reion}$ is a reasonable local probe of small-scale
bubble size statistics ($R_\text{eff} \lesssim 1\,\rm{cMpc}$). However, for
larger bubbles, the correspondence between $z_{\rm reion}$ and size statistics
weakens due to the time delay between the onset of reionization and the
expansion of large bubbles, particularly at high redshifts.
