First-order phase transitions as branch cuts First, I apologize for the long preface to my question. The actual question is in bold below.
In Kardar's Statistical Physics of Fields, on page 10, the author writes, "the free energy of the system is an analytical function in the $(P,T)$ plane, except for a branch cut along the phase boundary." While the above comment is regarding the liquid-gas transition, the same idea is described in the para/ferromagnetic phase transition. In each case, a first order phase transition corresponds to passing through a branch cut, while a second order phase transition occurs exactly at the branch point.
In some ways, describing the phase transition as passing through a branch cut in parameter space makes a lot of sense: namely, you expect your equilibrium observables to vary discontinuously when approaching the branch cut from each side. On the other hand, in complex analysis, branch cuts are usually used when discussing functions which are multivalued (ie, must be defined on multiple Riemann sheets). In this sense, describing the free energy (or other observables) as having branch cuts does not make sense: equilibrium observables are necessarily single-valued, and therefore one could not "push up" the branch cut and explore an additional Riemann sheet.
However, this line of thought does raise an interesting question. Suppose I start in some equilibrium state, and adiabatically vary my parameters so as to explore different equilibrium configurations -- for concreteness, let's say I'm working with an Ising model, so that I can vary magnetic field $h$ and temperature $T$. If I were to start with $h > 0$ and $T < T_c$, and I adiabatically lowered $h$ to pass through zero, Is there some possibility to explore some sort of metastable state with free energy, magnetization, etc, an analytic function of $h$ as $h$ passes through zero? Heuristically, I'm imagining that if we start in the completely ordered spin-up state, for a small field in the "down" direction, each individual spin would rather continue to align upwards with its neighbors instead of flipping downwards towards the field. Of course, much is missing to make this argument rigorous.
 A: First, about your assertion that "equilibrium observables are necessarily single-valued". This is only true for extremal states. At a first-order phase transition, if you consider a non-extremal state, then there always exist
macroscopic observables that do not take deterministic values. Let us, for instance, consider the state obtained by taking the thermodynamic limit with periodic boundary condition. Then the average magnetization (averaged over the whole infinite lattice) would yield either $m^*(\beta)$ or $-m^*(\beta)$, each possibility with probability $1/2$ (here, I used $m^*(\beta)$ to denote the spontaneous magnetization at inverse temperature $\beta$).
In this sense, if one takes into account all possible macroscopic states, then the magnetization can be considered as a multivalued function, taking values $\pm m^*(\beta)$ or even any value in the interval $[-m^*(\beta),m^*(\beta)]$ if you don't look at the magnetization in different samples, but average over samples. This is actually closely related to the fact that the free energy is not differentiable with respect to $h$ at $h=0$ (when $\beta>\beta_{\rm c}$): there are infinitely many "tangents" to its graph at this point. The slope of each of these tangents gives you one possible value of the magnetization. (This can be made very precise.)
This uncertainty is characteristic of a first-order phase transition: fixing the thermodynamic parameters (here, the inverse temperature $\beta>\beta_{\rm c}$ and the magnetic field $h=0$) is not sufficient to determine the macroscopic state (here, whether the system is in the plus or minus phase).

Now, to your main question: "Is there some possibility to explore some sort of metastable state with free energy, magnetization, etc, an analytic function of h as h passes through zero?"
Yes, it is possible to observe (dynamically) a metastable regime. This has been investigated in detail in this paper (see also this related, more recent paper).
The system will indeed remain for some large time in the "plus phase" if you start in the plus phase and set the magnetic field to a small negative value. It will decay to the (true equilibrium) minus phase after a time of the order $\exp(\lambda_{\rm c}/|h|)$, where $\lambda_{\rm c}$ is an explicit constant. The claim is roughly that, for $|h|\ll 1$, at any time of order $\exp(\lambda/|h|)$ with $\lambda<\lambda_c$, the system will be in the "plus phase", while at any time of order $\exp(\lambda/|h|)$ with $\lambda>\lambda_c$, the system will be in the (equilibrium) minus phase.
At time of order $\exp(\lambda/|h|)$ with $\lambda<\lambda_c$, they also derive an asymptotic expansion for the expectation of local observables as a series in $h$, which shows that the "metastable plus phase" can be described as a $C^\infty$ continuation of the equilibrium $+$ phase that exists when $h\geq 0$ to (small) negative $h$.
This analysis is particularly interesting, because it is known that one cannot describe metastability in the Ising model using an analytic continuation of the free energy from positive values of $h$ to negative values of $h$. Indeed, the free energy has an essential singularity at $h=0$ which prevents any such analytic continuation. In this sense, metastability can really be addressed only as a dynamical phenomenon.
(The fact that such a simple description of metastability is possible in the van der Waals-Maxwell theory is a consequence of the mean-field-type assumptions underlying this theory, which prevents spatial segregation of phases).
