While the universe becomes more and more homogeneous at large scales , statistical analysis of galaxy catalogs have revealed a fractal structure at small-scales ( \lambda < 100 h ^ { -1 } Mpc ) , with a fractal dimension D = 1.5 - 2 ( Sylos Labini et al 1996 ) . We study the thermodynamics of a self-gravitating system with the theory of critical phenomena and finite-size scaling and show that gravity provides a dynamical mechanism to produce this fractal structure . We develop a field theoretical approach to compute the galaxy distribution , assuming them to be in quasi-isothermal equilibrium . Only a limited , ( although large ) , range of scales is involved , between a short-distance cut-off below which other physics intervene , and a large-distance cut-off , where the thermodynamic equilibrium is not satisfied . The galaxy ensemble can be considered at critical conditions , with large density fluctuations developping at any scale . From the theory of critical phenomena , we derive the two independent critical exponents \nu and \eta and predict the fractal dimension D = 1 / \nu to be either 1.585 or 2 , depending on whether the long-range behaviour is governed by the Ising or the mean field fixed points , respectively . Both set of values are compatible with present observations . In addition , we predict the scaling behaviour of the gravitational potential to be r ^ { - \frac { 1 } { 2 } ( 1 + \eta ) } . That is , r ^ { -0.5 } for mean field or r ^ { -0.519 } for the Ising fixed point . The theory allows to compute the three and higher density correlators without any assumption or Ansatz . We find that the N -points density scales as r _ { 1 } ^ { ( N - 1 ) ( D - 3 ) } , when r _ { 1 } > > r _ { i } , 2 \leq i \leq N . There are no free parameters in this theory .