  Despite the success of adaptive time-stepping in ODE simulation, it has so
far seen few applications for Stochastic Differential Equations (SDEs). To
simulate SDEs adaptively, methods such as the Virtual Brownian Tree (VBT) have
been developed, which can generate Brownian motion (BM) non-chronologically.
However, in most applications, knowing only the values of Brownian motion is
not enough to achieve a high order of convergence; for that, we must compute
time-integrals of BM such as $\int_s^t W_r \, dr$. With the aim of using high
order SDE solvers adaptively, we extend the VBT to generate these integrals of
BM in addition to the Brownian increments. A JAX-based implementation of our
construction is included in the popular Diffrax library
(https://github.com/patrick-kidger/diffrax).
  Since the entire Brownian path produced by VBT is uniquely determined by a
single PRNG seed, previously generated samples need not be stored, which
results in a constant memory footprint and enables experiment repeatability and
strong error estimation. Based on binary search, the VBT's time complexity is
logarithmic in the tolerance parameter $\varepsilon$. Unlike the original VBT
algorithm, which was only precise at some dyadic times, we prove that our
construction exactly matches the joint distribution of the Brownian motion and
its time integrals at any query times, provided they are at least $\varepsilon$
apart.
  We present two applications of adaptive high order solvers enabled by our new
VBT. Using adaptive solvers to simulate a high-volatility CIR model, we achieve
more than twice the convergence order of constant stepping. We apply an
adaptive third order underdamped or kinetic Langevin solver to an MCMC problem,
where our approach outperforms the No U-Turn Sampler, while using only a tenth
of its function evaluations.
