  A detection of the 21 cm signal can provide a unique window of opportunity
for uncovering complex astrophysical phenomena at the epoch of reionization and
placing constraints on cosmology at high redshifts, which are usually elusive
to large-scale structure surveys. In this work, we provide a theoretical model
based on a quadratic bias expansion capable of recovering the 21 cm power
spectrum with high accuracy sufficient for upcoming ground-based radio
interferometer experiments. In particular, we develop a hybrid effective field
theory (HEFT) model in redshift space that leverages the accuracy of $N$-body
simulations with the predictive power of analytical bias expansion models, and
test it against the Thesan suite of radiative transfer hydrodynamical
simulations. We make predictions of the 21 cm brightness temperature field at
several distinct redshifts, ranging between $z = 6.5$ and 11, thus probing a
large fraction of the reionization history of the Universe ($x_{\rm HI} = 0.3
\sim 0.9$), and compare our model to the `true' 21 cm brightness in terms of
the correlation coefficient, power spectrum and modeling error. We find
percent-level agreement at large and intermediate scales, $k \lesssim 0.5
h/{\rm Mpc}$, and favorable behavior down to small scales, $k \sim 1 h/{\rm
Mpc}$, outperforming pure perturbation-theory-based models. To put our findings
into context, we show that even in the absence of any foreground contamination
the thermal noise of a futuristic HERA-like experiment is comparable with the
theoretical uncertainty in our model in the allowed `wedge' of observations,
providing further evidence in support of using HEFT-based models to approximate
a range of cosmological observables.
