The Potts model has been widely explored in the literature for the last few decades. 

While many analytical and numerical results concern with the traditional two site interaction model in various geometries and dimensions, little 

is yet known about models where more than two spins simultaneously interact. We consider a Ferromagnetic Four site interaction Potts model on the Square lattice, 

where the four spins reside in the corners of an elementary square.

Each spin can take an integer value $1,2,...,q$.

We write the partition function as a sum over clusters consisting 

of monochromatic faces.

When the number of faces becomes large, tracing out spin configurations is equivalent 

to enumerating large lattice animals. 

This, together with the assumption that typically, in large animals the number of faces and sites is (to leading order) equal,

implies that systems with $q\leq 4$ 

and $q>4$ exhibit a second and first order phase transitions, respectively.

However, higher order terms can make the borderline $q=4$ systems fall into the first order regime.

We find ${1}/{\log q}$ to be an upper bound on $T_c$, the exact critical point.

Using a low temperature expansion, we show that ${1}/{\theta\log q}$, where $\theta>1$ is a $q$ dependent geometrical term, is an improved upper bound on $T_c$.

Moreover, since large animals uniquely control long range order, we expect that $T_c=1/\theta\log q$.

This expression is used to estimate the finite correlation length in the first order transition case.

These results can be extended to other lattices.

Our analytical predictions are confirmed numerically by an extensive study

of the four site interaction model using the Wang-Landau entropic sampling method for $q=3,4,5$.

In particular, the $q=4$ model shows an ambiguous finite size pseudo-critical behaviour.