The Bessel process models the local eigenvalue statistics near $0$ of certain large positive definite matrices. In this work, we consider the probability \begin{align*} \mathbb{P}\Big( \mbox{there are no points in the Bessel process on } (0,x_{1})\cup(x_{2},x_{3})\cup\cdots\cup(x_{2g},x_{2g+1}) \Big), \end{align*} where $0<x_{1}<\cdots<x_{2g+1}$ and $g \geq 0$ is any non-negative integer. We obtain asymptotics for this probability as the size of the intervals becomes large, up to and including the oscillations of order $1$. In these asymptotics, the most intricate term is a one-dimensional integral along a linear flow on a $g$-dimensional torus, whose integrand involves ratios of Riemann $θ$-functions associated to a genus $g$ Riemann surface. We simplify this integral in two generic cases: (a) If the flow is ergodic, we compute the leading term in the asymptotics of this integral explicitly using Birkhoff's ergodic theorem. (b) If the linear flow has certain "good Diophantine properties", we obtain improved estimates on the error term in the asymptotics of this integral. In the case when the flow is both ergodic and has "good Diophantine properties" (which is always the case for $g=1$, and "almost always" the case for $g \geq 2$), these results can be combined, yielding particularly precise and explicit large gap asymptotics.