We present a paradigm for developing arbitrarily high order, linear, unconditionally energy stable numerical algorithms for gradient flow models. We apply the energy quadratization (EQ) technique to reformulate the general gradient flow model into an equivalent gradient flow model with a quadratic free energy and a modified mobility. Given solutions up to $t_n=n Delta t$ with $Delta t$ the time step size, we linearize the EQ-reformulated gradient flow model in $(t_n, t_{n+1}]$ by extrapolation. Then we employ an algebraically stable Runge-Kutta method to discretize the linearized model in $(t_n, t_{n+1}]$. Then we use the Fourier pseudo-spectral method for the spatial discretization to match the order of accuracy in time. The resulting fully discrete scheme is linear, unconditionally energy stable, uniquely solvable, and may reach arbitrarily high order. Furthermore, we present a family of linear schemes based on prediction-correction methods to complement the new linear schemes. Some benchmark numerical examples are given to demonstrate the accuracy and efficiency of the schemes.