We construct a model for double parton distribution functions (dPDFs) based on the notion of self-similarity, pursued earlier for small x physics at HERA. The most general form of dPDFs contains total 13 parameters to be fitted from data of proton–proton collision at LHC. It is shown that the constructed dPDF does not factorize into two single PDFs in conformity with QCD expectation, and it satisfies the condition that at the kinematic boundary x1+x2 = 1 (where x1 and x2 are the longitudinal fractional momenta of two partons), the dPDF vanishes.