The valence double parton distribution of the pion is analyzed in the framework of chiral quark models, where in the chiral limit factorization between the longitudinal and transverse degrees of freedom occurs. This feature leads, at the quark-model scale, to a particularly simple distribution of the form D(x1, x2, y) = δ(1−x1 −x2)F (y), where x1,2 are the longitudinal momentum fractions carried the valence quark and antiquark, and y is the separation of their transverse positions. For y = 0 this result complies immediately to the Gaunt-Sterling sum rules. The DGLAP evolution to higher scales is carried out in terms of the Mellin moments. We then explore its role on the longitudinal correlation quantified with the ratio of the double distribution to the product of single distributions, D(x1, x2, y)/D(x1)D(x2). We point out that the ratios of moments x n 1 x m 2 / x n 1 x m 2 are independent of the evolution, providing particularly suitable measures to be tested in the upcoming lattice simulations. The transverse form factor F (y) is obtained in variants of the Nambu-Jona-Lasinio model with the spectral and Pauli-Villars regularizations. Interestingly, with the spectral regularization of the model, the effective cross section for the double parton scattering of pions is exactly equal to the geometric cross section, σ eff = π y 2 and yields about 20 mb.