  During reionization, neutral hydrogen in the intergalactic medium (IGM)
imprints a damping wing absorption feature on the spectrum of high-redshift
quasars. A detection of this signature provides compelling evidence for a
significantly neutral Universe, and enables measurements of the hydrogen
neutral fraction $x_{\rm HI}(z)$ at that epoch. Obtaining reliable quantitative
constraints from this technique, however, is challenging due to stochasticity
induced by the patchy inside-out topology of reionization, degeneracies with
quasar lifetime, and the unknown unabsorbed quasar spectrum close to rest-frame
Ly$\alpha$. We combine a large-volume semi-numerical simulation of reionization
topology with 1D radiative transfer through high-resolution hydrodynamical
simulations of the high-redshift Universe to construct models of quasar
transmission spectra during reionization. Our state-of-the-art approach
captures the distribution of damping wing strengths in biased quasar halos that
should have reionized earlier, as well as the erosion of neutral gas in the
quasar environment caused by its own ionizing radiation. Combining this
detailed model with our new technique for predicting the quasar continuum and
its associated uncertainty, we introduce a Bayesian statistical method to
jointly constrain the neutral fraction of the Universe and the quasar lifetime
from individual quasar spectra. We apply this methodology to the spectra of the
two highest redshift quasars known, ULAS J1120+0641 and ULAS J1342+0928, and
measured volume-averaged neutral fractions $\langle x_{\rm HI}
\rangle(z=7.09)=0.48^{+0.26}_{-0.26}$ and $\langle x_{\rm HI}
\rangle(z=7.54)=0.60^{+0.20}_{-0.23}$ (posterior medians and 68% credible
intervals) when marginalized over quasar lifetimes of $10^3 \leq t_{\rm q} \leq
10^8$ years.
