An introduction to Landau theory can be found here.
Based on Ginzburg–Landau theory, the free energy of a ferroelectric material, in the absence of an electric field and applied stress may be written as a Taylor expansion in terms of the order parameter, P.  If a sixth order expansion is used (i.e. 8th order and higher terms truncated), the free energy is given by:
where P, P, and P are the components of the polarization vector in the x, y, and z directions respectively, and the coefficients,  must be consistent with the crystal symmetry.  To investigate domain formation and other phenomena in ferroelectrics, these equations are often used in the context of a phase field model.  Typically, this involves adding a gradient term, an electrostatic term and an elastic term to the free energy.  The equations are then discretized onto a grid using the finite difference method or finite element method and solved subject to the constraints of Gauss's law and Linear elasticity.
In all known ferroelectrics,  and .  These coefficients may be obtained experimentally or from ab-initio simulations.  For ferroelectrics with a first order phase transition, ,  whereas  for a second order phase transition.
The spontaneous polarization, P of a ferroelectric for a cubic to tetragonal phase transition may be obtained by considering the 1D expression of the free energy which is:
This free energy has the shape of a double well potential with two free energy minima at , the spontaneous polarization.  We find the derivative of the free energy, and set it equal to zero in order to solve for  :
Since the P = 0 solution of this equation rather corresponds to a free energy maxima in the ferroelectric phase, the desired solutions for P correspond to setting the remaining factor to zero:
whose solution is:
and eliminating solutions which take the square root of a negative number (for either the first or second order phase transitions) gives:
If , the solution for the spontaneous polarization reduces to:
The hysteresis loop (P versus E) may be obtained from the free energy expansion by including the term -E P corresponding to the energy due to an external electric field E interacting with the polarization P, as follows:
We find the stable polarization values of P under the influence of the external field, now denoted as P, again by setting the derivative of the energy with respect to P to zero:
Plotting E (on the X axis) as a function of P (but on the Y axis) gives an S-shaped curve which is multi-valued in  P for some values of E.  The central part of the S corresponds to a free energy local maximum (since  ).  Elimination of this region, and connection of the top and bottom portions of the S curve by vertical lines at the discontinuities gives the hysteresis loop of internal polarization due to an external electric field.