A bstract The production cross-section of J/ψ pairs in proton-proton collisions at a centre-of-mass energy of $$ \sqrt{s} $$ <mml:math xmlns:mml="http://www.w3.org/1998/Math/MathML"> <mml:msqrt> <mml:mi>s</mml:mi> </mml:msqrt> </mml:math> = 13 TeV is measured using a data sample corresponding to an integrated luminosity of 4 . 2 fb − 1 collected by the LHCb experiment. The measurement is performed with both J/ψ mesons in the transverse momentum range 0 &lt; p T &lt; 14 GeV/ c and rapidity range 2 . 0 &lt; y &lt; 4 . 5. The cross-section of this process is measured to be 16 . 36 ± 0 . 28 (stat) ± 0 . 88 (syst) nb. The contributions from single-parton scattering and double-parton scattering are separated based on the dependence of the cross-section on the absolute rapidity difference ∆ y between the two J/ψ mesons. The effective cross-section of double-parton scattering is measured to be σ eff = 13 . 1 ± 1 . 8 (stat) ± 2 . 3 (syst) mb. The distribution of the azimuthal angle ϕ CS of one of the J/ψ mesons in the Collins-Soper frame and the p T -spectrum of the J/ψ pairs are also measured for the study of the gluon transverse-momentum dependent distributions inside protons. The extracted values of ⟨cos 2 ϕ CS ⟩ and ⟨cos 4 ϕ CS ⟩ are consistent with zero, but the presence of azimuthal asymmetry at a few percent level is allowed.