  We introduce the THESAN project, a suite of large volume (L = 95.5 cMpc)
radiation-magnetohydrodynamic simulations that simultaneously model the
large-scale statistical properties of the IGM during reionization and the
resolved characteristics of the galaxies responsible for it. The flagship
simulation has dark matter (DM) and baryonic mass resolutions of $3.1 \times
10^6$ $M_\odot$ and $5.8 \times 10^5 $ $M_\odot$, respectively. The
gravitational forces are softened on scales of 2.2 ckpc with the smallest cell
sizes reaching 10 pc at z=5.5, enabling predictions down to the atomic cooling
limit. The simulations use an efficient radiation hydrodynamics solver
(AREPO-RT) that precisely captures the interaction between ionizing photons and
gas, coupled to well-tested galaxy formation (IllustrisTNG) and dust models to
accurately predict the properties of galaxies. Through a complementary set of
medium resolution simulations we investigate the changes to reionization
introduced by different assumptions for ionizing escape fractions, varying DM
models, and numerical convergence. The simulations produce realistic
reionization histories that match the observed evolution of the global neutral
hydrogen fraction and electron scattering optical depth to reionization. They
also match a wealth of high-redshift observationally inferred data, including
the stellar-to-halo-mass relation, stellar mass function, star formation rate
density, and the mass-metallicity relation, despite the galaxy formation model
being mainly calibrated at z=0. We demonstrate that different reionization
models give rise to varied bubble size distributions that imprint unique
signatures on the 21cm emission, especially on the slope of the power spectrum
at large spatial scales, enabling current and upcoming 21cm experiments to
accurately characterise the sources that dominate the ionizing photon budget.
[abridged]
