Mean field theory applied to spin-chains, where the ground state is assumed a product state and the energy is variationally minimized, is one of the tools in my toolkit. The examples I give in this question make me wonder if it behaves inappropriately under the addition of quartic terms and the like. I'm tempted that a resolution could be that one needs to do mean field theory on the effective theory picked by one's favorite method (symmetries a la Landau, path integral via resolutions of the identity or a Hubbard-Stratonovich transformation or the like, etc.), and that mean field theory applied via variational product states is not generally valid except for Hamiltonians with "low-order" terms in some sense. However, my main question is given below in terms of verifying or contradicting the predictions of this sort of mean field theory for a specific model.
I find it very enjoyable that mean field theory techniques, while not getting the critical behavior exactly in every dimension, still yield some important qualitative features of phase transitions. For example, the quantum transverse Ising model,
$$H=-J\left(\sum_{<ij>} \sigma^z_i \sigma^z_j + h/J \sum_{i} \sigma_i^x \right) =-J\left(\sum_{<ij>} \sigma^z_i \sigma^z_j + k \sum_{i} \sigma_i^x \right)$$ has a phase transition as one varies $k$ between a ferromagnetic phase with two nearly degenerate ground states followed by a paramagnetic phase at larger $k$ characterized by a single ground state. The ferromagnetic phase has a net magnetization in the $z$ direction, but the paramagnetic phase does not.
One variant of mean field theory that captures this is to assume the ground state is a translationally invariant product state without entanglement, $\otimes_i |\cos(\theta)\rangle_i$, where $\theta$ captures the directions of the spins between the $x$ and $z$ axis. Then one variationally minimizes with respect to $\theta$. Alternatively, one can characterize the product state by $m_z$ = $\langle \sigma^z\rangle$, as for spin-$1/2$ particles one has $m_z^2+m_y^2+m_z^2 = 1$, so that for this case with spins between the $x$ and $z$ axis, $m_x = \sqrt{1-m_z^2}$.
Then, taking $1d$ for ease (the mean field theory answer will only shift the location of the transition according to the coordination number of lattice), the expected energy per site is
$$\langle H \rangle/L = -J(m_z^2 + k m_x) =-J(m_z^2 + k \sqrt{1-m_z^2})$$
The magnetization $m_z$ that minimizes the energy above is the following, as a function of $k$:
The magnetization continuously goes to zero, with a clear nonanalyticity when it hits zero at $k=2$. While mean field theory gets the critical exponents and location of the phase transition wrong in $1d$, it still shows the qualitative behavior expected above of the $x$ field ultimately beating the $z$ interaction at strong enough fields.
Now, for my question. Let's add a quartic $z$ interaction: $$H=-J\left(\sum_{ijkl {\rm \: in\: a\: line}} \sigma^z_i \sigma^z_j \sigma^z_k \sigma^z_l + \sum_{<ij>} \sigma^z_i \sigma^z_j + k \sum_{i} \sigma_i^x \right)$$
Here the quartic term corresponds to four consecutive spins all in a row or a column.
Using the same mean field theory as above, the problem becomes to minimize (up to numerical prefactors like the coordination number, which will shift the location of the transition)
$$\langle H \rangle/L = -J(m_z^4 + m_z^2 + k \sqrt{1-m_z^2})$$
The resulting minimization problem has a discontinuity! See below:
So mean field theory predicts a discontinuity for this model with an additional quartic term in all dimensions. My question is: For which dimensions, if any, does this model actually show a discontinuity in the magnetic field across the phase transition? (Also, does this model have a name?)
