  The joint moments of the derivatives of the characteristic polynomial of a
random unitary matrix, and also a variant of the characteristic polynomial that
is real on the unit circle, in the large matrix size limit, have been studied
intensively in the past twenty five years, partly in relation to conjectural
connections to the Riemann zeta-function and Hardy's function. We completely
settle the most general version of the problem of convergence of these joint
moments, after they are suitably rescaled, for an arbitrary number of
derivatives and with arbitrary positive real exponents. Our approach relies on
a hidden, higher-order exchangeable structure, that of an exchangeable array.
Using these probabilistic techniques, we then give a combinatorial formula for
the leading order coefficient in the asymptotics of the joint moments, when the
power on the characteristic polynomial itself is a positive real number and the
exponents of the derivatives are integers, in terms of a finite number of
finite-dimensional integrals which are explicitly computable. Finally, we
develop a method, based on a class of Hankel determinants shifted by
partitions, that allows us to give an exact representation of all these joint
moments, for finite matrix size, in terms of derivatives of Painlev\'e V
transcendents, and then for the leading order coefficient in the large-matrix
limit in terms of derivatives of the $\sigma$-Painlev\'e III' equation. This
gives an efficient way to compute all these moments and coefficients
explicitly. Our methods can be used to obtain analogous results for a number of
other models sharing the same features.
