  We present a systematical approach to developing arbitrarily high order,
unconditionally energy stable numerical schemes for thermodynamically
consistent gradient flow models that satisfy energy dissipation laws. Utilizing
the energy quadratization (EQ) method, We formulate the gradient flow model
into an equivalent form with a corresponding quadratic free energy functional.
Based on the equivalent form with a quadratic energy, we propose two classes of
energy stable numerical approximations. In the first approach, we use a
prediction-correction strategy to improve the accuracy of linear numerical
schemes. In the second approach, we adopt the Gaussian collocation method to
discretize the equivalent form with a quadratic energy, arriving at an
arbitrarily high-order scheme for gradient flow models. Schemes derived using
both approaches are proved rigorously to be unconditionally energy stable. The
proposed schemes are then implemented in four gradient flow models numerically
to demonstrate their accuracy and effectiveness. Detailed numerical comparisons
among these schemes are carried out as well. These numerical strategies are
rather general so that they can be readily generalized to solve any
thermodynamically consistent PDE models.
