  In recent years, a number of methods have been developed to infer complex
demographic histories, especially historical population size changes, from
genomic sequence data. Coalescent Hidden Markov Models have proven to be
particularly useful for this type of inference. Due to the Markovian structure
of these models, an essential building block is the joint distribution of local
genealogical trees, or statistics of these genealogies, at two neighboring loci
in populations of variable size. Here, we present a novel method to compute the
marginal and the joint distribution of the total length of the genealogical
trees at two loci separated by at most one recombination event for samples of
arbitrary size. To our knowledge, no method to compute these distributions has
been presented in the literature to date. We show that they can be obtained
from the solution of certain hyperbolic systems of partial differential
equations. We present a numerical algorithm, based on the method of
characteristics, that can be used to efficiently and accurately solve these
systems and compute the marginal and the joint distributions. We demonstrate
its utility to study the properties of the joint distribution. Our flexible
method can be straightforwardly extended to handle an arbitrary fixed number of
recombination events, to include the distributions of other statistics of the
genealogies as well, and can also be applied in structured populations.
