  The production cross-section of $J/\psi$ pairs in proton-proton collisions at
a centre-of-mass energy of $\sqrt{s}=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/\psi$ mesons in the
transverse momentum range $0<p_{\text{T}}<14$ GeV/$c$ and rapidity range
$2.0<y<4.5$. The cross-section of this process is measured to be
16.36$\pm$0.28(stat)$\pm$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 $\Delta y$ between the
two $J/\psi$ mesons. The effective cross-section of double-parton scattering is
measured to be $\sigma_{\text{eff}}=$13.1$\pm$1.8(stat)$\pm$2.3(syst) mb. The
distribution of the azimuthal angle $\phi_{\text{CS}}$ of one of the $J/\psi$
mesons in the Collins-Soper frame and the $p_{\text{T}}$-spectrum of the
$J/\psi$ pairs are also measured for the study of the gluon transverse-momentum
dependent distributions inside protons. The extracted values of
$\langle\cos2\phi_{\text{CS}}\rangle$ and $\langle\cos4\phi_{\text{CS}}\rangle$
are consistent with zero, but the presence of azimuthal asymmetry at a few
percent level is allowed.
