Free energy functions are analytic or non-analytic in phase transitions?

Landau free energy is just an approximation to the real free energy in the thermodynamic limit. For that reason, Landau free energy can be analytic, while the real one is not. Let me show how the approximation works.

As you may know, the Landau free energy is defined in the following way, assuming the Ising model:

$$Z\left(h,T\right)=\sum_{\left\{ s_{i}\right\} }\exp\left(-\beta H\left(\left\{ s_{i}\right\} \right)\right)=\sum_{m}\exp\left(-\beta F_{L}\left(m,h,T\right)\right)$$

where $Z$ is the partition function, $\{s_i\}$ stands for all possible spin configurations, and the sum in $m$ stands for every possible magnetization.

In Landau theory, the critical point the point $m^*$ such that $F_L(m^*,h,T)$ is minimum. Then, note that $ \exp\left(-\beta F_{L}\left(m^*,h,T\right)\right) $ is a maximum because of the minus sign. Then, we write

$$\log Z= \log \left [\sum_{m}\exp\left(-\beta F_{L}\left(m,h,T\right)\right) \right ] \geq \log \left [\sum_{m}\exp\left(-\beta F_{L}\left(m^*,h,T\right)\right) \right ],$$

In addition to this inequality, we can get an upper bound to this expression. The sum includes many different values for the magnetization, from -1 to +1; you can convince yourself (this is the hardest part of the demonstration) that if we replace the sum by a product of the $N$ at the minimum configuration, this quantity will be bigger than the original one:

$$\log Z \leq \log \left [N \exp\left(-\beta F_{L}\left(m^*,h,T\right)\right) \right ]$$

Now we almost have it. The quantity is bounded,

$$\log\left[N\exp\left(-\beta NF_{L}\left(m^{*},h,T\right)\right)\right]\geq\log Z\geq\log\left[\exp\left(-\beta F_{L}\left(m^{*},h,T\right)\right)\right]$$

$$\log N-\beta F_{L}\left(m^{*},h,T\right)\geq \log Z \geq-\beta F_{L}\left(m^{*},h,T\right).$$

Now we use the definition of the real, free energy, $F=-kT\log Z$. Multiplying by $1/\beta$ the expression above, we get that the real, non-analytic free energy, is bounded by the Landau free energy at the minimum:

$$k T\log N-F_{L}\left(m^{*},h,T\right)\geq -F\left(h,T\right)\geq-\beta F_{L}\left(m^{*},h,T\right).$$

Then, we change into intensive variables $F_L=Nf_L$ and divide by the number of spins $N$, to get:

$$\frac{\log N}{N}-f_{L}\left(m^{*},h,T\right)\geq f\left(h,T\right)\geq- f_{L}\left(m^{*},h,T\right)$$.

Notice that once we do the thermodynamic limit, the term $\log(N)/N \rightarrow 0$ and then we have that $f_{L}\left(m^{*},h,T\right) = f\left(m^{*},h,T\right)$, so in the thermodynamic limit, the analytic Landau free energy is the same as the real one. However, real system have real not an infinite number of spins, meaning that this is only an approximation. In experiments, there is an huge number of spins, and this is why Landau theory works very well, but if you work with little $N$ the difference between the two is noticeable.

About your second question, I think you are confusing things a bit: the transition is second order because the magnetization, which is the first derivative of the free energy with respect to the external magnetic field, is continuous at the critical point. However, susceptibility, which is the second derivative, is not continuous, meaning it is a second order phase transition.

Maybe what confuses you is that Landau theory gives you a piecewise defined function for magnetization. However, you precisely compute $T_c$ by the constraint that magnetization has to be continuous.

Edit: the non-analyticity of the free energy, if I am not mistaken, is respect to the magnetization -but not with respect to the external field $h$.

The Wikipedia page is being sloppy. They mean that the free energy density is an analytic function of the mean-field order parameter, whereas at a thermal phase transition the free energy density is a non-analytic function of the temperature (or for a zero-temperature phase transition, of the external parameter being tuned across the transition).

The order parameter is a derivative of the free energy density with respect to a symmetry-breaking external parameter. In the paramagnet-ferromagnet transition, the external parameter is the applied field $h$ that enters in the Hamiltonian as a symmetry-breaking term $-h\ M$, in which case the order parameter $M = T\ \partial f / \partial h$ is indeed a derivative of the free energy density.

I'll provide one possible answer to your question. I use the language of the Ising ferromagnet/paramagnet phase transition, but this is just for convenience as the same holds much more generally. Assume that the system is below the Curie temperature. Then, the free energy is an analytic function of the magnetic field both when $h>0$ and when $h<0$. At $h=0$, the free energy is not even differentiable, so it is obviously not analytic. This is not the relevant question.

Namely, the real point is whether the free energy can be analytically continued from $h<0$ to $h>0$. It can be shown that the free energy possesses directional derivatives of all orders at $h=0$, that is, $$ \frac{{\rm d}^n}{{\rm d}_-h^n} f(\beta,h)\vert_{h=0} $$ exists for all $n$, where ${\rm d}/{\rm d}_-h$ means the left-derivative. Okay, so the question reduces to whether the Taylor series converges in a small disk around $0$. This is actually not the case, at least for finite-range models at low enough temperatures: there is no analytic continuation of the free energy beyond the transition point. This result was originally proved by Isakov in this paper. His analysis was restricted to the Ising model, but was more recently extended to a large class of 2-phases models; see the discussion in this paper.

One of the conceptually most interesting consequences of Isakov's result is that the usual description of metastable states as analytic continuation of the free energy across the transition point is incorrect, at least for finite-range interactions (and probably short-range ones too).