Background: A familiar behaviour of independent and identically distributed (i.i.d.) random variables $X_1, X_2,\ldots X_n$ is concentration: the probability that the sum $X_1+X_2+\ldots X_n$ exceeds $\mathbb{E}(X_1+X_2+\ldots X_n)$ by $n\varepsilon$ decreases as $e^{-n\varepsilon^2}$ (see here). This behaviour is well studied in other settings where similar `independence' is expected, such as for the Gibbs distribution of locally interacting spins.
Surface-order large deviation is a remarkable phenomena that arises in two dimensional spin systems. Consider a $\sqrt{n}\times \sqrt{n}$ square lattice with the set of edges $E$ and the Ising ferromagnetic hamiltonian: $H=-\sum_{(i,j)\in E}\sigma^z_i\sigma^z_j$ ($\sigma^z$ is the `Z' Pauli matrix). The Gibbs distribution of the hamiltonian is $e^{-\beta H}$ (up to normalization). If there is no boundary condition imposed, the expected value of magnetization $M=\sum_i \sigma^z_i$ is $0$, since flipping all the spins in any spin configuration does not change the energy, but changes the value of magnetization from $m$ to $-m$. If the inverse temperature $\beta$ is smaller than a critical value $\beta_c$, the probability (under the Gibbs distribution) that magnetization is larger than $n\varepsilon$ decays exponentially in $n$, similar to the i.i.d. random variables. But for $\beta>\beta_c$, there is a constant $q$, such that the probability that the magnetization is larger than $nq$ decays as $e^{-O(\sqrt{n})}$, which is significantly larger than $e^{-O(n)}$. The source of this behaviour is that if we set the boundary spins to all $1$'s (which has probability $e^{-O(\text{number of boundary spins})}=e^{-O(\sqrt{n})}$), the expected magnetization within the system jumps from $0$ to a value larger than $nq$. Two excellent references on the topic are this and this
Question: Does surface order large deviation persist for $\beta>\beta_c$, when an external magnetic field is present? More precisely, consider the hamiltonian $H=-\sum_{(i,j)\in E}\sigma^z_i\sigma^z_j - h\sum_i \sigma^z_i$ with $h>0$ and the Gibbs state $e^{-\beta H}$ (up to normalization). Let the expected magnetization be $\langle M\rangle$. Is there a constant magnetic field strengh $h$ (constant means there is no dependence on $n$) and a constant $q$, such that the probability that the magnetization exceeds $\langle M\rangle + qn$ is $e^{-O(\sqrt{n})}$?