  Double parton distributions are the nonperturbative ingredients needed for
computing double parton scattering processes in hadron-hadron collisions. They
describe a variety of correlations between two partons in a hadron and depend
on a large number of variables, including two independent renormalization
scales. This makes it challenging to compute their scale evolution with
satisfactory numerical accuracy while keeping computational costs at a
manageable level. We show that this problem can be solved using interpolation
on Chebyshev grids, extending the methods we previously developed for ordinary
single-parton distributions. Using an implementation of these methods in the
C++ library ChiliPDF, we study for the first time the evolution of double
parton distributions beyond leading order in perturbation theory.
